library(spdep)
## Warning: package 'spdep' was built under R version 4.3.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.3.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.3.3
library(RColorBrewer)
library(lattice)
library(gstat)
## Warning: package 'gstat' was built under R version 4.3.3
library(raster)
## Warning: package 'raster' was built under R version 4.3.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.2
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
## 
##     extract
library(tmap)
## Warning: package 'tmap' was built under R version 4.3.3
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(sarima)
## Warning: package 'sarima' was built under R version 4.3.3
## Loading required package: stats4
## 
## Attaching package: 'sarima'
## The following object is masked from 'package:stats':
## 
##     spectrum
library(extrafont)
## Registering fonts with R
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.2
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(spgwr)
## Warning: package 'spgwr' was built under R version 4.3.3
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.3.3
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
require(splancs)
## Loading required package: splancs
## Warning: package 'splancs' was built under R version 4.3.3
## 
## 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:tidyr':
## 
##     tribble
## The following object is masked from 'package:dplyr':
## 
##     tribble
## The following object is masked from 'package:raster':
## 
##     zoom

#Membuat plot Jawa Barat

setwd("C:/Users/deni ramon siregar/Downloads/peta jabar-20240523T044012Z-001/peta jabar")
Jabar=shapefile("Jabar.shp")
## Warning in .local(x, ...): .prj file is missing
plot(Jabar)
KAB = factor(c("Kab Bogor", "Kab Sukabumi", "Kab Cianjur", "Kab Bandung", "Kab Garut",
               "Kab Tasikmalaya", "Kab Ciamis", "Kab Kuningan", "Kab Cirebon", "Kab Majalengka",
               "Kab Sumedang", "Kab Indramayu", "Kab Subang", "Kab Purwakarta", "Kab Karawang",
               "Kab Bekasi", "Kab Bandung Barat", "Kab Pangandaran", "Kota Bogor", "Kota Sukabumi",
               "Kota Bandung", "Kota Cirebon", "Kota Bekasi", "Kota Depok", "Kota Cimahi",
               "Kota Tasikmalaya", "Kota Banjar"))

Jabar$KAB = KAB
text(Jabar, "KAB", cex=0.5)

id <- c(1:27)
Jabar$id <- c(1:27)
CoordK <- coordinates(Jabar)

#Matriks Pembobotan Spasial Matriks Bobot Rook

W = poly2nb(Jabar, row.names = id, queen = F)
WB = nb2mat(W, style = "B", zero.policy = T)
WL = nb2listw(W)

Matriks Bobot Queen

W2 = poly2nb(Jabar, row.names = id, queen = T)
WB2 = nb2mat(W2, style = "B", zero.policy = T)
WL2 = nb2listw(W2)

K-Nearest (k=4)

W3 <- knn2nb(knearneigh(CoordK, k = 4), row.names = id)
WB3 <- nb2mat(W3, style='B', zero.policy = TRUE)
WL3<-nb2listw(W3)

#Memanggil data

