PEMODELAN CRIME RATE DI JAWA TIMUR DENGAN GEOGRAPHICALLY WEIGHTED REGRESSION (GWR)
Maharani Rizki F
email: maharani22001@mail.unpad.ac.id
NPM: 140610220005
department: Program Studi Statistika, FMIPA Universitas Padjadjaran
dosen pengampu: I Gede Nyoman Mindra Jaya, M.Si., Ph.D.
GEOGRAPHICALLY WEIGHTED REGRESSION
GWR adalah teknis non-stationer yang memodelkan hubungan yang bervariasi secara spasial. Dibandingkan dengan regresi dasar (global), koefisien dalam GWR berupa fungsi dari lokasi spasial. Fotheringham dkk. (1998, 2002) memberikan bentuk model umum dari GWR sebagai berikut :
\[ y_i = \beta_{i0} + \sum_{k=1}^{m} \beta_{ik} x_{ik} + \epsilon_i \]
Dengan:
- \(y_i\): variabel dependen di lokasi ke-\(i\)
- \(x_{ik}\): variabel independen ke-\(k\) di lokasi ke-\(i\)
- \(\beta_{ik}\): koefisien regresi lokal untuk variabel independen ke-\(k\) di lokasi ke-\(i\)
- \(\epsilon_i\): kesalahan acak di lokasi ke-\(i\)
GWR memungkinkan koefisien berubah secara kontinu di seluruh wilayah studi, sehingga menghasilkan estimasi lokal yang berbeda-beda pada setiap lokasi. Dengan kalibrasi di setiap titik, GWR menghitung hubungan lokal dengan mempertimbangkan “area pengaruh” di sekitar titik tersebut—pengamatan yang lebih dekat memiliki pengaruh lebih besar dibandingkan pengamatan yang jauh (Fotheringham et al., 1998). Estimasi koefisien diperoleh menggunakan metode weighted least squares (kuadrat terkecil berbobot), dengan rumus:
\[ \hat{\beta}_i = (X^T W_i X)^{-1} X^T W_i y \]
Di mana:
- \(X\): matriks variabel independen dengan satu kolom intercept
- \(y\): vektor variabel dependen
- \(\hat{\beta}_i\): vektor koefisien regresi lokal di lokasi \(i\)
- \(W_i\): matriks diagonal yang menunjukkan pembobotan geografis di lokasi \(i\)
Langkah-langkah Analisis GWR
Langkah-langkah analisis dalam penelitian ini adalah sebagai berikut : 1. Memetakan penyebaran Crime rate di Jawa timur 2. Mengidentifikasi pola hubungan antar variabel dengan regresi linear (OLS) 3. Menetapkan matriks pembobot spasial (W) 4. Melakukan diagnosa statistik untuk memeriksa apakah model memiliki efek spasial atau heterogenitas spasial yang diinginkan 5. Memeriksa apakah ada efek spasial dengan uji Moran’s I 6. Melakukan pemodelan GWR 7. Melakukan pemilihan bandwidth optimum dengan Cross Validation (CV) minimum 8. Membentuk matriks pembobot (W) untuk setiap lokasi pengamatan 9. Melakukan estimasi parameter model GWR 10. Menarik Kesimpulan
Analisis menggunakan R
library
library(fields)## Loading required package: spam## Spam version 2.10-0 (2023-10-23) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: viridisLite##
## Try help(fields) to get started.library(ggplot2)
library(viridisLite)
library(gridExtra)
library(knitr)
library(gt)
library(terra) ## terra 1.7.78##
## Attaching package: 'terra'## The following object is masked from 'package:knitr':
##
## spin## The following object is masked from 'package:fields':
##
## describelibrary(geodata)##
## Attaching package: 'geodata'## The following object is masked from 'package:fields':
##
## worldlibrary(leaflet)##
## Attaching package: 'leaflet'## The following object is masked from 'package:fields':
##
## addLegendlibrary(maps)
library(raster)## Loading required package: splibrary(stars)## Loading required package: abind## Loading required package: sf## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUElibrary(gtools)
library(sp)
library(spdep) ## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(splancs)##
## Spatial Point Pattern Analysis Code in S-Plus
##
## Version 2 - Spatial and Space-Time analysis##
## Attaching package: 'splancs'## The following object is masked from 'package:raster':
##
## zoom## The following objects are masked from 'package:terra':
##
## as.points, is.points, zoomlibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ tidyr::extract() masks raster::extract(), terra::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::map() masks maps::map()
## ✖ dplyr::select() masks raster::select()
## ✖ dplyr::tribble() masks tidyr::tribble(), tibble::tribble(), splancs::tribble()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsrequire(mgcv)## Loading required package: mgcv
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## The following object is masked from 'package:raster':
##
## getData
##
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
##
## Attaching package: 'mgcv'
##
## The following object is masked from 'package:gtools':
##
## scatrequire(dplyr)
require(DAAG)## Loading required package: DAAG
##
## Attaching package: 'DAAG'
##
## The following object is masked from 'package:maps':
##
## ozonerequire(reshape2)## Loading required package: reshape2
##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smithslibrary(gstat)
library(oce) ## Loading required package: gsw
##
## Attaching package: 'oce'
##
## The following object is masked from 'package:terra':
##
## rescalelibrary(chirps)
library(sf)Memetakan penyebaran Crime rate di Jawa timur
PLOT JAWA TIMUR
variabel Y disini adalah variabel yg ditransformasikan menggunakan transformasi logaritma
setwd("D:/Spatial/")
data <- read.csv("crime jatim.csv", header = TRUE)
Indo_Kec<-readRDS('gadm36_IDN_2_sp.rds')
Jatim<-readRDS('gadm36_IDN_2_sp.rds')
Jatim <-Jatim[Indo_Kec$NAME_1 == "Jawa Timur",]
IDs <- Jatim$NAME_2
coords<-coordinates(Jatim)
Jatim$Y <- log(data$risiko.penduduk.tindak.pidana)
Jatim$X1 <- data$TPT.2022
Jatim$X2 <- data$Pengeluaran.perkapita.2022
Jatim$X3 <- data$pertumbuhan.ekonomi.2022
Jatim$X4 <- data$persentase.kepadatan.2022
Jatim$X5 <- data$belum.tidak.sekolah.2022
Jatim$X6 <- data$penduduk.miskin.2022
Jatim$Kab.Kota <- data$Kab.Kota
# Ubah ke sf object (jika belum diubah)
Jatim_sf <- st_as_sf(Jatim)
# Fungsi untuk membuat peta untuk setiap variabel
plot_variable <- function(data, variable, title) {
ggplot(data) +
geom_sf(aes_string(fill = variable), color = "black") +
scale_fill_gradient(low = "yellow", high = "red", na.value = "white") +
theme_minimal() +
labs(title = title, fill = variable) +
theme(legend.position = "right",
axis.text.x = element_blank(), # Remove x-axis labels
axis.text.y = element_blank()) # Remove y-axis labels
}
# Plot sebaran risiko penduduk terhadap tindak pidana (variabel Y)
plot_Y <- plot_variable(Jatim_sf, "Y", "Risiko Penduduk Tindak Pidana")## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# Plot sebaran variabel independen
plot_X1 <- plot_variable(Jatim_sf, "X1", "Tingkat Pengangguran Terbuka (TPT) 2022")
plot_X2 <- plot_variable(Jatim_sf, "X2", "Pengeluaran Perkapita 2022")
plot_X3 <- plot_variable(Jatim_sf, "X3", "Pertumbuhan Ekonomi 2022")
plot_X4 <- plot_variable(Jatim_sf, "X4", "Persentase Kepadatan Penduduk 2022")
plot_X5 <- plot_variable(Jatim_sf, "X5", "Belum/Tidak Sekolah 2022")
plot_X6 <- plot_variable(Jatim_sf, "X6", "Persentase Penduduk Miskin 2022")Risiko Penduduk Tindak Pidana (Y)
plot_Y 
VARIABEL INDEPENDEN (X1-X6)
library(gridExtra)
grid.arrange(plot_X1, plot_X2, plot_X3, plot_X4, plot_X5, plot_X6, ncol = 2)
Mengidentifikasi pola hubungan antar variabel dengan regresi linear (OLS)
Menentukan variabel
y <- Jatim$Y
x1 <- Jatim$X1
x2 <- Jatim$X2
x3 <- Jatim$X3
x4 <- Jatim$X4
x5 <- Jatim$X5
x6 <- Jatim$X6REGRESI LINEAR MODEL OLS
#Linear model OLS
crimelm <- lm(y ~ (x1+x2+x3+x4+x5+x6))
summary(crimelm)##
## Call:
## lm(formula = y ~ (x1 + x2 + x3 + x4 + x5 + x6))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.60742 -0.19711 0.02819 0.13504 0.57466
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.025e+00 4.318e-01 -2.373 0.02405 *
## x1 -2.981e-02 3.674e-02 -0.811 0.42337
## x2 1.269e-04 3.289e-05 3.859 0.00054 ***
## x3 1.291e-03 2.171e-02 0.059 0.95296
## x4 4.262e-01 5.926e-02 7.192 4.36e-08 ***
## x5 2.866e-02 2.135e-02 1.343 0.18910
## x6 -4.103e-03 1.703e-03 -2.410 0.02209 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2763 on 31 degrees of freedom
## Multiple R-squared: 0.7954, Adjusted R-squared: 0.7558
## F-statistic: 20.09 on 6 and 31 DF, p-value: 1.963e-09Melakukan diagnosa statistik untuk memeriksa apakah model memiliki efek spasial atau heterogenitas spasial
#Model diagnostics
par(mfrow=c(2,2))
plot(crimelm) 
UJI MULTIKOLINEARITAS
Pemeriksaan multikolinieritas ditunjukkan dengan nilai VIF. Adanya hubungan antar variabel prediktor pada model regresi linier berganda ditunjukkan ketika nilai VIFl >10. Nilai VIF pada setiap variabel prediktor kurang dari 10 sehingga dapat disimpulkan bahwa tidak terjadi multikolinieritas antar variabel prediktor
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:DAAG':
##
## vif## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## some## The following object is masked from 'package:gtools':
##
## logitvif(crimelm)## x1 x2 x3 x4 x5 x6
## 2.047876 2.614428 1.348171 4.582783 3.196681 6.464971library(corrplot)## corrplot 0.95 loadedcordata = data[,3:11]
cordata = cor(cordata)
corrplot(cordata,method = 'number')
UJI HETEROKEDASTISITAS
Hipotesis
H0 : Homokedastisitas
H1 : Heterokedastisitas
Kriteria Uji
Tolak H0 jika p-value < alpha. (signifikan)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:terra':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(crimelm)##
## studentized Breusch-Pagan test
##
## data: crimelm
## BP = 6.0998, df = 6, p-value = 0.4121UJI NORMALITAS
Hipotesis
H0 : Data berdistribusi normal
H1 : Data tidak berdistribusi normal
Kriteria Uji
Tolak H0 jika p-value < 0.05
resids <- residuals(crimelm)
shapiro.test(resids)##
## Shapiro-Wilk normality test
##
## data: resids
## W = 0.98452, p-value = 0.8671UJI NON-AUTOKORELASI
Hipotesis
H0 : Data non autokorelasi
H1 : Data terdapat autokorelasi
Kriteria Uji
Tolak H0 jika p-value < alpha
durbinWatsonTest(crimelm)## lag Autocorrelation D-W Statistic p-value
## 1 -0.29815 2.578644 0.104
## Alternative hypothesis: rho != 0Menetapkan matriks pembobot spasial (W)
SPASIAL WEIGHTS
Regresi sebelumnya tidak mempertimbangkan unsur spasial dalam analisisnya. Oleh karena itu, kita akan melakukan uji autokorelasi spasial dengan menggunakan Moran Test. Untuk itu, perlu dicari dan dihitung terlebih dahulu bobot spasialnya.
K-nearest neighboard
W.knn<-knn2nb(knearneigh(coords, k=5 ,longlat=TRUE))
W.knn1 <- nb2listw(W.knn,style='W')
W.knn1## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 190
## Percentage nonzero weights: 13.15789
## Average number of links: 5
## Non-symmetric neighbours list
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 13.44 156.96Inverse distance
distance_matrix <- as.matrix(dist(coords, method="euclidean"))
inv_dist <- 1/distance_matrix
diag(inv_dist) <- 0
W.id<-mat2listw(inv_dist, style="W")
summary(W.id)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 1406
## Percentage nonzero weights: 97.36842
## Average number of links: 37
## Link number distribution:
##
## 37
## 38
## 38 least connected regions:
## 60 71 82 92 93 94 95 96 97 61 62 63 64 65 66 67 68 69 70 72 73 74 75 76 77 78 79 80 81 83 84 85 86 87 88 89 90 91 with 37 links
## 38 most connected regions:
## 60 71 82 92 93 94 95 96 97 61 62 63 64 65 66 67 68 69 70 72 73 74 75 76 77 78 79 80 81 83 84 85 86 87 88 89 90 91 with 37 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 3.384919 153.4568Rook Contiguity
W.rook <- poly2nb(Jatim, row.names=IDs, queen=FALSE)## Warning in poly2nb(Jatim, row.names = IDs, queen = FALSE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.W.rook.s <- nb2listw(W.rook,style='B', zero.policy = TRUE) ; W.rook.s## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 138
## Percentage nonzero weights: 9.556787
## Average number of links: 3.631579
## 2 disjoint connected subgraphs
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 38 1444 138 276 2712Queen Contiguity
W.queen <- poly2nb(Jatim, row.names=IDs, queen=TRUE)## Warning in poly2nb(Jatim, row.names = IDs, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.W.queen.s <- nb2listw(W.queen,style='B', zero.policy = TRUE) ; W.queen.s## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 138
## Percentage nonzero weights: 9.556787
## Average number of links: 3.631579
## 2 disjoint connected subgraphs
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 38 1444 138 276 2712Choosing spatial weights
MI.knn <- moran(y, W.knn1, n=length(W.knn1$neighbours), S0=Szero(W.knn1))
MI.id <- moran(y, W.id, n=length(W.id$neighbours), S0=Szero(W.id))
MI.rook <- moran.test(y, W.rook.s, randomisation = TRUE, zero.policy = TRUE)
MI.queen <- moran.test(y, W.queen.s, randomisation = TRUE, zero.policy = TRUE)
moran_index<-data.frame(
"Weight matrix"=c("KNN (k=5)", "Inverse distance weight", "Rook Contiguity", "Queen Contiguity"),
"Moran's I"=c(MI.knn$I, MI.id$I, 0.22785882,0.22785882 ))
moran_index## Weight.matrix Moran.s.I
## 1 KNN (k=5) 0.084730095
## 2 Inverse distance weight 0.009997162
## 3 Rook Contiguity 0.227858820
## 4 Queen Contiguity 0.227858820Memeriksa apakah ada efek spasial dengan uji Moran’s I
W Optimum
Memilih Weight Optimum didapatkan dari Statistik Moran’s I terbesarnya
W_optimum <- W.rook.smoran.test(y,W_optimum)##
## Moran I test under randomisation
##
## data: y
## weights: W_optimum
##
## Moran I statistic standard deviate = 2.4041, p-value = 0.008107
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.22785882 -0.02702703 0.01124070par(mar = c(5, 5, 2, 2))
moran.plot(y, W_optimum, labels=IDs)
Melakukan pemodelan GWR
Langkah pertama dalam analisis GWR adalah menentukan bandwidth yang akan digunakan.
Adaptive kernel model
library
library(spgwr)## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)Adaptive Gaussian
adapt_gaussian <- gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data,
coords = coords, adapt = T)## Adaptive q: 0.381966 CV score: 3.835687
## Adaptive q: 0.618034 CV score: 3.659894
## Adaptive q: 0.763932 CV score: 3.581449
## Adaptive q: 0.854102 CV score: 3.555415
## Adaptive q: 0.9459191 CV score: 3.529633
## Adaptive q: 0.9108481 CV score: 3.537693
## Adaptive q: 0.9665762 CV score: 3.521865
## Adaptive q: 0.9793429 CV score: 3.51718
## Adaptive q: 0.9872332 CV score: 3.514011
## Adaptive q: 0.9921097 CV score: 3.512258
## Adaptive q: 0.9951235 CV score: 3.511243
## Adaptive q: 0.9969862 CV score: 3.51064
## Adaptive q: 0.9981374 CV score: 3.510276
## Adaptive q: 0.9988488 CV score: 3.510055
## Adaptive q: 0.9992885 CV score: 3.509919
## Adaptive q: 0.9995603 CV score: 3.509836
## Adaptive q: 0.9997282 CV score: 3.509784
## Adaptive q: 0.999832 CV score: 3.509752
## Adaptive q: 0.9998962 CV score: 3.509733
## Adaptive q: 0.9999369 CV score: 3.50972
## Adaptive q: 0.9999369 CV score: 3.50972## Warning in gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords =
## coords, : Bandwidth converged to upper bound:1adapt_gaussian## [1] 0.9999369gwr.model1 = gwr(y ~ x1 + x2 + x3 + x4 + x5 + x6,
data = data,
coords=coords,
adapt=adapt_gaussian,
hatmatrix=TRUE,
se.fit=TRUE) ; gwr.model1 ## Call:
## gwr(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords = coords,
## adapt = adapt_gaussian, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Adaptive quantile: 0.9999369 (about 37 of 38 data points)
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -1.09022891 -1.06550164 -1.02038718 -0.96591665 -0.94162822
## x1 -0.03258835 -0.03204880 -0.03136498 -0.03026810 -0.02901729
## x2 0.00012185 0.00012357 0.00012649 0.00012926 0.00013072
## x3 0.00085772 0.00133087 0.00198858 0.00231946 0.00288038
## x4 0.41784653 0.42281351 0.43156239 0.43610508 0.43853789
## x5 0.02508723 0.02778556 0.03180873 0.03365560 0.03484647
## x6 -0.00444612 -0.00437408 -0.00427879 -0.00412478 -0.00396527
## Global
## X.Intercept. -1.0246
## x1 -0.0298
## x2 0.0001
## x3 0.0013
## x4 0.4262
## x5 0.0287
## x6 -0.0041
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 8.129071
## Effective degrees of freedom (residual: 2traceS - traceS'S): 29.87093
## Sigma (residual: 2traceS - traceS'S): 0.2769322
## Effective number of parameters (model: traceS): 7.595004
## Effective degrees of freedom (model: traceS): 30.405
## Sigma (model: traceS): 0.2744893
## Sigma (ML): 0.2455309
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 24.10671
## AIC (GWR p. 96, eq. 4.22): 8.705049
## Residual sum of squares: 2.290845
## Quasi-global R2: 0.8019674Adaptive bi-square
adapt_bisquare <- gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data,
coords = coords, gweight = gwr.bisquare, adapt = T)## Adaptive q: 0.381966 CV score: 12.25745
## Adaptive q: 0.618034 CV score: 5.799481
## Adaptive q: 0.763932 CV score: 4.674648
## Adaptive q: 0.7659283 CV score: 4.653013
## Adaptive q: 0.8553358 CV score: 4.179043
## Adaptive q: 0.8543961 CV score: 4.184318
## Adaptive q: 0.9105926 CV score: 3.930376
## Adaptive q: 0.9447432 CV score: 3.834364
## Adaptive q: 0.9658494 CV score: 3.812391
## Adaptive q: 0.971543 CV score: 3.809969
## Adaptive q: 0.9779528 CV score: 3.776084
## Adaptive q: 0.9863741 CV score: 3.721964
## Adaptive q: 0.9915787 CV score: 3.694482
## Adaptive q: 0.9947954 CV score: 3.679323
## Adaptive q: 0.9967833 CV score: 3.670562
## Adaptive q: 0.998012 CV score: 3.665362
## Adaptive q: 0.9987713 CV score: 3.662225
## Adaptive q: 0.9992407 CV score: 3.660315
## Adaptive q: 0.9995307 CV score: 3.659146
## Adaptive q: 0.99971 CV score: 3.658427
## Adaptive q: 0.9998207 CV score: 3.657985
## Adaptive q: 0.9998892 CV score: 3.657712
## Adaptive q: 0.9999315 CV score: 3.657543
## Adaptive q: 0.9999315 CV score: 3.657543## Warning in gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords =
## coords, : Bandwidth converged to upper bound:1adapt_bisquare## [1] 0.9999315gwr.model2 = gwr(y ~ x1 + x2 + x3 + x4 + x5 + x6,
data = data,
coords=coords,
adapt=adapt_bisquare,
hatmatrix=TRUE,
se.fit=TRUE) ; gwr.model2## Call:
## gwr(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords = coords,
## adapt = adapt_bisquare, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Adaptive quantile: 0.9999315 (about 37 of 38 data points)
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -1.09023014 -1.06550221 -1.02038714 -0.96591595 -0.94162460
## x1 -0.03258839 -0.03204886 -0.03136509 -0.03026812 -0.02901726
## x2 0.00012185 0.00012357 0.00012649 0.00012926 0.00013072
## x3 0.00085771 0.00133087 0.00198859 0.00231949 0.00288044
## x4 0.41784618 0.42281341 0.43156273 0.43610548 0.43853809
## x5 0.02508707 0.02778554 0.03180895 0.03365576 0.03484652
## x6 -0.00444612 -0.00437409 -0.00427880 -0.00412478 -0.00396526
## Global
## X.Intercept. -1.0246
## x1 -0.0298
## x2 0.0001
## x3 0.0013
## x4 0.4262
## x5 0.0287
## x6 -0.0041
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 8.129105
## Effective degrees of freedom (residual: 2traceS - traceS'S): 29.87089
## Sigma (residual: 2traceS - traceS'S): 0.2769323
## Effective number of parameters (model: traceS): 7.595023
## Effective degrees of freedom (model: traceS): 30.40498
## Sigma (model: traceS): 0.2744892
## Sigma (ML): 0.2455307
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 24.10674
## AIC (GWR p. 96, eq. 4.22): 8.705032
## Residual sum of squares: 2.290843
## Quasi-global R2: 0.8019676Adaptive Tricube
adapt_tricube <- gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 + x6 , data = data,
coords = coords, gweight = gwr.tricube, adapt = T)## Adaptive q: 0.381966 CV score: 14.96297
## Adaptive q: 0.618034 CV score: 5.967062
## Adaptive q: 0.763932 CV score: 4.655802
## Adaptive q: 0.7499279 CV score: 4.705405
## Adaptive q: 0.7998247 CV score: 4.405693
## Adaptive q: 0.8762849 CV score: 4.058567
## Adaptive q: 0.943083 CV score: 3.817575
## Adaptive q: 0.9175684 CV score: 3.891476
## Adaptive q: 0.9648234 CV score: 3.79211
## Adaptive q: 0.9699964 CV score: 3.791426
## Adaptive q: 0.9691233 CV score: 3.791517
## Adaptive q: 0.9814568 CV score: 3.730442
## Adaptive q: 0.9885396 CV score: 3.68241
## Adaptive q: 0.9929171 CV score: 3.656542
## Adaptive q: 0.9956225 CV score: 3.641941
## Adaptive q: 0.9972946 CV score: 3.633417
## Adaptive q: 0.998328 CV score: 3.628331
## Adaptive q: 0.9989666 CV score: 3.625255
## Adaptive q: 0.9993613 CV score: 3.623379
## Adaptive q: 0.9996053 CV score: 3.622229
## Adaptive q: 0.9997561 CV score: 3.621523
## Adaptive q: 0.9998492 CV score: 3.621087
## Adaptive q: 0.9999068 CV score: 3.620818
## Adaptive q: 0.9999475 CV score: 3.620629
## Adaptive q: 0.9999475 CV score: 3.620629## Warning in gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords =
## coords, : Bandwidth converged to upper bound:1adapt_tricube## [1] 0.9999475gwr.model3 = gwr(y ~ x1 + x2 + x3 + x4+ x5 + x6,
data = data,
coords=coords,
adapt=adapt_tricube,
hatmatrix=TRUE,
se.fit=TRUE) ; gwr.model3## Call:
## gwr(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords = coords,
## adapt = adapt_tricube, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Adaptive quantile: 0.9999475 (about 37 of 38 data points)
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -1.09022649 -1.06550052 -1.02038726 -0.96591805 -0.94163537
## x1 -0.03258827 -0.03204867 -0.03136478 -0.03026807 -0.02901735
## x2 0.00012185 0.00012357 0.00012649 0.00012926 0.00013072
## x3 0.00085773 0.00133086 0.00198857 0.00231941 0.00288026
## x4 0.41784722 0.42281371 0.43156172 0.43610429 0.43853751
## x5 0.02508753 0.02778560 0.03180829 0.03365527 0.03484639
## x6 -0.00444612 -0.00437406 -0.00427877 -0.00412478 -0.00396528
## Global
## X.Intercept. -1.0246
## x1 -0.0298
## x2 0.0001
## x3 0.0013
## x4 0.4262
## x5 0.0287
## x6 -0.0041
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 8.129003
## Effective degrees of freedom (residual: 2traceS - traceS'S): 29.871
## Sigma (residual: 2traceS - traceS'S): 0.2769322
## Effective number of parameters (model: traceS): 7.594967
## Effective degrees of freedom (model: traceS): 30.40503
## Sigma (model: traceS): 0.2744894
## Sigma (ML): 0.2455311
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 24.10666
## AIC (GWR p. 96, eq. 4.22): 8.705083
## Residual sum of squares: 2.290849
## Quasi-global R2: 0.8019671Fixed kernel model
Gaussian
bw4 <- gwr.sel(y ~ x1 + x2 + x3 + x4 +x5 + x6, data = data, coords = coords)## Bandwidth: 1.277611 CV score: 3.544682
## Bandwidth: 2.065154 CV score: 3.485895
## Bandwidth: 2.551882 CV score: 3.481357
## Bandwidth: 2.399479 CV score: 3.482166
## Bandwidth: 2.696807 CV score: 3.48087
## Bandwidth: 2.942265 CV score: 3.480438
## Bandwidth: 3.033769 CV score: 3.48036
## Bandwidth: 3.146203 CV score: 3.480306
## Bandwidth: 3.220007 CV score: 3.48029
## Bandwidth: 3.260134 CV score: 3.480286
## Bandwidth: 3.277209 CV score: 3.480286
## Bandwidth: 3.281796 CV score: 3.480286
## Bandwidth: 3.282538 CV score: 3.480286
## Bandwidth: 3.282599 CV score: 3.480286
## Bandwidth: 3.28264 CV score: 3.480286
## Bandwidth: 3.282599 CV score: 3.480286bw4## [1] 3.282599gwr.model4 = gwr(y ~ x1 + x2 + x3 + x4 +x5 + x6,
data = data,
coords=coords,
bandwidth = bw4,
hatmatrix=TRUE,
se.fit=TRUE) ; gwr.model4## Call:
## gwr(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords = coords,
## bandwidth = bw4, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Fixed bandwidth: 3.282599
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -1.08101775 -1.04305187 -1.02221780 -0.99961396 -0.96026927
## x1 -0.03168505 -0.03079823 -0.03035027 -0.02992680 -0.02918283
## x2 0.00012299 0.00012543 0.00012671 0.00012794 0.00013019
## x3 0.00091380 0.00129846 0.00159560 0.00175242 0.00226992
## x4 0.41966802 0.42469847 0.42797756 0.43099554 0.43636285
## x5 0.02588809 0.02825898 0.02976174 0.03112701 0.03359106
## x6 -0.00434042 -0.00423168 -0.00417150 -0.00410903 -0.00399582
## Global
## X.Intercept. -1.0246
## x1 -0.0298
## x2 0.0001
## x3 0.0013
## x4 0.4262
## x5 0.0287
## x6 -0.0041
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 7.577627
## Effective degrees of freedom (residual: 2traceS - traceS'S): 30.42237
## Sigma (residual: 2traceS - traceS'S): 0.2763181
## Effective number of parameters (model: traceS): 7.297397
## Effective degrees of freedom (model: traceS): 30.7026
## Sigma (model: traceS): 0.2750542
## Sigma (ML): 0.2472373
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 23.60664
## AIC (GWR p. 96, eq. 4.22): 8.933833
## Residual sum of squares: 2.3228
## Quasi-global R2: 0.7992051Bi-square
bw5 <- gwr.sel(y ~ x1 + x2 + x3 + x4+x5 + x6, data = data,
coords = coords, gweight = gwr.bisquare)## Bandwidth: 1.277611 CV score: NA## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value## Bandwidth: 2.065154 CV score: 3.876342
## Bandwidth: 2.551882 CV score: 3.606655
## Bandwidth: 2.308518 CV score: 3.722263
## Bandwidth: 2.852697 CV score: 3.510528
## Bandwidth: 3.03861 CV score: 3.481488
## Bandwidth: 3.178348 CV score: 3.472824
## Bandwidth: 3.215659 CV score: 3.471501
## Bandwidth: 3.262933 CV score: 3.470215
## Bandwidth: 3.292151 CV score: 3.469598
## Bandwidth: 3.310208 CV score: 3.469273
## Bandwidth: 3.321368 CV score: 3.469092
## Bandwidth: 3.328265 CV score: 3.468987
## Bandwidth: 3.332528 CV score: 3.468925
## Bandwidth: 3.335162 CV score: 3.468887
## Bandwidth: 3.336791 CV score: 3.468864
## Bandwidth: 3.337797 CV score: 3.46885
## Bandwidth: 3.338419 CV score: 3.468842
## Bandwidth: 3.338803 CV score: 3.468836
## Bandwidth: 3.339041 CV score: 3.468833
## Bandwidth: 3.339187 CV score: 3.468831
## Bandwidth: 3.339278 CV score: 3.46883
## Bandwidth: 3.339334 CV score: 3.468829
## Bandwidth: 3.339375 CV score: 3.468829
## Bandwidth: 3.339375 CV score: 3.468829## Warning in gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords =
## coords, : Bandwidth converged to upper bound:3.33942504384475bw5## [1] 3.339375gwr.model5<-gwr(y ~ x1 + x2 + x3 + x4+x5 + x6,
data = data,
coords=coords,
bandwidth = bw5,
gweight = gwr.bisquare,
se.fit=T, hatmatrix=T) ; gwr.model5## Call:
## gwr(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords = coords,
## bandwidth = bw5, gweight = gwr.bisquare, hatmatrix = T, se.fit = T)
## Kernel function: gwr.bisquare
## Fixed bandwidth: 3.339375
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -1.34691921 -1.11650541 -1.01883611 -0.91938050 -0.69501667
## x1 -0.04826314 -0.03612184 -0.03279481 -0.03095697 -0.02707438
## x2 0.00010654 0.00012058 0.00012632 0.00013192 0.00014571
## x3 0.00106453 0.00185140 0.00257714 0.00322717 0.00575466
## x4 0.39688474 0.42030610 0.43442469 0.44949447 0.49065612
## x5 0.01546553 0.02690887 0.03378031 0.04111422 0.06411439
## x6 -0.00560927 -0.00469972 -0.00440058 -0.00412965 -0.00364457
## Global
## X.Intercept. -1.0246
## x1 -0.0298
## x2 0.0001
## x3 0.0013
## x4 0.4262
## x5 0.0287
## x6 -0.0041
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 9.445521
## Effective degrees of freedom (residual: 2traceS - traceS'S): 28.55448
## Sigma (residual: 2traceS - traceS'S): 0.2765467
## Effective number of parameters (model: traceS): 8.386073
## Effective degrees of freedom (model: traceS): 29.61393
## Sigma (model: traceS): 0.2715549
## Sigma (ML): 0.2397253
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 25.12412
## AIC (GWR p. 96, eq. 4.22): 7.677519
## Residual sum of squares: 2.183792
## Quasi-global R2: 0.8112216Tri-cube
bw6 <- gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 +x6, data = data,
coords = coords, gweight = gwr.tricube)## Bandwidth: 1.277611 CV score: NA## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value## Bandwidth: 2.065154 CV score: 3.841964
## Bandwidth: 2.551882 CV score: 3.591722
## Bandwidth: 2.308518 CV score: 3.718769
## Bandwidth: 2.852697 CV score: 3.484216
## Bandwidth: 3.03861 CV score: 3.44378
## Bandwidth: 3.153511 CV score: 3.430226
## Bandwidth: 3.274335 CV score: 3.422965
## Bandwidth: 3.228184 CV score: 3.425061
## Bandwidth: 3.299197 CV score: 3.422124
## Bandwidth: 3.314563 CV score: 3.421693
## Bandwidth: 3.324059 CV score: 3.421458
## Bandwidth: 3.329929 CV score: 3.421325
## Bandwidth: 3.333556 CV score: 3.421247
## Bandwidth: 3.335798 CV score: 3.4212
## Bandwidth: 3.337183 CV score: 3.421172
## Bandwidth: 3.33804 CV score: 3.421155
## Bandwidth: 3.338569 CV score: 3.421144
## Bandwidth: 3.338896 CV score: 3.421138
## Bandwidth: 3.339098 CV score: 3.421134
## Bandwidth: 3.339223 CV score: 3.421131
## Bandwidth: 3.3393 CV score: 3.42113
## Bandwidth: 3.339348 CV score: 3.421129
## Bandwidth: 3.339348 CV score: 3.421129## Warning in gwr.sel(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords =
## coords, : Bandwidth converged to upper bound:3.33942504384475bw6## [1] 3.339348gwr.model6<-gwr(y ~ x1 + x2 + x3 + x4 + x5 + x6,
data = data,
coords=coords,
bandwidth = bw6,
gweight = gwr.tricube,
se.fit=T, hatmatrix=T) ; gwr.model6## Call:
## gwr(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = data, coords = coords,
## bandwidth = bw6, gweight = gwr.tricube, hatmatrix = T, se.fit = T)
## Kernel function: gwr.tricube
## Fixed bandwidth: 3.339348
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -1.35921452 -1.11541025 -1.01395150 -0.91003863 -0.71660631
## x1 -0.05099429 -0.03641636 -0.03213486 -0.02987671 -0.02642712
## x2 0.00010730 0.00011954 0.00012613 0.00013236 0.00014641
## x3 0.00132534 0.00193495 0.00243316 0.00297555 0.00612822
## x4 0.39894869 0.41980377 0.43141486 0.44824670 0.49290275
## x5 0.01667125 0.02625218 0.03212489 0.04008481 0.06820361
## x6 -0.00571571 -0.00463953 -0.00430993 -0.00409892 -0.00370199
## Global
## X.Intercept. -1.0246
## x1 -0.0298
## x2 0.0001
## x3 0.0013
## x4 0.4262
## x5 0.0287
## x6 -0.0041
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 8.589249
## Effective degrees of freedom (residual: 2traceS - traceS'S): 29.41075
## Sigma (residual: 2traceS - traceS'S): 0.2748753
## Effective number of parameters (model: traceS): 7.920993
## Effective degrees of freedom (model: traceS): 30.07901
## Sigma (model: traceS): 0.2718047
## Sigma (ML): 0.2418226
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 24.09946
## AIC (GWR p. 96, eq. 4.22): 7.874457
## Residual sum of squares: 2.222171
## Quasi-global R2: 0.807904BFC02.gwr.test(gwr.model1)##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA
##
## data: gwr.model1
## F = 1.0331, df1 = 31.000, df2 = 29.871, p-value = 0.4654
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 2.366717 2.290845BFC02.gwr.test(gwr.model2)##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA
##
## data: gwr.model2
## F = 1.0331, df1 = 31.000, df2 = 29.871, p-value = 0.4654
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 2.366717 2.290843BFC02.gwr.test(gwr.model3)##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA
##
## data: gwr.model3
## F = 1.0331, df1 = 31.000, df2 = 29.871, p-value = 0.4654
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 2.366717 2.290849BFC02.gwr.test(gwr.model4)##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA
##
## data: gwr.model4
## F = 1.0189, df1 = 31.000, df2 = 30.422, p-value = 0.4799
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 2.366717 2.322800BFC02.gwr.test(gwr.model5)##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA
##
## data: gwr.model5
## F = 1.0838, df1 = 31.000, df2 = 28.554, p-value = 0.416
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 2.366717 2.183792BFC02.gwr.test(gwr.model6)##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA
##
## data: gwr.model6
## F = 1.065, df1 = 31.000, df2 = 29.411, p-value = 0.4331
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 2.366717 2.222171Plotting residuals model GWR
results2 <-as.data.frame(gwr.model2$SDF)
head(results2)## sum.w X.Intercept. x1 x2 x3 x4
## 1 32.65265 -0.9548101 -0.02971345 0.0001222592 0.002880444 0.4224349
## 2 30.76056 -0.9611221 -0.02958802 0.0001235342 0.001917716 0.4187790
## 3 33.51984 -1.0264778 -0.03205498 0.0001268705 0.002002818 0.4330901
## 4 33.35846 -1.0657559 -0.03252041 0.0001293131 0.001306284 0.4374622
## 5 34.37289 -1.0521564 -0.03139705 0.0001282209 0.001626588 0.4337353
## 6 31.83697 -0.9530801 -0.02950313 0.0001229271 0.002140946 0.4184205
## x5 x6 X.Intercept._se x1_se x2_se x3_se
## 1 0.02744854 -0.004080680 0.4335868 0.03670999 3.296123e-05 0.02160274
## 2 0.02576067 -0.004012704 0.4334143 0.03672743 3.300801e-05 0.02163867
## 3 0.03303662 -0.004367822 0.4317156 0.03661257 3.284925e-05 0.02158917
## 4 0.03484652 -0.004446124 0.4310174 0.03665412 3.281577e-05 0.02158555
## 5 0.03243686 -0.004279567 0.4299761 0.03659679 3.277167e-05 0.02158775
## 6 0.02559239 -0.004004059 0.4340663 0.03673469 3.303077e-05 0.02163345
## x4_se x5_se x6_se gwr.e pred pred.se localR2
## 1 0.05928811 0.02134347 0.001698083 0.05920939 0.6339378 0.21797881 0.8052679
## 2 0.05916470 0.02132108 0.001700964 -0.05485056 1.8466100 0.09293542 0.8045983
## 3 0.05927952 0.02135645 0.001700337 -0.60894746 0.6089475 0.10497478 0.7964849
## 4 0.05925091 0.02137251 0.001701485 0.57599331 1.2157662 0.06452976 0.7935800
## 5 0.05926277 0.02129173 0.001699945 0.03691897 1.0616933 0.23403510 0.7985192
## 6 0.05911050 0.02132964 0.001699730 -0.20235941 0.8955066 0.07285156 0.8051239
## X.Intercept._se_EDF x1_se_EDF x2_se_EDF x3_se_EDF x4_se_EDF x5_se_EDF
## 1 0.4374459 0.03703671 3.325459e-05 0.02179501 0.05981578 0.02153343
## 2 0.4372718 0.03705431 3.330179e-05 0.02183126 0.05969128 0.02151084
## 3 0.4355580 0.03693843 3.314161e-05 0.02178132 0.05980712 0.02154652
## 4 0.4348535 0.03698035 3.310784e-05 0.02177766 0.05977826 0.02156273
## 5 0.4338030 0.03692251 3.306334e-05 0.02177988 0.05979022 0.02148123
## 6 0.4379296 0.03706164 3.332475e-05 0.02182599 0.05963660 0.02151947
## x6_se_EDF pred.se.1 coord.x coord.y
## 1 0.001713196 0.21991888 112.9302 -7.043867
## 2 0.001716103 0.09376257 114.2054 -8.364818
## 3 0.001715470 0.10590909 112.5311 -7.833109
## 4 0.001716629 0.06510409 112.2368 -8.129433
## 5 0.001715075 0.23611807 111.8098 -7.255433
## 6 0.001714858 0.07349995 113.9476 -7.944023resid_gwr2 <- results2$gwr.e
Jatim$resid_gwr2 <- resid_gwr2
spplot(Jatim,zcol="resid_gwr2",main="Peta Sebaran Residual Model Adaptive bi-square")
Residuals diagnostics
Sama seperti OLS, pengujian menggunakan shapiro-wilk, Breusch-pagan dan autokorelasi residual dengan three moment approximation untuk Moran’s I, hasil ketiganya harus tidak signifikan untuk asumsi residual terpenuhi
NORMALITY
shapiro.test(resid_gwr2) #Normality##
## Shapiro-Wilk normality test
##
## data: resid_gwr2
## W = 0.9844, p-value = 0.8636HETEROKEDASTICIY
bptest(gwr.model2$lm, weights = gwr.model2$gweight) #Heteroskedasticity##
## studentized Breusch-Pagan test
##
## data: gwr.model2$lm
## BP = 6.0998, df = 6, p-value = 0.4121NON AUTOKORELASI
gwr.morantest(gwr.model2, W_optimum, zero.policy = TRUE) #Non-autocorrelation## Warning in gwr.morantest(gwr.model2, W_optimum, zero.policy = TRUE):
## W_optimumnot row standardised##
## Leung et al. 2000 three moment approximation for Moran's I
##
## data: GWR residuals
## statistic = 426.15, df = 419.43, p-value = 0.3997
## sample estimates:
## I
## -0.1277619Parameter Significance Test
pengujian terhadap masing-masing parameter juga dilakukan untuk melihat apakah setiap parameter yang diestimasi signifikan terhadap model 𝑦𝑖. Pengujian tersebut dilakukan dengan menggunakan uji F(3) yang juga dikemukakan oleh Leung dkk.
LMZ.F3GWR.test(gwr.model2)##
## Leung et al. (2000) F(3) test
##
## F statistic Numerator d.f. Denominator d.f. Pr(>)
## (Intercept) 1.39606 8.30521 30.899 0.235746
## x1 0.12071 7.67975 30.899 0.997549
## x2 0.87378 13.12551 30.899 0.587163
## x3 0.27469 9.78830 30.899 0.981332
## x4 2.28621 10.53953 30.899 0.036317 *
## x5 3.39439 11.78552 30.899 0.003147 **
## x6 1.23411 11.60969 30.899 0.306004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparing AIC
data.frame("MODEL" = c("GWR","OLS"),
"AIC" = c(gwr.model2$results$AICh,AIC(crimelm)),
"R2" = c(0.8019676 ,0.7954))%>% arrange(AIC)## MODEL AIC R2
## 1 GWR 8.705032 0.8019676
## 2 OLS 18.348201 0.7954000library(tmap)## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')results2 <-as.data.frame(gwr.model2$SDF)
gwr.map <- cbind(Jatim, as.matrix(results2))
qtm(gwr.map, fill = "localR2")
gwr.map2 <- st_as_sf(gwr.map)#Attach coefficients to original dataframe
Jatim$coefx1<-results2$x1
Jatim$coefx2<-results2$x2
Jatim$coefx3<-results2$x3
Jatim$coefx4<-results2$x4
Jatim$coefx5<-results2$x5
Jatim$coefx6<-results2$x6#Spatial distribution of x1
map1 <- tm_shape(gwr.map2) +
tm_fill("x1",
n = 5,
style = "quantile",
title = "x1 Coefficient (tingkat pengangguran terbuka)") +
tm_borders(col = "black", lwd = 0.5)+
tm_layout(frame = FALSE,
legend.text.size = 0.5,
legend.title.size = 0.6)
map1
#Spatial distribution of x2
map2 <- tm_shape(gwr.map2) +
tm_fill("x2",
n = 5,
style = "quantile",
title = "x2 Coefficient (pengeluaran per kapita)") +
tm_borders(col = "black", lwd = 0.5)+
tm_layout(frame = FALSE,
legend.text.size = 0.5,
legend.title.size = 0.6)
map2
#Spatial distribution of x3
map3 <- tm_shape(gwr.map2) +
tm_fill("x3",
n = 5,
style = "quantile",
title = "x3 Coefficient (pertumbuhan ekonomi)") +
tm_borders(col = "black", lwd = 0.5)+
tm_layout(frame = FALSE,
legend.text.size = 0.5,
legend.title.size = 0.6)
map3
#Spatial distribution of x4
map4 <- tm_shape(gwr.map2) +
tm_fill("x4",
n = 5,
style = "quantile",
title = "x4 Coefficient (kepadatan penduduk)",
group = IDs) +
tm_borders(col = "black", lwd = 0.5)+
tm_layout(frame = FALSE,
legend.text.size = 0.5,
legend.title.size = 0.6)
map4
#Spatial distribution of x5
map5 <- tm_shape(gwr.map2) +
tm_fill("x5",
n = 5,
style = "quantile",
title = "x5 Coefficient (% belum / tidak sekolah",
group = IDs) +
tm_borders(col = "black", lwd = 0.5)+
tm_layout(frame = FALSE,
legend.text.size = 0.5,
legend.title.size = 0.6)
map5
#Spatial distribution of x4
map6 <- tm_shape(gwr.map2) +
tm_fill("x4",
n = 5,
style = "quantile",
title = "x6 koefisien (jumlah penduduk miskin)",
group = IDs) +
tm_borders(col = "black", lwd = 0.5)+
tm_layout(frame = FALSE,
legend.text.size = 0.5,
legend.title.size = 0.6)
map6
#Local significance test
t1 = gwr.model2$SDF$x1 / gwr.model2$SDF$x1_se
t2 = gwr.model2$SDF$x2 / gwr.model2$SDF$x2_se
t3 = gwr.model2$SDF$x3 / gwr.model2$SDF$x3_se
t4 = gwr.model2$SDF$x4 / gwr.model2$SDF$x4_se
t5 = gwr.model2$SDF$x5 / gwr.model2$SDF$x5_se
t6 = gwr.model2$SDF$x6 / gwr.model2$SDF$x6_se
pvalue1 = 2*pt(q=t1, df=37, lower.tail=TRUE)
pvalue2 = 2*pt(q=t2, df=37, lower.tail=TRUE)
pvalue3 = 2*pt(q=t3, df=37, lower.tail=TRUE)
pvalue4 = 2*pt(q=t4, df=37, lower.tail=TRUE)
pvalue5 = 2*pt(q=t5, df=37, lower.tail=TRUE)
pvalue6 = 2*pt(q=t6, df=37, lower.tail=TRUE)
data.frame(IDs, t1, pvalue1) %>% mutate(ket=ifelse(pvalue1<0.05,"Reject", "Not Reject"))## IDs t1 pvalue1 ket
## 1 Bangkalan -0.8094105 0.4234539 Not Reject
## 2 Banyuwangi -0.8056108 0.4256146 Not Reject
## 3 Batu -0.8755183 0.3869404 Not Reject
## 4 Blitar -0.8872239 0.3806894 Not Reject
## 5 Bojonegoro -0.8579181 0.3964612 Not Reject
## 6 Bondowoso -0.8031409 0.4270227 Not Reject
## 7 Gresik -0.8364115 0.4082930 Not Reject
## 8 Jember -0.8150057 0.4202844 Not Reject
## 9 Jombang -0.8638412 0.3932408 Not Reject
## 10 Kediri -0.8760964 0.3866302 Not Reject
## 11 Kota Blitar -0.8852682 0.3817292 Not Reject
## 12 Kota Kediri -0.8756923 0.3868470 Not Reject
## 13 Kota Madiun -0.8679670 0.3910073 Not Reject
## 14 Kota Malang -0.8841617 0.3823184 Not Reject
## 15 Kota Mojokerto -0.8566187 0.3971699 Not Reject
## 16 Kota Pasuruan -0.8378570 0.4074910 Not Reject
## 17 Kota Probolinggo -0.8219312 0.4163814 Not Reject
## 18 Lamongan -0.8449525 0.4035682 Not Reject
## 19 Lumajang -0.8384256 0.4071757 Not Reject
## 20 Madiun -0.8684167 0.3907644 Not Reject
## 21 Magetan -0.8677061 0.3911483 Not Reject
## 22 Malang -0.8878178 0.3803740 Not Reject
## 23 Mojokerto -0.8595331 0.3955814 Not Reject
## 24 Nganjuk -0.8681358 0.3909161 Not Reject
## 25 Ngawi -0.8631755 0.3936019 Not Reject
## 26 Pacitan -0.8729328 0.3883299 Not Reject
## 27 Pamekasan -0.7943914 0.4320336 Not Reject
## 28 Pasuruan -0.8514931 0.3999731 Not Reject
## 29 Ponorogo -0.8738542 0.3878343 Not Reject
## 30 Probolinggo -0.8190816 0.4179846 Not Reject
## 31 Sampang -0.7987359 0.4295410 Not Reject
## 32 Sidoarjo -0.8449275 0.4035819 Not Reject
## 33 Situbondo -0.7993725 0.4291766 Not Reject
## 34 Sumenep -0.7902367 0.4344253 Not Reject
## 35 Surabaya -0.8327815 0.4103114 Not Reject
## 36 Trenggalek -0.8791269 0.3850064 Not Reject
## 37 Tuban -0.8479167 0.4019364 Not Reject
## 38 Tulungagung -0.8819151 0.3835164 Not Rejectdata.frame(IDs, t2, pvalue2) %>% mutate(ket=ifelse(pvalue2<0.05,"Reject", "Not Reject"))## IDs t2 pvalue2 ket
## 1 Bangkalan 3.709183 1.999321 Not Reject
## 2 Banyuwangi 3.742553 1.999383 Not Reject
## 3 Batu 3.862205 1.999564 Not Reject
## 4 Blitar 3.940577 1.999653 Not Reject
## 5 Bojonegoro 3.912555 1.999623 Not Reject
## 6 Bondowoso 3.721593 1.999344 Not Reject
## 7 Gresik 3.800858 1.999478 Not Reject
## 8 Jember 3.742801 1.999383 Not Reject
## 9 Jombang 3.886925 1.999594 Not Reject
## 10 Kediri 3.931499 1.999643 Not Reject
## 11 Kota Blitar 3.944833 1.999657 Not Reject
## 12 Kota Kediri 3.937248 1.999649 Not Reject
## 13 Kota Madiun 3.949239 1.999662 Not Reject
## 14 Kota Malang 3.846132 1.999543 Not Reject
## 15 Kota Mojokerto 3.849444 1.999547 Not Reject
## 16 Kota Pasuruan 3.748423 1.999393 Not Reject
## 17 Kota Probolinggo 3.731739 1.999363 Not Reject
## 18 Lamongan 3.848928 1.999546 Not Reject
## 19 Lumajang 3.767192 1.999425 Not Reject
## 20 Madiun 3.944629 1.999657 Not Reject
## 21 Magetan 3.955881 1.999668 Not Reject
## 22 Malang 3.862023 1.999563 Not Reject
## 23 Mojokerto 3.845788 1.999542 Not Reject
## 24 Nganjuk 3.926181 1.999638 Not Reject
## 25 Ngawi 3.945215 1.999658 Not Reject
## 26 Pacitan 3.978415 1.999689 Not Reject
## 27 Pamekasan 3.693090 1.999288 Not Reject
## 28 Pasuruan 3.768116 1.999427 Not Reject
## 29 Ponorogo 3.966056 1.999678 Not Reject
## 30 Probolinggo 3.731464 1.999363 Not Reject
## 31 Sampang 3.696079 1.999295 Not Reject
## 32 Sidoarjo 3.771845 1.999433 Not Reject
## 33 Situbondo 3.715555 1.999333 Not Reject
## 34 Sumenep 3.693913 1.999290 Not Reject
## 35 Surabaya 3.751886 1.999399 Not Reject
## 36 Trenggalek 3.974795 1.999686 Not Reject
## 37 Tuban 3.887753 1.999595 Not Reject
## 38 Tulungagung 3.964077 1.999676 Not Rejectdata.frame(IDs, t3, pvalue3) %>% mutate(ket=ifelse(pvalue3<0.05,"Reject", "Not Reject"))## IDs t3 pvalue3 ket
## 1 Bangkalan 0.13333698 1.105351 Not Reject
## 2 Banyuwangi 0.08862448 1.070142 Not Reject
## 3 Batu 0.09276955 1.073413 Not Reject
## 4 Blitar 0.06051663 1.047930 Not Reject
## 5 Bojonegoro 0.07534773 1.059656 Not Reject
## 6 Bondowoso 0.09896461 1.078299 Not Reject
## 7 Gresik 0.11606311 1.091770 Not Reject
## 8 Jember 0.09136704 1.072306 Not Reject
## 9 Jombang 0.08739476 1.069171 Not Reject
## 10 Kediri 0.06781150 1.053699 Not Reject
## 11 Kota Blitar 0.05941412 1.047058 Not Reject
## 12 Kota Kediri 0.06507208 1.051533 Not Reject
## 13 Kota Madiun 0.05762920 1.045646 Not Reject
## 14 Kota Malang 0.09256935 1.073255 Not Reject
## 15 Kota Mojokerto 0.10085235 1.079788 Not Reject
## 16 Kota Pasuruan 0.11652233 1.092131 Not Reject
## 17 Kota Probolinggo 0.10830976 1.085665 Not Reject
## 18 Lamongan 0.10074794 1.079705 Not Reject
## 19 Lumajang 0.09286442 1.073488 Not Reject
## 20 Madiun 0.06046635 1.047890 Not Reject
## 21 Magetan 0.05342938 1.042323 Not Reject
## 22 Malang 0.08228573 1.065137 Not Reject
## 23 Mojokerto 0.10161905 1.080392 Not Reject
## 24 Nganjuk 0.07035128 1.055707 Not Reject
## 25 Ngawi 0.05858487 1.046402 Not Reject
## 26 Pacitan 0.03971868 1.031469 Not Reject
## 27 Pamekasan 0.12163432 1.096153 Not Reject
## 28 Pasuruan 0.11292444 1.089299 Not Reject
## 29 Ponorogo 0.04849380 1.038416 Not Reject
## 30 Probolinggo 0.10415343 1.082390 Not Reject
## 31 Sampang 0.12633621 1.099850 Not Reject
## 32 Sidoarjo 0.12262348 1.096931 Not Reject
## 33 Situbondo 0.10154824 1.080336 Not Reject
## 34 Sumenep 0.11739507 1.092818 Not Reject
## 35 Surabaya 0.12873145 1.101733 Not Reject
## 36 Trenggalek 0.04340844 1.034391 Not Reject
## 37 Tuban 0.08479485 1.067118 Not Reject
## 38 Tulungagung 0.05007118 1.039665 Not Rejectdata.frame(IDs, t4, pvalue4) %>% mutate(ket=ifelse(pvalue4<0.05,"Reject", "Not Reject"))## IDs t4 pvalue4 ket
## 1 Bangkalan 7.125120 2 Not Reject
## 2 Banyuwangi 7.078190 2 Not Reject
## 3 Batu 7.305897 2 Not Reject
## 4 Blitar 7.383214 2 Not Reject
## 5 Bojonegoro 7.318849 2 Not Reject
## 6 Bondowoso 7.078615 2 Not Reject
## 7 Gresik 7.211736 2 Not Reject
## 8 Jember 7.102292 2 Not Reject
## 9 Jombang 7.307840 2 Not Reject
## 10 Kediri 7.359669 2 Not Reject
## 11 Kota Blitar 7.384096 2 Not Reject
## 12 Kota Kediri 7.363762 2 Not Reject
## 13 Kota Madiun 7.363121 2 Not Reject
## 14 Kota Malang 7.303419 2 Not Reject
## 15 Kota Mojokerto 7.271336 2 Not Reject
## 16 Kota Pasuruan 7.183578 2 Not Reject
## 17 Kota Probolinggo 7.140544 2 Not Reject
## 18 Lamongan 7.254802 2 Not Reject
## 19 Lumajang 7.173035 2 Not Reject
## 20 Madiun 7.359630 2 Not Reject
## 21 Magetan 7.368876 2 Not Reject
## 22 Malang 7.319347 2 Not Reject
## 23 Mojokerto 7.272883 2 Not Reject
## 24 Nganjuk 7.343795 2 Not Reject
## 25 Ngawi 7.353417 2 Not Reject
## 26 Pacitan 7.397361 2 Not Reject
## 27 Pamekasan 7.087630 2 Not Reject
## 28 Pasuruan 7.212272 2 Not Reject
## 29 Ponorogo 7.386447 2 Not Reject
## 30 Probolinggo 7.128185 2 Not Reject
## 31 Sampang 7.101491 2 Not Reject
## 32 Sidoarjo 7.204422 2 Not Reject
## 33 Situbondo 7.072961 2 Not Reject
## 34 Sumenep 7.075479 2 Not Reject
## 35 Surabaya 7.176266 2 Not Reject
## 36 Trenggalek 7.402022 2 Not Reject
## 37 Tuban 7.286639 2 Not Reject
## 38 Tulungagung 7.396039 2 Not Rejectdata.frame(IDs, t5, pvalue5) %>% mutate(ket=ifelse(pvalue5<0.05,"Reject", "Not Reject"))## IDs t5 pvalue5 ket
## 1 Bangkalan 1.286040 1.793581 Not Reject
## 2 Banyuwangi 1.208226 1.765374 Not Reject
## 3 Batu 1.546916 1.869605 Not Reject
## 4 Blitar 1.630437 1.888506 Not Reject
## 5 Bojonegoro 1.523449 1.863851 Not Reject
## 6 Bondowoso 1.199851 1.762176 Not Reject
## 7 Gresik 1.410519 1.833261 Not Reject
## 8 Jember 1.243422 1.778466 Not Reject
## 9 Jombang 1.528672 1.865149 Not Reject
## 10 Kediri 1.590544 1.879779 Not Reject
## 11 Kota Blitar 1.627533 1.887889 Not Reject
## 12 Kota Kediri 1.592968 1.880324 Not Reject
## 13 Kota Madiun 1.576067 1.876477 Not Reject
## 14 Kota Malang 1.558803 1.872443 Not Reject
## 15 Kota Mojokerto 1.489235 1.855099 Not Reject
## 16 Kota Pasuruan 1.368705 1.820655 Not Reject
## 17 Kota Probolinggo 1.288238 1.794339 Not Reject
## 18 Lamongan 1.454954 1.845886 Not Reject
## 19 Lumajang 1.344753 1.813109 Not Reject
## 20 Madiun 1.574679 1.876156 Not Reject
## 21 Magetan 1.579422 1.877248 Not Reject
## 22 Malang 1.574920 1.876212 Not Reject
## 23 Mojokerto 1.495461 1.856724 Not Reject
## 24 Nganjuk 1.562915 1.873414 Not Reject
## 25 Ngawi 1.557947 1.872241 Not Reject
## 26 Pacitan 1.612830 1.884720 Not Reject
## 27 Pamekasan 1.198346 1.761598 Not Reject
## 28 Pasuruan 1.420233 1.836088 Not Reject
## 29 Ponorogo 1.606517 1.883338 Not Reject
## 30 Probolinggo 1.270436 1.788140 Not Reject
## 31 Sampang 1.225954 1.772038 Not Reject
## 32 Sidoarjo 1.421046 1.836323 Not Reject
## 33 Situbondo 1.187663 1.757465 Not Reject
## 34 Sumenep 1.176078 1.752925 Not Reject
## 35 Surabaya 1.380496 1.824282 Not Reject
## 36 Trenggalek 1.630486 1.888516 Not Reject
## 37 Tuban 1.479968 1.852653 Not Reject
## 38 Tulungagung 1.631216 1.888671 Not Rejectdata.frame(IDs, t6, pvalue6) %>% mutate(ket=ifelse(pvalue6<0.05,"Reject", "Not Reject"))## IDs t6 pvalue6 ket
## 1 Bangkalan -2.403110 0.02138512 Reject
## 2 Banyuwangi -2.359076 0.02371057 Reject
## 3 Batu -2.568798 0.01437336 Reject
## 4 Blitar -2.613084 0.01289540 Reject
## 5 Bojonegoro -2.517474 0.01627983 Reject
## 6 Bondowoso -2.355703 0.02389781 Reject
## 7 Gresik -2.466677 0.01839182 Reject
## 8 Jember -2.386655 0.02222884 Reject
## 9 Jombang -2.536655 0.01554181 Reject
## 10 Kediri -2.574557 0.01417277 Reject
## 11 Kota Blitar -2.607105 0.01308642 Reject
## 12 Kota Kediri -2.573623 0.01420512 Reject
## 13 Kota Madiun -2.550052 0.01504460 Reject
## 14 Kota Malang -2.587750 0.01372279 Reject
## 15 Kota Mojokerto -2.517874 0.01626411 Reject
## 16 Kota Pasuruan -2.468048 0.01833166 Reject
## 17 Kota Probolinggo -2.418643 0.02061546 Reject
## 18 Lamongan -2.485103 0.01759841 Reject
## 19 Lumajang -2.458171 0.01876895 Reject
## 20 Madiun -2.550755 0.01501894 Reject
## 21 Magetan -2.550698 0.01502102 Reject
## 22 Malang -2.601314 0.01327391 Reject
## 23 Mojokerto -2.525922 0.01595091 Reject
## 24 Nganjuk -2.548735 0.01509282 Reject
## 25 Ngawi -2.535150 0.01559856 Reject
## 26 Pacitan -2.575270 0.01414808 Reject
## 27 Pamekasan -2.352512 0.02407617 Reject
## 28 Pasuruan -2.501687 0.01691121 Reject
## 29 Ponorogo -2.572562 0.01424195 Reject
## 30 Probolinggo -2.406949 0.02119253 Reject
## 31 Sampang -2.369510 0.02313978 Reject
## 32 Sidoarjo -2.488973 0.01743582 Reject
## 33 Situbondo -2.346818 0.02439739 Reject
## 34 Sumenep -2.336891 0.02496668 Reject
## 35 Surabaya -2.460045 0.01868525 Reject
## 36 Trenggalek -2.593421 0.01353345 Reject
## 37 Tuban -2.489559 0.01741132 Reject
## 38 Tulungagung -2.599505 0.01333298 Reject