#Mapping dan data
setwd("/Users/rayfarekafitri/Desktop/SPASIAL")
Jabar=shapefile("Jabar.shp")
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)

data_spasial <- read.csv("Jabar 2022.csv",sep = ";", dec = ",")
head(data_spasial)
##   id    Kab.Kota     UMK  PDRB   IPM
## 1  1       Bogor 4217206 48096 71.20
## 2  2    Sukabumi 3125445 27165 67.64
## 3  3     Cianjur 2699814 21232 65.94
## 4  4     Bandung 3241930 38455 73.16
## 5  5       Garut 1975221 25346 67.41
## 6  6 Tasikmalaya 2326772 22378 66.84
names(data_spasial)
## [1] "id"       "Kab.Kota" "UMK"      "PDRB"     "IPM"
#statistik deskriptif
deskriptif <- data_spasial %>%
  summarize(
    Min_Y = min(UMK, na.rm = TRUE),
    Max_Y = max(UMK, na.rm = TRUE),
    Mean_Y = mean(UMK, na.rm = TRUE),
    StdDev_Y = sd(UMK, na.rm = TRUE),
    Variance_Y = var(UMK, na.rm = TRUE),
    
    Min_X1 = min(IPM, na.rm = TRUE),
    Max_X1 = max(IPM, na.rm = TRUE),
    Mean_X1 = mean(IPM, na.rm = TRUE),
    StdDev_X1 = sd(IPM, na.rm = TRUE),
    Variance_X1 = var(IPM, na.rm = TRUE),
    
    Min_X2 = min(PDRB, na.rm = TRUE),
    Max_X2 = max(PDRB, na.rm = TRUE),
    Mean_X2 = mean(PDRB, na.rm = TRUE),
    StdDev_X2 = sd(PDRB, na.rm = TRUE),
    Variance_X2 = var(PDRB, na.rm = TRUE),
  )