setwd("C:/Users/deni ramon siregar/Downloads")
data_spasial <- read.csv("datauts1.csv",header = TRUE,sep = ";",dec = ",")
head(data_spasial)
##   id                kabkot     y    x1  x2  x3
## 1  1       Kabupaten Bogor 92.73 34.42 211 600
## 2  2    Kabupaten Sukabumi 92.99 19.87 130  87
## 3  3     Kabupaten Cianjur 89.78 15.70 156 431
## 4  4     Kabupaten Bandung 93.60 38.01 313 279
## 5  5       Kabupaten Garut 90.23 24.60 226 151
## 6  6 Kabupaten Tasikmalaya 89.72 24.74 117  60
y <- data_spasial$y
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   87.87   89.75   91.54   91.82   93.47   97.62
x1 <- data_spasial$x1
summary(x1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.70   24.67   31.62   35.85   47.74   66.45
x2 <- data_spasial$x2
summary(x2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     9.0    37.5    83.0   109.4   158.0   389.0
x3 <- data_spasial$x3
summary(x3)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    60.0   166.0   279.0   650.9   861.5  2991.0

#Moran test Rook

moran.test(y,WL)
## 
##  Moran I test under randomisation
## 
## data:  y  
## weights: WL    
## 
## Moran I statistic standard deviate = 3.6571, p-value = 0.0001275
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.47743307       -0.03846154        0.01989955

Queen

moran.test(y,WL2)
## 
##  Moran I test under randomisation
## 
## data:  y  
## weights: WL2    
## 
## Moran I statistic standard deviate = 3.6571, p-value = 0.0001275
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.47743307       -0.03846154        0.01989955
moran.test(y,WL3)
## 
##  Moran I test under randomisation
## 
## data:  y  
## weights: WL3    
## 
## Moran I statistic standard deviate = 3.6221, p-value = 0.0001461
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.40399702       -0.03846154        0.01492223

#Membuat plot data Plot Y

tmap_options(check.and.fix = TRUE)
par(mfrow=c(2,2))
Jabar$y.y =data_spasial$y
tm_shape(Jabar) + tm_polygons("y.y",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)
## Warning: Currect projection of shape Jabar unknown. Long-lat (WGS84) is
## assumed.
## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).

Plot X1

Jabar$x1 =data_spasial$x1
tm_shape(Jabar) + tm_polygons("x1",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)
## Warning: Currect projection of shape Jabar unknown. Long-lat (WGS84) is
## assumed.
## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).

Plot X2

Jabar$x2 =data_spasial$x2
tm_shape(Jabar) + tm_polygons("x2",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)
## Warning: Currect projection of shape Jabar unknown. Long-lat (WGS84) is
## assumed.
## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).

Plot X3

Jabar$x3 =data_spasial$x3
tm_shape(Jabar) + tm_polygons("x3",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("KAB",size=0.6)
## Warning: Currect projection of shape Jabar unknown. Long-lat (WGS84) is
## assumed.
## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).

#Model OLS

modelols=lm(y~x1+x2+x3, data=data_spasial)
summary(modelols)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = data_spasial)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9440 -0.8501 -0.0378  1.0575  3.4700 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.593e+01  1.407e+00  61.057  < 2e-16 ***
## x1          1.478e-01  3.549e-02   4.164 0.000374 ***
## x2          4.997e-03  4.339e-03   1.151 0.261365    
## x3          7.966e-05  6.227e-04   0.128 0.899309    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.883 on 23 degrees of freedom
## Multiple R-squared:  0.5618, Adjusted R-squared:  0.5047 
## F-statistic:  9.83 on 3 and 23 DF,  p-value: 0.0002309

##Uji Asumsi Asumsi Heterogenitas

bptest(modelols)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelols
## BP = 3.9529, df = 3, p-value = 0.2666

Asumsi Multikolinearitas

vif(modelols)
##       x1       x2       x3 
## 1.916636 1.248737 1.612728

#Model GWR ##Mencari Bandwith Optimal

gwr.b1 <- gwr.sel(y~x1+x2+x3, data_spasial, coords = CoordK)
## Bandwidth: 0.894596 CV score: 105.2858 
## Bandwidth: 1.446042 CV score: 104.6353 
## Bandwidth: 1.786854 CV score: 104.7619 
## Bandwidth: 1.509593 CV score: 104.6507 
## Bandwidth: 1.235408 CV score: 104.65 
## Bandwidth: 1.371466 CV score: 104.6265 
## Bandwidth: 1.366135 CV score: 104.6264 
## Bandwidth: 1.356341 CV score: 104.6263 
## Bandwidth: 1.357864 CV score: 104.6263 
## Bandwidth: 1.357931 CV score: 104.6263 
## Bandwidth: 1.357972 CV score: 104.6263 
## Bandwidth: 1.357931 CV score: 104.6263
gwr.fit1 = gwr(y~x1+x2+x3, data_spasial, coords = CoordK,
               bandwidth=gwr.b1,se.fit=T,hatmatrix=T)
gwr.fit1
## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b1, hatmatrix = T, se.fit = T)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 1.357931 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept.  8.5567e+01  8.5713e+01  8.5808e+01  8.5881e+01  8.6173e+01
## x1            1.4701e-01  1.5042e-01  1.5098e-01  1.5151e-01  1.5208e-01
## x2            4.4388e-03  4.9935e-03  5.1995e-03  5.4469e-03  5.6937e-03
## x3           -2.8650e-04 -2.8201e-05  4.3690e-05  1.1148e-04  1.5789e-04
##               Global
## X.Intercept. 85.9272
## x1            0.1478
## x2            0.0050
## x3            0.0001
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 5.603196 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 21.3968 
## Sigma (residual: 2traceS - traceS'S): 1.846307 
## Effective number of parameters (model: traceS): 4.880241 
## Effective degrees of freedom (model: traceS): 22.11976 
## Sigma (model: traceS): 1.815885 
## Sigma (ML): 1.643602 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 119.2369 
## AIC (GWR p. 96, eq. 4.22): 108.335 
## Residual sum of squares: 72.93852 
## Quasi-global R2: 0.6080315
gwr.b2 <- gwr.sel(y~x1+x2+x3, data_spasial, coords = CoordK, gweight=gwr.bisquare)
## Bandwidth: 0.894596 CV score: 263.0043 
## Bandwidth: 1.446042 CV score: 117.7826 
## Bandwidth: 1.786854 CV score: 107.4628 
## Bandwidth: 1.674408 CV score: 109.4273 
## Bandwidth: 1.88636 CV score: 106.3524 
## Bandwidth: 2.058985 CV score: 105.4788 
## Bandwidth: 2.08559 CV score: 105.4143 
## Bandwidth: 2.163834 CV score: 105.2574 
## Bandwidth: 2.230474 CV score: 105.148 
## Bandwidth: 2.271659 CV score: 105.0894 
## Bandwidth: 2.297113 CV score: 105.0563 
## Bandwidth: 2.312845 CV score: 105.0369 
## Bandwidth: 2.322568 CV score: 105.0253 
## Bandwidth: 2.328577 CV score: 105.0183 
## Bandwidth: 2.33229 CV score: 105.014 
## Bandwidth: 2.334585 CV score: 105.0114 
## Bandwidth: 2.336004 CV score: 105.0098 
## Bandwidth: 2.336881 CV score: 105.0088 
## Bandwidth: 2.337422 CV score: 105.0082 
## Bandwidth: 2.337757 CV score: 105.0078 
## Bandwidth: 2.337964 CV score: 105.0075 
## Bandwidth: 2.338092 CV score: 105.0074 
## Bandwidth: 2.338171 CV score: 105.0073 
## Bandwidth: 2.33822 CV score: 105.0072 
## Bandwidth: 2.33822 CV score: 105.0072
## Warning in gwr.sel(y ~ x1 + x2 + x3, data_spasial, coords = CoordK, gweight =
## gwr.bisquare): Bandwidth converged to upper bound:2.33829918297788
gwr.fit2 <- gwr(y~x1+x2+x3, data_spasial, coords = CoordK,
                bandwidth=gwr.b2,se.fit=T,hatmatrix=T)
gwr.fit2
## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b2, hatmatrix = T, se.fit = T)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 2.33822 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept.  8.5798e+01  8.5847e+01  8.5878e+01  8.5901e+01  8.5986e+01
## x1            1.4817e-01  1.4887e-01  1.4904e-01  1.4919e-01  1.4939e-01
## x2            4.8181e-03  5.0025e-03  5.0719e-03  5.1699e-03  5.2758e-03
## x3           -1.1100e-05  4.8746e-05  6.8289e-05  9.4222e-05  1.1673e-04
##               Global
## X.Intercept. 85.9272
## x1            0.1478
## x2            0.0050
## x3            0.0001
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 4.562993 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 22.43701 
## Sigma (residual: 2traceS - traceS'S): 1.870989 
## Effective number of parameters (model: traceS): 4.292003 
## Effective degrees of freedom (model: traceS): 22.708 
## Sigma (model: traceS): 1.859791 
## Sigma (ML): 1.705579 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 119.2534 
## AIC (GWR p. 96, eq. 4.22): 109.7455 
## Residual sum of squares: 78.54296 
## Quasi-global R2: 0.5779135
gwr.b3 <- gwr.sel(y~x1+x2+x3, data_spasial, coords = CoordK, adapt=T)
## Adaptive q: 0.381966 CV score: 105.7782 
## Adaptive q: 0.618034 CV score: 104.3396 
## Adaptive q: 0.763932 CV score: 103.8271 
## Adaptive q: 0.950804 CV score: 104.127 
## Adaptive q: 0.805191 CV score: 103.8001 
## Adaptive q: 0.8056439 CV score: 103.7989 
## Adaptive q: 0.8610901 CV score: 103.9998 
## Adaptive q: 0.8268225 CV score: 103.8508 
## Adaptive q: 0.8137334 CV score: 103.7782 
## Adaptive q: 0.818733 CV score: 103.8008 
## Adaptive q: 0.8106435 CV score: 103.7857 
## Adaptive q: 0.8156431 CV score: 103.7811 
## Adaptive q: 0.8137741 CV score: 103.7782 
## Adaptive q: 0.8143274 CV score: 103.7769 
## Adaptive q: 0.8148299 CV score: 103.7758 
## Adaptive q: 0.8151405 CV score: 103.7778 
## Adaptive q: 0.8146766 CV score: 103.7761 
## Adaptive q: 0.8149486 CV score: 103.7766 
## Adaptive q: 0.814778 CV score: 103.7758 
## Adaptive q: 0.8147373 CV score: 103.7759 
## Adaptive q: 0.814778 CV score: 103.7758
gwr.fit3 <- gwr(y~x1+x2+x3, data_spasial, coords = CoordK,
                bandwidth=gwr.b3,se.fit=T,hatmatrix=T)
gwr.fit3
## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b3, hatmatrix = T, se.fit = T)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 0.814778 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept.  8.5075e+01  8.5606e+01  8.5790e+01  8.6203e+01  8.6870e+01
## x1            1.3701e-01  1.4724e-01  1.5235e-01  1.5675e-01  1.5898e-01
## x2            3.1198e-03  4.7583e-03  5.3460e-03  5.5838e-03  6.1531e-03
## x3           -1.4934e-03 -3.1083e-04 -1.7068e-05  8.8972e-05  2.1946e-04
##               Global
## X.Intercept. 85.9272
## x1            0.1478
## x2            0.0050
## x3            0.0001
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 8.129885 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 18.87011 
## Sigma (residual: 2traceS - traceS'S): 1.756134 
## Effective number of parameters (model: traceS): 6.508086 
## Effective degrees of freedom (model: traceS): 20.49191 
## Sigma (model: traceS): 1.685208 
## Sigma (ML): 1.468125 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 119.283 
## AIC (GWR p. 96, eq. 4.22): 103.866 
## Residual sum of squares: 58.19555 
## Quasi-global R2: 0.6872596
gwr.b4 <- gwr.sel(y~x1+x2+x3, data_spasial, coords = CoordK, gweight=gwr.tricube)
## Bandwidth: 0.894596 CV score: 286.8821 
## Bandwidth: 1.446042 CV score: 121.9812 
## Bandwidth: 1.786854 CV score: 108.9068 
## Bandwidth: 1.682104 CV score: 111.2711 
## Bandwidth: 1.903186 CV score: 107.1225 
## Bandwidth: 2.069385 CV score: 105.8862 
## Bandwidth: 2.119325 CV score: 105.7082 
## Bandwidth: 2.193759 CV score: 105.5232 
## Bandwidth: 2.248969 CV score: 105.4213 
## Bandwidth: 2.28309 CV score: 105.3673 
## Bandwidth: 2.304178 CV score: 105.3366 
## Bandwidth: 2.317211 CV score: 105.3185 
## Bandwidth: 2.325266 CV score: 105.3076 
## Bandwidth: 2.330244 CV score: 105.301 
## Bandwidth: 2.333321 CV score: 105.2969 
## Bandwidth: 2.335222 CV score: 105.2944 
## Bandwidth: 2.336398 CV score: 105.2929 
## Bandwidth: 2.337124 CV score: 105.292 
## Bandwidth: 2.337573 CV score: 105.2914 
## Bandwidth: 2.33785 CV score: 105.291 
## Bandwidth: 2.338022 CV score: 105.2908 
## Bandwidth: 2.338128 CV score: 105.2907 
## Bandwidth: 2.338193 CV score: 105.2906 
## Bandwidth: 2.338234 CV score: 105.2905 
## Bandwidth: 2.338234 CV score: 105.2905
## Warning in gwr.sel(y ~ x1 + x2 + x3, data_spasial, coords = CoordK, gweight =
## gwr.tricube): Bandwidth converged to upper bound:2.33829918297788
gwr.fit4 <- gwr(y~x1+x2+x3, data_spasial, coords = CoordK,
                bandwidth=gwr.b4,se.fit=T,hatmatrix=T)
gwr.fit4
## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b4, hatmatrix = T, se.fit = T)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 2.338234 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept.  8.5798e+01  8.5847e+01  8.5878e+01  8.5901e+01  8.5986e+01
## x1            1.4817e-01  1.4887e-01  1.4904e-01  1.4919e-01  1.4939e-01
## x2            4.8181e-03  5.0025e-03  5.0719e-03  5.1699e-03  5.2758e-03
## x3           -1.1098e-05  4.8746e-05  6.8289e-05  9.4222e-05  1.1673e-04
##               Global
## X.Intercept. 85.9272
## x1            0.1478
## x2            0.0050
## x3            0.0001
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 4.562986 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 22.43701 
## Sigma (residual: 2traceS - traceS'S): 1.870989 
## Effective number of parameters (model: traceS): 4.291999 
## Effective degrees of freedom (model: traceS): 22.708 
## Sigma (model: traceS): 1.859792 
## Sigma (ML): 1.705579 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 119.2534 
## AIC (GWR p. 96, eq. 4.22): 109.7455 
## Residual sum of squares: 78.54299 
## Quasi-global R2: 0.5779133
gwr.b5 <- gwr.sel(y~x1+x2+x3, data_spasial, coords = CoordK, gweight=gwr.bisquare, adapt=T)
## Adaptive q: 0.381966 CV score: 323.3525 
## Adaptive q: 0.618034 CV score: 138.852 
## Adaptive q: 0.763932 CV score: 109.936 
## Adaptive q: 0.755868 CV score: 110.5396 
## Adaptive q: 0.8017185 CV score: 107.4483 
## Adaptive q: 0.8774553 CV score: 106.3327 
## Adaptive q: 0.8559472 CV score: 107.0191 
## Adaptive q: 0.9242632 CV score: 106.0968 
## Adaptive q: 0.9072666 CV score: 106.0649 
## Adaptive q: 0.9117225 CV score: 106.0725 
## Adaptive q: 0.8958797 CV score: 106.0481 
## Adaptive q: 0.8888422 CV score: 106.0404 
## Adaptive q: 0.8844928 CV score: 106.147 
## Adaptive q: 0.8915303 CV score: 106.0425 
## Adaptive q: 0.8871809 CV score: 106.0804 
## Adaptive q: 0.889869 CV score: 106.0404 
## Adaptive q: 0.8892824 CV score: 106.0397 
## Adaptive q: 0.8893434 CV score: 106.0398 
## Adaptive q: 0.8891997 CV score: 106.0396 
## Adaptive q: 0.8890632 CV score: 106.0395 
## Adaptive q: 0.8889788 CV score: 106.0394 
## Adaptive q: 0.8889266 CV score: 106.0393 
## Adaptive q: 0.8888859 CV score: 106.0393 
## Adaptive q: 0.8889266 CV score: 106.0393
gwr.fit5 <- gwr(y~x1+x2+x3, data_spasial, coords = CoordK,
                bandwidth=gwr.b5,se.fit=T,hatmatrix=T)
gwr.fit5
## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b5, hatmatrix = T, se.fit = T)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 0.8889266 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept.  8.5184e+01  8.5631e+01  8.5776e+01  8.6119e+01  8.6692e+01
## x1            1.3997e-01  1.4910e-01  1.5230e-01  1.5572e-01  1.5729e-01
## x2            3.4764e-03  4.8238e-03  5.3096e-03  5.5706e-03  6.0913e-03
## x3           -1.1709e-03 -2.2991e-04 -2.8361e-06  1.0027e-04  1.9283e-04
##               Global
## X.Intercept. 85.9272
## x1            0.1478
## x2            0.0050
## x3            0.0001
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 7.531969 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 19.46803 
## Sigma (residual: 2traceS - traceS'S): 1.781439 
## Effective number of parameters (model: traceS): 6.100115 
## Effective degrees of freedom (model: traceS): 20.89989 
## Sigma (model: traceS): 1.719333 
## Sigma (ML): 1.51269 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 119.2589 
## AIC (GWR p. 96, eq. 4.22): 105.0728 
## Residual sum of squares: 61.78226 
## Quasi-global R2: 0.6679848
gwr.b6 <- gwr.sel(y~x1+x2+x3, data_spasial, coords = CoordK, gweight=gwr.tricube, adapt=T)
## Adaptive q: 0.381966 CV score: 395.6541 
## Adaptive q: 0.618034 CV score: 149.3318 
## Adaptive q: 0.763932 CV score: 112.3934 
## Adaptive q: 0.7521692 CV score: 113.4382 
## Adaptive q: 0.7942934 CV score: 109.8069 
## Adaptive q: 0.8728663 CV score: 107.7016 
## Adaptive q: 0.858572 CV score: 108.3562 
## Adaptive q: 0.9214271 CV score: 107.1975 
## Adaptive q: 0.9063575 CV score: 107.159 
## Adaptive q: 0.9105842 CV score: 107.17 
## Adaptive q: 0.893565 CV score: 107.1244 
## Adaptive q: 0.8856588 CV score: 107.2198 
## Adaptive q: 0.8984513 CV score: 107.1378 
## Adaptive q: 0.8905451 CV score: 107.1159 
## Adaptive q: 0.8886787 CV score: 107.1182 
## Adaptive q: 0.8903427 CV score: 107.1154 
## Adaptive q: 0.8898597 CV score: 107.114 
## Adaptive q: 0.8894086 CV score: 107.1127 
## Adaptive q: 0.8891298 CV score: 107.1119 
## Adaptive q: 0.8889575 CV score: 107.1114 
## Adaptive q: 0.888851 CV score: 107.1125 
## Adaptive q: 0.8890125 CV score: 107.1116 
## Adaptive q: 0.8889168 CV score: 107.1113 
## Adaptive q: 0.8889168 CV score: 107.1113
gwr.fit6 <- gwr(y~x1+x2+x3, data_spasial, coords = CoordK,
                bandwidth=gwr.b6,se.fit=T,hatmatrix=T)
gwr.fit6
## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b6, hatmatrix = T, se.fit = T)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 0.8889168 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept.  8.5184e+01  8.5631e+01  8.5776e+01  8.6119e+01  8.6692e+01
## x1            1.3996e-01  1.4910e-01  1.5230e-01  1.5572e-01  1.5729e-01
## x2            3.4763e-03  4.8238e-03  5.3096e-03  5.5706e-03  6.0913e-03
## x3           -1.1709e-03 -2.2991e-04 -2.8378e-06  1.0027e-04  1.9283e-04
##               Global
## X.Intercept. 85.9272
## x1            0.1478
## x2            0.0050
## x3            0.0001
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 7.53204 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 19.46796 
## Sigma (residual: 2traceS - traceS'S): 1.781436 
## Effective number of parameters (model: traceS): 6.100162 
## Effective degrees of freedom (model: traceS): 20.89984 
## Sigma (model: traceS): 1.719329 
## Sigma (ML): 1.512685 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 119.2589 
## AIC (GWR p. 96, eq. 4.22): 105.0727 
## Residual sum of squares: 61.78183 
## Quasi-global R2: 0.667987

Perbandingan AIC

AIC(modelols)
## [1] 116.4639
gwr.fit1$results$AICh
## [1] 108.335
gwr.fit2$results$AICh
## [1] 109.7455
gwr.fit3$results$AICh
## [1] 103.866
gwr.fit4$results$AICh
## [1] 109.7455
gwr.fit5$results$AICh
## [1] 105.0728
gwr.fit6$results$AICh
## [1] 105.0727

##Uji Asumsi Model GWR Uji Autokorelasi Spasial

residualgwr=gwr.fit3$SDF$gwr.e
coords <- CoordK
nb <- knn2nb(knearneigh(coords, k=4))
listw <- nb2listw(nb, style="W")
moran.test(residualgwr, listw)
## 
##  Moran I test under randomisation
## 
## data:  residualgwr  
## weights: listw    
## 
## Moran I statistic standard deviate = 0.61068, p-value = 0.2707
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.03561950       -0.03846154        0.01471569

Uji Normalitas

shapiro.test(residualgwr)
## 
##  Shapiro-Wilk normality test
## 
## data:  residualgwr
## W = 0.97057, p-value = 0.6169

#Uji Signifikansi Lokal

t1 = gwr.fit3$SDF$x1 / gwr.fit3$SDF$x1_se
t2 = gwr.fit3$SDF$x2 / gwr.fit3$SDF$x2_se
t3 = gwr.fit3$SDF$x3 / gwr.fit3$SDF$x3_se
pvalue1 = 2*pt(q=t1, df=26, lower.tail=FALSE)
pvalue2 = 2*pt(q=t2, df=26, lower.tail=FALSE)
pvalue3 = 2*pt(q=t3, df=26, lower.tail=FALSE)
data.frame(id, t1, pvalue1) %>% mutate(ket=ifelse(pvalue1<0.05,"Reject", "Not
Reject"))
##    id       t1      pvalue1    ket
## 1   1 4.044870 4.158440e-04 Reject
## 2   2 3.848315 6.935204e-04 Reject
## 3   3 4.383862 1.708536e-04 Reject
## 4   4 4.752444 6.460261e-05 Reject
## 5   5 4.704976 7.323118e-05 Reject
## 6   6 4.502716 1.249025e-04 Reject
## 7   7 4.135127 3.284031e-04 Reject
## 8   8 3.894880 6.146029e-04 Reject
## 9   9 3.926792 5.656927e-04 Reject
## 10 10 4.501066 1.254473e-04 Reject
## 11 11 4.793784 5.792041e-05 Reject
## 12 12 4.631022 8.902455e-05 Reject
## 13 13 4.858636 4.880368e-05 Reject
## 14 14 4.735257 6.760269e-05 Reject
## 15 15 4.600554 9.648139e-05 Reject
## 16 16 4.364303 1.798831e-04 Reject
## 17 17 4.694058 7.537347e-05 Reject
## 18 18 3.903863 6.004302e-04 Reject
## 19 19 4.084447 3.749820e-04 Reject
## 20 20 4.179555 2.923073e-04 Reject
## 21 21 4.811061 5.533697e-05 Reject
## 22 22 3.921832 5.730347e-04 Reject
## 23 23 4.234348 2.531605e-04 Reject
## 24 24 4.083757 3.756592e-04 Reject
## 25 25 4.777627 6.044556e-05 Reject
## 26 26 4.452704 1.425099e-04 Reject
## 27 27 3.838546 7.112978e-04 Reject
data.frame(id, t2, pvalue2) %>% mutate(ket=ifelse(pvalue2<0.05,"Reject", "Not
Reject"))
##    id        t2   pvalue2         ket
## 1   1 0.9954712 0.3286821 Not\nReject
## 2   2 0.9356394 0.3580685 Not\nReject
## 3   3 1.1047910 0.2793696 Not\nReject
## 4   4 1.2549310 0.2206672 Not\nReject
## 5   5 1.1567067 0.2579098 Not\nReject
## 6   6 0.9818888 0.3352038 Not\nReject
## 7   7 0.9625735 0.3446292 Not\nReject
## 8   8 1.0591122 0.2992894 Not\nReject
## 9   9 1.2067427 0.2383965 Not\nReject
## 10 10 1.3201249 0.1982983 Not\nReject
## 11 11 1.3781194 0.1799119 Not\nReject
## 12 12 1.4687636 0.1538930 Not\nReject
## 13 13 1.4372948 0.1625609 Not\nReject
## 14 14 1.3125159 0.2008150 Not\nReject
## 15 15 1.2939917 0.2070451 Not\nReject
## 16 16 1.1615904 0.2559550 Not\nReject
## 17 17 1.2547663 0.2207260 Not\nReject
## 18 18 0.7068122 0.4859701 Not\nReject
## 19 19 1.0097571 0.3219171 Not\nReject
## 20 20 1.0405321 0.3076726 Not\nReject
## 21 21 1.3172129 0.1992585 Not\nReject
## 22 22 1.2071901 0.2382271 Not\nReject
## 23 23 1.0859406 0.2874715 Not\nReject
## 24 24 1.0124960 0.3206312 Not\nReject
## 25 25 1.3008812 0.2047109 Not\nReject
## 26 26 1.0296538 0.3126565 Not\nReject
## 27 27 0.8251889 0.4167716 Not\nReject
data.frame(id, t3, pvalue3) %>% mutate(ket=ifelse(pvalue3<0.05,"Reject", "Not
Reject"))
##    id          t3   pvalue3         ket
## 1   1  0.13973709 0.8899446 Not\nReject
## 2   2  0.08613946 0.9320155 Not\nReject
## 3   3 -0.03259587 1.0257542 Not\nReject
## 4   4 -0.21630234 1.1695598 Not\nReject
## 5   5 -0.73064182 1.5284688 Not\nReject
## 6   6 -1.28593284 1.7901986 Not\nReject
## 7   7 -1.16870844 1.7468745 Not\nReject
## 8   8 -0.76465353 1.5486373 Not\nReject
## 9   9 -0.27392684 1.2136946 Not\nReject
## 10 10 -0.20287958 1.1591915 Not\nReject
## 11 11 -0.06006080 1.0474335 Not\nReject
## 12 12  0.38197756 0.7055823 Not\nReject
## 13 13  0.32775654 0.7457204 Not\nReject
## 14 14  0.22383601 0.8246344 Not\nReject
## 15 15  0.31282357 0.7569102 Not\nReject
## 16 16  0.23405962 0.8167720 Not\nReject
## 17 17  0.05080784 0.9598670 Not\nReject
## 18 18 -1.66917724 1.8929187 Not\nReject
## 19 19  0.14236772 0.8878871 Not\nReject
## 20 20  0.09381129 0.9259783 Not\nReject
## 21 21 -0.03008113 1.0237679 Not\nReject
## 22 22 -0.26855726 1.2096101 Not\nReject
## 23 23  0.18664078 0.8533910 Not\nReject
## 24 24  0.14835271 0.8832090 Not\nReject
## 25 25  0.02939939 0.9767706 Not\nReject
## 26 26 -1.13562175 1.7335248 Not\nReject
## 27 27 -1.37018014 1.8176532 Not\nReject

#Plot Koefisien GWR X1

spplot(gwr.fit3$SDF,"x1", main = "Koefisien GWR untuk x1", col.regions = terrain.colors(100))

X2

spplot(gwr.fit3$SDF,"x2", main = "Koefisien GWR untuk x2", col.regions = terrain.colors(100))

X3

spplot(gwr.fit3$SDF,"x3", main = "Koefisien GWR untuk x3", col.regions = terrain.colors(100))

#Peta Sebaran Penduga Parameter GWR di Setiap Lokasi

Jabar$bwadapt = gwr.fit3$bandwidth
labels <- c("Very Low", "Low", "High", "Very High")
jb_sf<-st_as_sf(Jabar)
jb_merged <- jb_sf %>% left_join(data_spasial, by = "id")

##Penduga Parameter beta1(X1)

b.x1<- gwr.fit3$SDF$x1
Jabar@data$b1<- b.x1
summary(b.x1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1370  0.1472  0.1524  0.1517  0.1567  0.1590
breaks_b.x1 <- c(-Inf, 0.1472, 0.1524, 0.1567,Inf)

jb_merged$b.x1_disc <-
  cut(b.x1, breaks = breaks_b.x1, labels = labels, right = TRUE)

ggplot() +
  geom_sf(data=jb_merged, aes(fill = b.x1_disc),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(title = "Peta Sebaran Penduga Parameter beta1 (x1)",
       fill="b.x1")

##Penduga Parameter beta2(X2)

b.x2<- gwr.fit3$SDF$x2
Jabar@data$b2<- b.x2
summary(b.x2)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.003120 0.004758 0.005346 0.005103 0.005584 0.006153
breaks_b.x2 <- c(-Inf, 0.004758, 0.005346, 0.005584,Inf)
jb_merged$b.x2_disc <-
  cut(b.x2, breaks = breaks_b.x2, labels = labels, right = TRUE)

ggplot() +
  geom_sf(data=jb_merged, aes(fill = b.x2_disc),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(title = "Peta Sebaran Penduga Parameter beta2 (x2)",
       fill="b.x2")

##Penduga Parameter beta3(X3)

b.x3<- gwr.fit3$SDF$x3
Jabar@data$b3<- b.x3
summary(b.x3)
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.493e-03 -3.108e-04 -1.707e-05 -2.025e-04  8.897e-05  2.195e-04
breaks_b.x3 <- c(-Inf, -3.108e-04, -1.707e-05, 8.897e-05,Inf)

jb_merged$b.x3_disc <-
  cut(b.x3, breaks = breaks_b.x3, labels = labels, right = TRUE)

ggplot() +
  geom_sf(data=jb_merged, aes(fill = b.x3_disc),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(title = "Peta Sebaran Penduga Parameter beta3 (x3)",
       fill="b.x3")

#Plot Nilai Aktual Y

summary(data_spasial$y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   87.87   89.75   91.54   91.82   93.47   97.62
breaks <- c(-Inf, 89.75, 91.54, 93.47,Inf)
jb_merged$y_Discrete <-
  cut(jb_merged$y.y, breaks =
        breaks, labels = labels, right = TRUE)

ggplot() +
  geom_sf(data=jb_merged, aes(fill = y_Discrete),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(title = "Peta Sebaran Y",
       fill = "Y")

#Plot Sebaran Dugaan Y (GWR)

gwr.prediksi <- gwr.fit3$SDF$pred
Jabar@data$pred <- gwr.prediksi
summary(gwr.prediksi)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   89.33   90.20   90.92   91.82   93.41   96.79
breaks_gwrpred <- c(-Inf, 90.20, 90.92, 93.41,Inf)

jb_merged$gwrpred_disc <-
  cut(gwr.prediksi, breaks = breaks_gwrpred, labels = labels, right = TRUE)

ggplot() +
  geom_sf(data=jb_merged, aes(fill = gwrpred_disc),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(title = "Peta Sebaran Prediksi Y Model GWR",
       fill="gwrpred")

#Plot Sebaran Dugaan Y (OLS)

modelols.pred <- predict(modelols,data=data_spasial)
Jabar@data$rlbpred <- modelols.pred
summary(modelols.pred)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   89.06   90.31   91.19   91.82   93.27   96.29
breaks_olspred <- c(-Inf, 90.31, 91.19, 93.27,Inf)

jb_merged$ols_disc <-
  cut(modelols.pred, breaks =
        breaks_olspred, labels = labels, right = TRUE)

ggplot() +
  geom_sf(data=jb_merged, aes(fill = ols_disc),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(title = "Peta Sebaran Prediksi Y Model OLS",
       fill = "Pred Y OLS")

#Plot Sisaan Model GWR

rtg.sisaan <- data_spasial$y-gwr.fit3$SDF$pred
Jabar@data$rtg.sisaan <- rtg.sisaan
summary(rtg.sisaan)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -2.810694 -1.054542  0.053079  0.003024  0.828803  2.840422
breaks_rtg <- c(-Inf, -1.054542, 0.053079, 0.828803,Inf)

jb_merged$rtg_disc <-
  cut(rtg.sisaan, breaks =
        breaks_rtg, labels = labels, right = TRUE)

ggplot() +
  geom_sf(data=jb_merged, aes(fill = rtg_disc),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(title = "Peta Sebaran Sisaan Model GWR",
       fill = "sisaan")