print(deskriptif)
##     Min_Y   Max_Y  Mean_Y StdDev_Y   Variance_Y Min_X1 Max_X1  Mean_X1
## 1 1852100 4816921 3072179  1009197 1.018478e+12  65.94   82.5 72.61296
##   StdDev_X1 Variance_X1 Min_X2 Max_X2  Mean_X2 StdDev_X2 Variance_X2
## 1  4.720933    22.28721  21232 133378 46685.11  30059.02   903544662
#MATRIKS BOBOT
## Rook
W <- poly2nb(Jabar, row.names=id, queen=FALSE)
WB <- nb2mat(W, style='B', zero.policy = TRUE)
WB[c(1:5),c(1:5)]
##   [,1] [,2] [,3] [,4] [,5]
## 1    0    1    1    0    0
## 2    1    0    1    0    0
## 3    1    1    0    1    1
## 4    0    0    1    0    1
## 5    0    0    1    1    0
WL<-nb2listw(W)
WL
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 106 
## Percentage nonzero weights: 14.54047 
## Average number of links: 3.925926 
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 27 729 27 16.73552 121.5371
moran.test(data_spasial$UMK, WL)
## 
##  Moran I test under randomisation
## 
## data:  data_spasial$UMK  
## weights: WL    
## 
## Moran I statistic standard deviate = 5.7843, p-value = 3.64e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.78440760       -0.03846154        0.02023736
## Queen
W2 <- poly2nb(Jabar, row.names=id, queen=TRUE) #Mendapatkan W
WB2 <- nb2mat(W2, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"
WB2[1:5,1:5]
##   [,1] [,2] [,3] [,4] [,5]
## 1    0    1    1    0    0
## 2    1    0    1    0    0
## 3    1    1    0    1    1
## 4    0    0    1    0    1
## 5    0    0    1    1    0
WL2 <-nb2listw(W2) ;WL2
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 106 
## Percentage nonzero weights: 14.54047 
## Average number of links: 3.925926 
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 27 729 27 16.73552 121.5371
moran.test(data_spasial$UMK, WL2)
## 
##  Moran I test under randomisation
## 
## data:  data_spasial$UMK  
## weights: WL2    
## 
## Moran I statistic standard deviate = 5.7843, p-value = 3.64e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.78440760       -0.03846154        0.02023736
## K-nearest
W3 <- knn2nb(knearneigh(CoordK, k = 5), row.names = id) ;W
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 106 
## Percentage nonzero weights: 14.54047 
## Average number of links: 3.925926
WB3 <- nb2mat(W3, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"
WB3[1:5,1:5]
##   [,1] [,2] [,3] [,4] [,5]
## 1    0    0    0    0    0
## 2    1    0    1    0    0
## 3    0    1    0    1    0
## 4    0    0    1    0    1
## 5    0    0    0    1    0
WL3<-nb2listw(W3)
moran.test(data_spasial$UMK, WL3)
## 
##  Moran I test under randomisation
## 
## data:  data_spasial$UMK  
## weights: WL3    
## 
## Moran I statistic standard deviate = 7.4313, p-value = 5.377e-14
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.76428955       -0.03846154        0.01166897
summary(data_spasial$UMK)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 1852100 2292463 3064218 3072179 3974215 4816921
y = data_spasial$UMK
x1 = data_spasial$IPM
x2 = data_spasial$PDRB
#Bandwidth optimal
bandwidth_optimal <- gwr.sel(
  y ~ x1 + x2,
  data = data_spasial,
  coords = CoordK,
  adapt = TRUE
)
## Adaptive q: 0.381966 CV score: 1.297736e+13 
## Adaptive q: 0.618034 CV score: 1.607583e+13 
## Adaptive q: 0.236068 CV score: 1.002845e+13 
## Adaptive q: 0.145898 CV score: 7.790378e+12 
## Adaptive q: 0.09016994 CV score: 7.297545e+12 
## Adaptive q: 0.07765572 CV score: 7.401549e+12 
## Adaptive q: 0.1004436 CV score: 7.128993e+12 
## Adaptive q: 0.1178057 CV score: 7.057578e+12 
## Adaptive q: 0.1137482 CV score: 7.036849e+12 
## Adaptive q: 0.1120917 CV score: 7.042434e+12 
## Adaptive q: 0.1140558 CV score: 7.036797e+12 
## Adaptive q: 0.1139544 CV score: 7.036781e+12 
## Adaptive q: 0.1139951 CV score: 7.036784e+12 
## Adaptive q: 0.1139137 CV score: 7.036784e+12 
## Adaptive q: 0.1139544 CV score: 7.036781e+12
#model GWR
model_gwr <- gwr(
  y ~ x1 + x2,
  data = data_spasial,
  coords = CoordK,
  adapt = bandwidth_optimal,
  hatmatrix = TRUE,
  se.fit = TRUE
)
#koef pada setiap lokasi
coef_gwr <- model_gwr$SDF@data[, c("x1", "x2")]
summary(model_gwr)
##           Length Class                  Mode     
## SDF        27    SpatialPointsDataFrame S4       
## lhat      729    -none-                 numeric  
## lm         11    -none-                 list     
## results    14    -none-                 list     
## bandwidth  27    -none-                 numeric  
## adapt       1    -none-                 numeric  
## hatmatrix   1    -none-                 logical  
## gweight     1    -none-                 character
## gTSS        1    -none-                 numeric  
## this.call   7    -none-                 call     
## fp.given    1    -none-                 logical  
## timings    12    -none-                 numeric
#model ols
model_ols <- lm(y ~ x1 + x2, data = data_spasial)
summary(model_ols)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = data_spasial)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1436659  -616065     4109   680379  1225993 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -2.748e+06  2.675e+06  -1.027   0.3144  
## x1           7.121e+04  3.865e+04   1.843   0.0778 .
## x2           1.392e+01  6.070e+00   2.294   0.0308 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 795800 on 24 degrees of freedom
## Multiple R-squared:  0.426,  Adjusted R-squared:  0.3782 
## F-statistic: 8.907 on 2 and 24 DF,  p-value: 0.001278
##UJI ASUMSI
#hetergenitas
bptest(model_ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 8.4484, df = 2, p-value = 0.01464
#AIC model ols
AIC(model_ols)
## [1] 815.1453
#multikolinearitas
vif(model_ols)
##       x1       x2 
## 1.366582 1.366582
#model GWR
#Gaussian
gwr.b1 = gwr.sel(y ~ x1 + x2, data = data_spasial, coords = CoordK) # GWR
## Bandwidth: 0.894596 CV score: 1.40175e+13 
## Bandwidth: 1.446042 CV score: 1.797911e+13 
## Bandwidth: 0.5537839 CV score: 1.015825e+13 
## Bandwidth: 0.3431504 CV score: 7.710901e+12 
## Bandwidth: 0.2129718 CV score: 9.029513e+12 
## Bandwidth: 0.3574277 CV score: 7.836414e+12 
## Bandwidth: 0.3167292 CV score: 7.565956e+12 
## Bandwidth: 0.2770974 CV score: 7.736796e+12 
## Bandwidth: 0.3114455 CV score: 7.556492e+12 
## Bandwidth: 0.306402 CV score: 7.555284e+12 
## Bandwidth: 0.3081271 CV score: 7.554788e+12 
## Bandwidth: 0.3081698 CV score: 7.554788e+12 
## Bandwidth: 0.3082105 CV score: 7.554789e+12 
## Bandwidth: 0.3081698 CV score: 7.554788e+12
gwr.b1
## [1] 0.3081698
gwr.fit1 = gwr(y ~ x1 + x2, data = data_spasial,coords = CoordK,
               bandwidth=gwr.b1,se.fit=TRUE,hatmatrix=TRUE)
gwr.fit1
## Call:
## gwr(formula = y ~ x1 + x2, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b1, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 0.3081698 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -3.4392e+06 -1.6331e+06  8.1718e+05  2.5629e+06  4.6457e+06
## x1           -2.9136e+04 -7.6156e+03  1.8199e+04  6.7734e+04  8.2009e+04
## x2            3.5582e-01  6.3895e+00  1.1213e+01  1.3831e+01  2.8097e+01
##                   Global
## X.Intercept. -2.7483e+06
## x1            7.1206e+04
## x2            1.3924e+01
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 16.12423 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 10.87577 
## Sigma (residual: 2traceS - traceS'S): 404625.4 
## Effective number of parameters (model: traceS): 12.84091 
## Effective degrees of freedom (model: traceS): 14.15909 
## Sigma (model: traceS): 354621.9 
## Sigma (ML): 256803.8 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 810.7195 
## AIC (GWR p. 96, eq. 4.22): 762.0912 
## Residual sum of squares: 1.780601e+12 
## Quasi-global R2: 0.9327579
#bisquare
gwr.b2 <- gwr.sel(y ~ x1 + x2, data = data_spasial, coords = CoordK, gweight=gwr.bisquare)
## Bandwidth: 0.894596 CV score: 9.856972e+12 
## Bandwidth: 1.446042 CV score: 1.105687e+13 
## Bandwidth: 0.5537839 CV score: 9.456118e+12 
## Bandwidth: 0.1993288 CV score: NA
## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value
## Bandwidth: 0.4183941 CV score: NA
## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value
## Bandwidth: 0.6839625 CV score: 1.121812e+13 
## Bandwidth: 0.5020696 CV score: 1.264932e+13 
## Bandwidth: 0.6035077 CV score: 9.754935e+12 
## Bandwidth: 0.5340308 CV score: 9.446588e+12 
## Bandwidth: 0.5408752 CV score: 9.44596e+12 
## Bandwidth: 0.5384846 CV score: 9.445079e+12 
## Bandwidth: 0.5378968 CV score: 9.444995e+12 
## Bandwidth: 0.5372582 CV score: 9.444981e+12 
## Bandwidth: 0.537464 CV score: 9.444977e+12 
## Bandwidth: 0.5375047 CV score: 9.444977e+12 
## Bandwidth: 0.5374233 CV score: 9.444977e+12 
## Bandwidth: 0.537464 CV score: 9.444977e+12
gwr.fit2 <- gwr(y ~ x1 + x2, data = data_spasial, coords = CoordK ,bandwidth=gwr.b2, se.fit=TRUE, hatmatrix=TRUE)
gwr.fit2
## Call:
## gwr(formula = y ~ x1 + x2, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b2, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 0.537464 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -1.8107e+06 -8.3545e+05  4.0744e+05  9.9939e+05  1.5231e+06
## x1            2.7614e+03  1.2168e+04  2.8876e+04  5.3006e+04  6.8198e+04
## x2            9.9004e+00  1.1681e+01  1.1945e+01  1.3593e+01  1.4869e+01
##                   Global
## X.Intercept. -2.7483e+06
## x1            7.1206e+04
## x2            1.3924e+01
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 8.868996 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 18.131 
## Sigma (residual: 2traceS - traceS'S): 500464.2 
## Effective number of parameters (model: traceS): 6.844831 
## Effective degrees of freedom (model: traceS): 20.15517 
## Sigma (model: traceS): 474668.9 
## Sigma (ML): 410111.6 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 797.862 
## AIC (GWR p. 96, eq. 4.22): 781.3735 
## Residual sum of squares: 4.541172e+12 
## Quasi-global R2: 0.8285084
#adaptive
gwr.b3 <- gwr.sel(y ~ x1 + x2, data = data_spasial, coords = CoordK, adapt = TRUE)
## Adaptive q: 0.381966 CV score: 1.297736e+13 
## Adaptive q: 0.618034 CV score: 1.607583e+13 
## Adaptive q: 0.236068 CV score: 1.002845e+13 
## Adaptive q: 0.145898 CV score: 7.790378e+12 
## Adaptive q: 0.09016994 CV score: 7.297545e+12 
## Adaptive q: 0.07765572 CV score: 7.401549e+12 
## Adaptive q: 0.1004436 CV score: 7.128993e+12 
## Adaptive q: 0.1178057 CV score: 7.057578e+12 
## Adaptive q: 0.1137482 CV score: 7.036849e+12 
## Adaptive q: 0.1120917 CV score: 7.042434e+12 
## Adaptive q: 0.1140558 CV score: 7.036797e+12 
## Adaptive q: 0.1139544 CV score: 7.036781e+12 
## Adaptive q: 0.1139951 CV score: 7.036784e+12 
## Adaptive q: 0.1139137 CV score: 7.036784e+12 
## Adaptive q: 0.1139544 CV score: 7.036781e+12
gwr.fit3 <- gwr(y ~ x1 + x2, data = data_spasial, coords = CoordK,adapt=gwr.b3,se.fit=TRUE,hatmatrix=TRUE)
gwr.fit3
## Call:
## gwr(formula = y ~ x1 + x2, data = data_spasial, coords = CoordK, 
##     adapt = gwr.b3, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Adaptive quantile: 0.1139544 (about 3 of 27 data points)
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -3.3989e+06 -4.3354e+05  7.5355e+05  1.6895e+06  6.2428e+06
## x1           -4.8906e+04  2.8261e+03  1.9889e+04  4.7325e+04  8.1026e+04
## x2            4.7358e-01  9.8356e+00  1.1549e+01  1.3019e+01  1.5155e+01
##                   Global
## X.Intercept. -2.7483e+06
## x1            7.1206e+04
## x2            1.3924e+01
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 13.88439 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 13.11561 
## Sigma (residual: 2traceS - traceS'S): 430160.1 
## Effective number of parameters (model: traceS): 11.02541 
## Effective degrees of freedom (model: traceS): 15.97459 
## Sigma (model: traceS): 389771 
## Sigma (ML): 299807.6 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 804.0791 
## AIC (GWR p. 96, eq. 4.22): 768.6365 
## Residual sum of squares: 2.426883e+12 
## Quasi-global R2: 0.9083518
#tricube
gwr.b4 = gwr.sel(y ~ x1 + x2, data = data_spasial, coords = CoordK, gweight = gwr.tricube)
## Bandwidth: 0.894596 CV score: 1.038942e+13 
## Bandwidth: 1.446042 CV score: 1.159096e+13 
## Bandwidth: 0.5537839 CV score: 9.528785e+12 
## Bandwidth: 0.3431504 CV score: NA
## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value
## Bandwidth: 0.6839625 CV score: 1.116352e+13 
## Bandwidth: 0.473329 CV score: NA
## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value
## Bandwidth: 0.6035077 CV score: 9.68769e+12 
## Bandwidth: 0.5230529 CV score: 1.024866e+13 
## Bandwidth: 0.5738166 CV score: 9.594639e+12 
## Bandwidth: 0.5420457 CV score: 9.506025e+12 
## Bandwidth: 0.5347911 CV score: 9.563155e+12 
## Bandwidth: 0.5460385 CV score: 9.506599e+12 
## Bandwidth: 0.5437322 CV score: 9.504691e+12 
## Bandwidth: 0.5438648 CV score: 9.504697e+12 
## Bandwidth: 0.5437729 CV score: 9.504691e+12 
## Bandwidth: 0.5436915 CV score: 9.504692e+12 
## Bandwidth: 0.5437322 CV score: 9.504691e+12
gwr.fit4 = gwr(y ~ x1 + x2, data = data_spasial,coords = CoordK,bandwidth=gwr.b4,se.fit=TRUE,hatmatrix=TRUE)
gwr.fit4
## Call:
## gwr(formula = y ~ x1 + x2, data = data_spasial, coords = CoordK, 
##     bandwidth = gwr.b4, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 0.5437322 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -1.8198e+06 -8.5546e+05  3.4288e+05  1.0092e+06  1.5085e+06
## x1            3.0016e+03  1.2055e+04  2.9751e+04  5.3273e+04  6.8188e+04
## x2            9.9396e+00  1.1686e+01  1.2116e+01  1.3632e+01  1.4905e+01
##                   Global
## X.Intercept. -2.7483e+06
## x1            7.1206e+04
## x2            1.3924e+01
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 8.755985 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 18.24401 
## Sigma (residual: 2traceS - traceS'S): 502645.8 
## Effective number of parameters (model: traceS): 6.761141 
## Effective degrees of freedom (model: traceS): 20.23886 
## Sigma (model: traceS): 477231.7 
## Sigma (ML): 413181.1 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 797.9098 
## AIC (GWR p. 96, eq. 4.22): 781.6924 
## Residual sum of squares: 4.609402e+12 
## Quasi-global R2: 0.8259318
#model residual GWR
residual_gwr <- model_gwr$SDF$gwr.e

#membuat matriks bobot spasial berdasarkan tetangga terdekat
coords <- CoordK
nb <- knn2nb(knearneigh(coords, k=4))
listw <- nb2listw(nb, style="W")

#moran residual GWR
moran.test(residual_gwr, listw)
## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: listw    
## 
## Moran I statistic standard deviate = -0.81026, p-value = 0.7911
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.13116466       -0.03846154        0.01308985
#normalitas residual GWR
shapiro.test(residual_gwr)
## 
##  Shapiro-Wilk normality test
## 
## data:  residual_gwr
## W = 0.93313, p-value = 0.08249
##PLOT SPPLOT
#koefisien lokal untuk variabel x1 dan x2
coef_x1 <- gwr.fit1$SDF$x1
coef_x2 <- gwr.fit1$SDF$x2
#plot koefisien regresi lokal untuk variabel
#X1
spplot(gwr.fit1$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))

#Konversi SpatialPointsDataFrame menjadi data frame biasa
gwr_df <- as.data.frame(gwr.fit1$SDF)
#Peta Sebaran Penduga Parameter GWR di Setiap Lokasi
b.y<- gwr.fit1$SDF$Intercept
Jabar@data$y<- CoordK[,2]
spplot(Jabar,zcol="y",main="Peta Sebaran Y")

#Penduga Parameter beta1(X1)
b.x1<- gwr.fit1$SDF$x1
Jabar@data$b1<- b.x1
spplot(Jabar,zcol="b1",main="Peta Sebaran Penduga Parameter beta1 (x1)")

#beta2 (X2)
b.x2<- gwr.fit1$SDF$x2
Jabar@data$b2<- b.x2
spplot(Jabar,zcol="b2",main="Peta Sebaran Penduga Parameter beta2 (x2)")

CoordK <- coordinates(Jabar)
Jabar_sf <- st_as_sf(Jabar)
model_ols <- lm(UMK ~ PDRB + IPM, data = data_spasial)
model_ols <- predict(model_ols,data=data_spasial)
Jabar@data$rlbpred <- model_ols
Jabar@data$y <- data_spasial$UMK
Jabar@data$pred <- gwr.fit1$SDF$pred
Jabar_sf <- st_as_sf(Jabar)
Jabar_merged <- Jabar_sf %>% left_join(data_spasial, by = "id")
Jabar_merged <- Jabar_merged %>% filter(!is.na(pred) & !is.na(UMK))

# Buat plot dengan ggplot2 menggunakan geom_sf()
ggplot(Jabar_merged) + 
  geom_sf(aes(fill = UMK)) +  # Menggunakan geom_sf() untuk membuat peta
  scale_fill_gradient(low = "yellow", high = "red", name = "UMK") + 
  labs(title = "Peta Sebaran Aktual UMK") + 
  theme_minimal()

Jabar_merged <- Jabar_merged %>% filter(!is.na(pred) & !is.na(UMK))
# Buat plot dengan ggplot2
ggplot(Jabar_merged) + 
  geom_sf(aes(fill = pred)) +  # Gunakan 'pred' dari Jabar_merged
  scale_fill_gradient(low = "yellow", high = "red", name = "Prediksi UMK GWR") + 
  labs(title = "Peta Sebaran Prediksi UMK GWR") +
  theme_minimal()

# Buat plot dengan ggplot2
ggplot(Jabar_merged) + 
  geom_sf(aes(fill = rlbpred)) +  # Gunakan 'rlbpred' dari Jabar_merged
  scale_fill_gradient(low = "yellow", high = "red", name = "Prediksi UMK OLS") + 
  labs(title = "Peta Sebaran OLS UMK") +
  theme_minimal()

#Local significance test
t1 = gwr.fit1$SDF$x1 / gwr.fit1$SDF$x1_se
t2 = gwr.fit1$SDF$x2 / gwr.fit1$SDF$x2_se

pvalue1 = 2*pt(q=t1, df=26, lower.tail=FALSE)
pvalue2 = 2*pt(q=t2, df=26, lower.tail=FALSE)
#Variabel x1: Indeks Pengembangan Manusia (IPM)
data.frame(id, t1, pvalue1) %>% mutate(ket=ifelse(pvalue1<0.05,"Reject", "Not
Reject"))
##    id          t1    pvalue1         ket
## 1   1  2.58457090 0.01571664      Reject
## 2   2 -0.22626074 1.17723202 Not\nReject
## 3   3 -0.10743039 1.08472774 Not\nReject
## 4   4  0.64760888 0.52291575 Not\nReject
## 5   5  1.54540891 0.13433359 Not\nReject
## 6   6 -0.09447819 1.07454634 Not\nReject
## 7   7 -0.05878294 1.04642552 Not\nReject
## 8   8  0.08149176 0.93567497 Not\nReject
## 9   9  0.28059976 0.78123813 Not\nReject
## 10 10  1.48344473 0.14997831 Not\nReject
## 11 11  1.67702305 0.10552210 Not\nReject
## 12 12  1.14113516 0.26421569 Not\nReject
## 13 13 -0.95853350 1.65337693 Not\nReject
## 14 14 -0.32814124 1.25456710 Not\nReject
## 15 15  2.64524509 0.01366479      Reject
## 16 16  2.48219924 0.01983963      Reject
## 17 17 -0.69397760 1.50614935 Not\nReject
## 18 18 -0.11292098 1.08903943 Not\nReject
## 19 19  2.74133855 0.01092031      Reject
## 20 20  1.13951343 0.26487885 Not\nReject
## 21 21 -0.34005556 1.26345329 Not\nReject
## 22 22  0.27393027 0.78630276 Not\nReject
## 23 23  2.25538800 0.03275201      Reject
## 24 24  2.16489238 0.03975781      Reject
## 25 25 -0.77171592 1.55275952 Not\nReject
## 26 26 -0.02000397 1.01580707 Not\nReject
## 27 27 -0.10772874 1.08496210 Not\nReject
#Variabel x2: PDRD
data.frame(id, t2, pvalue2) %>% mutate(ket=ifelse(pvalue1<0.05,"Reject", "Not
Reject"))
##    id         t2      pvalue2         ket
## 1   1 3.63888082 0.0011900029      Reject
## 2   2 3.90407835 0.0006000937 Not\nReject
## 3   3 2.25278933 0.0329365092 Not\nReject
## 4   4 1.02891022 0.3129992015 Not\nReject
## 5   5 0.47706643 0.6373004504 Not\nReject
## 6   6 2.10258490 0.0453355146 Not\nReject
## 7   7 1.18141519 0.2481317559 Not\nReject
## 8   8 0.63100610 0.5335448233 Not\nReject
## 9   9 0.25813753 0.7983331722 Not\nReject
## 10 10 0.04878144 0.9614662791 Not\nReject
## 11 11 0.61596745 0.5432714212 Not\nReject
## 12 12 0.14403111 0.8865864909 Not\nReject
## 13 13 4.06545397 0.0003940713 Not\nReject
## 14 14 3.95329942 0.0005279938 Not\nReject
## 15 15 3.54599145 0.0015088935      Reject
## 16 16 2.56694407 0.0163644700      Reject
## 17 17 2.71244010 0.0116858489 Not\nReject
## 18 18 0.37827938 0.7082941345 Not\nReject
## 19 19 3.84922518 0.0006918859      Reject
## 20 20 3.75345065 0.0008863176 Not\nReject
## 21 21 1.94954532 0.0620934485 Not\nReject
## 22 22 0.25782215 0.7985739463 Not\nReject
## 23 23 2.34716481 0.0268080445      Reject
## 24 24 2.57211864 0.0161717637      Reject
## 25 25 2.44631961 0.0215070149 Not\nReject
## 26 26 1.99140991 0.0570393039 Not\nReject
## 27 27 0.82309710 0.4179386723 Not\nReject