This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

Memanggil Library

## 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')`
## Loading required package: sf
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
## Loading required package: sp
## corrplot 0.95 loaded
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
## 
##     Recode
## Loading required package: Matrix
## 
## 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
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()  masks Matrix::expand()
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ tidyr::pack()    masks Matrix::pack()
## ✖ dplyr::recode()  masks car::recode()
## ✖ dplyr::select()  masks raster::select()
## ✖ purrr::some()    masks car::some()
## ✖ tidyr::unpack()  masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
## 
## Loading required package: robustbase
## 
## Loading required package: Rcpp
## 
## Welcome to GWmodel version 2.4-2.

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

Set Working Directory

# Set working directory
setwd("C:/Users/suhaila/Downloads")

Load Data

# Load data
data <- read.csv("spasial  - Sheet1 (6).csv", dec=",")
jabar <- read_sf("jabar_modified.shp")

Mapping Variabel Dependen

#####MAPPING####
jabar$id<-c(1:27)
data$id<-c(1:27)
jabar_sf<-st_as_sf(jabar)
jabar_merged <- jabar_sf %>%
  left_join(data, by = "id")
# Menghitung centroid dari setiap kabupaten
jabar_centroid <- st_centroid(jabar_merged)
## Warning: st_centroid assumes attributes are constant over geometries
ggplot() +
  geom_sf(data = jabar_merged, aes(fill = Balita.Gizi.Buruk), color = "white") +
  theme_bw() +
  scale_fill_gradient(low = "pink", high = "brown") +
  # Menambahkan nama kabupaten
  geom_sf_text(data = jabar_centroid, aes(label = Kabupaten.kota), size = 3, color = "white") +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "right",
    axis.text.x = element_blank(),  # Remove x-axis labels
    axis.text.y = element_blank()   # Remove y-axis labels
  ) +
  labs(title = "", fill = "Gizi Buruk")

summary(data$Balita.Gizi.Buruk)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    55.0   125.0   242.7   380.5   846.0

Perbandingan Moran’s I Variabel Dependen

#QUEEN
# Menggunakan koordinat dari data untuk menghitung matriks bobot
coords <- cbind(data$longitude, data$latitude)

# 1. Matriks Bobot KNN (contoh 5 dan 6 tetangga terdekat)
k_neighbors1 <- 5
knn_weights1 <- knearneigh(coords, k = k_neighbors1) 
knn_nb1 <- knn2nb(knn_weights1)
knn_listw1 <- nb2listw(knn_nb1, style = "W")  # Menggunakan bobot standar

k_neighbors2 <- 6
knn_weights2 <- knearneigh(coords, k = k_neighbors2) 
knn_nb2 <- knn2nb(knn_weights2)
knn_listw2 <- nb2listw(knn_nb2, style = "W")  # Menggunakan bobot standar


# 2. Matriks Bobot Queen
queen_nb <- poly2nb(jabar_sf)  # Pastikan jabar_sf adalah data geospasial
queen_listw <- nb2listw(queen_nb, style = "W")  # Menggunakan bobot standar

# 3. Matriks Bobot Rook
rook_nb <- poly2nb(jabar_sf, queen = FALSE)  # Rook menggunakan queen = FALSE
rook_listw <- nb2listw(rook_nb, style = "W")  # Menggunakan bobot standar

# Load library
library(spdep)

# Misalkan 'variabel_target' adalah nama kolom dari data yang ingin diuji
variabel_target <- data$Balita.Gizi.Buruk  # Ganti 'nama_kolom_target' sesuai dengan nama kolom Anda

# Hitung nilai Moran's I untuk KNN
knn_moran1 <- moran.test(variabel_target, knn_listw1)
print(knn_moran1)
## 
##  Moran I test under randomisation
## 
## data:  variabel_target  
## weights: knn_listw1    
## 
## Moran I statistic standard deviate = -0.86937, p-value = 0.8077
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.123347931      -0.038461538       0.009533903
knn_moran2 <- moran.test(variabel_target, knn_listw2)
print(knn_moran2)
## 
##  Moran I test under randomisation
## 
## data:  variabel_target  
## weights: knn_listw2    
## 
## Moran I statistic standard deviate = -1.2181, p-value = 0.8884
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.143005613      -0.038461538       0.007366245
# Hitung nilai Moran's I untuk Queen
queen_moran <- moran.test(variabel_target, queen_listw)
print(queen_moran)
## 
##  Moran I test under randomisation
## 
## data:  variabel_target  
## weights: queen_listw    
## 
## Moran I statistic standard deviate = -2.1286, p-value = 0.9834
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.33447439       -0.03846154        0.01933863
# Hitung nilai Moran's I untuk Rook
rook_moran <- moran.test(variabel_target, rook_listw)
print(rook_moran)
## 
##  Moran I test under randomisation
## 
## data:  variabel_target  
## weights: rook_listw    
## 
## Moran I statistic standard deviate = -2.1286, p-value = 0.9834
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.33447439       -0.03846154        0.01933863
# Bandingkan hasil Moran's I
moran_results <- data.frame(
  Method = c("KNN 5","KNN 6", "Queen", "Rook"),
  Moran_I = c(knn_moran1$estimate[1],knn_moran2$estimate[1], queen_moran$estimate[1], rook_moran$estimate[1]),
  Expectation = c(knn_moran1$estimate[2],knn_moran2$estimate[2], queen_moran$estimate[2], rook_moran$estimate[2]),
  Variance = c(knn_moran1$estimate[3],knn_moran2$estimate[3], queen_moran$estimate[3], rook_moran$estimate[3]),
  P_Value = c(knn_moran1$p.value, knn_moran2$p.value, queen_moran$p.value, rook_moran$p.value)
)

print(moran_results)
##   Method    Moran_I Expectation    Variance   P_Value
## 1  KNN 5 -0.1233479 -0.03846154 0.009533903 0.8076765
## 2  KNN 6 -0.1430056 -0.03846154 0.007366245 0.8884035
## 3  Queen -0.3344744 -0.03846154 0.019338634 0.9833570
## 4   Rook -0.3344744 -0.03846154 0.019338634 0.9833570

Mapping Variabel Independen

# Daftar variabel untuk dipetakan
vars <- c("PDRB.per.kapita", "tingkat.pengangguran.terbuka", 
          "Indeks.keparahan.kemiskinan", "kepadatan.penduduk",
          "jumlah.dokter", "cakupan.imunisasi.dasar", 
          "keluarga.dengan.akses.fasilitas.sanitasi",
          "balita.berat.badan.kurang")

# Loop untuk membuat peta setiap variabel
for (var in vars) {
  # Plot untuk setiap variabel
  p <- ggplot() +
    geom_sf(data = jabar_merged, aes_string(fill = var), color = "white") +
    scale_fill_gradient(low = "lightblue", high = "darkblue", na.value = "grey90") +
    theme_bw() +
    geom_sf_text(data = jabar_centroid, aes(label = Kabupaten.kota), 
                 size = 3, color = "white", fontface = "bold") +
    theme(
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      legend.position = "right",
      axis.text.x = element_blank(),
      axis.text.y = element_blank()
    ) +
    labs(title = paste("Peta Distribusi", var), fill = var)
  
  # Menampilkan plot
  print(p)
}
## 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.

## Deskriptif Data

summary(data)
##  Kabupaten.kota     Balita.Gizi.Buruk PDRB.per.kapita
##  Length:27          Min.   :  0.0     Min.   :14205  
##  Class :character   1st Qu.: 55.0     1st Qu.:18882  
##  Mode  :character   Median :125.0     Median :23973  
##                     Mean   :242.7     Mean   :32041  
##                     3rd Qu.:380.5     3rd Qu.:35051  
##                     Max.   :846.0     Max.   :88554  
##  tingkat.pengangguran.terbuka Indeks.keparahan.kemiskinan kepadatan.penduduk
##  Min.   : 1.520               Min.   :0.0300              Min.   :  382.0   
##  1st Qu.: 6.535               1st Qu.:0.2000              1st Qu.:  824.5   
##  Median : 7.650               Median :0.2600              Median : 1449.0   
##  Mean   : 7.186               Mean   :0.2807              Mean   : 3873.3   
##  3rd Qu.: 8.500               3rd Qu.:0.3450              3rd Qu.: 5749.0   
##  Max.   :10.520               Max.   :0.6200              Max.   :15047.0   
##  jumlah.puskesmas jumlah.dokter   cakupan.imunisasi.dasar
##  Min.   : 10.00   Min.   : 45.0   Min.   : 82.67         
##  1st Qu.: 23.50   1st Qu.: 76.0   1st Qu.: 95.51         
##  Median : 38.00   Median : 96.0   Median : 97.14         
##  Mean   : 40.96   Mean   :114.1   Mean   : 97.52         
##  3rd Qu.: 50.50   3rd Qu.:135.5   3rd Qu.: 99.67         
##  Max.   :101.00   Max.   :237.0   Max.   :123.83         
##  keluarga.dengan.akses.fasilitas.sanitasi
##  Min.   :  54951                         
##  1st Qu.: 221323                         
##  Median : 468556                         
##  Mean   : 481529                         
##  3rd Qu.: 625727                         
##  Max.   :1119445                         
##  jumlah.penduduk.dengan.akses.berkelanjutan.air.minum.layak
##  Min.   : 182879                                           
##  1st Qu.: 917626                                           
##  Median :1164650                                           
##  Mean   :1422445                                           
##  3rd Qu.:2026029                                           
##  Max.   :2966480                                           
##  balita.berat.badan.kurang   longitude           latitude        
##  Min.   :  638             Min.   :-7701397   Min.   :  1081971  
##  1st Qu.: 2482             1st Qu.:-6994542   1st Qu.:106878230  
##  Median : 4482             Median :-6836531   Median :107512302  
##  Mean   : 5845             Mean   :-6605407   Mean   : 92565328  
##  3rd Qu.: 6771             3rd Qu.:-6515690   3rd Qu.:108016295  
##  Max.   :15799             Max.   : -687121   Max.   :108557818  
##        id      
##  Min.   : 1.0  
##  1st Qu.: 7.5  
##  Median :14.0  
##  Mean   :14.0  
##  3rd Qu.:20.5  
##  Max.   :27.0

Model OLS

#OLS#
y <- as.numeric(data$Balita.Gizi.Buruk)
x1 <- as.numeric(data$PDRB.per.kapita)
x2 <- data$tingkat.pengangguran.terbuka
x3 <- data$Indeks.keparahan.kemiskinan
x4 <- as.numeric(data$kepadatan.penduduk)
x5 <- data$jumlah.puskesmas
x6 <- data$jumlah.dokter
x7 <- data$cakupan.imunisasi.dasar
x8 <- data$keluarga.dengan.akses.fasilitas.sanitasi
x9 <- data$jumlah.penduduk.dengan.akses.berkelanjutan
x10 <- data$balita.berat.badan.kurang
model_ols <- lm(y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8  +x10, data = data)
summary(model_ols)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 + x10, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -211.60  -65.80  -32.12   71.38  215.83 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -9.022e+02  5.013e+02  -1.800  0.08871 . 
## x1           2.510e-03  1.633e-03   1.537  0.14173   
## x2          -3.486e+01  1.693e+01  -2.059  0.05424 . 
## x3           7.271e+02  2.447e+02   2.971  0.00819 **
## x4           1.170e-02  9.355e-03   1.251  0.22693   
## x6           1.229e+00  7.466e-01   1.646  0.11709   
## x7           6.988e+00  5.213e+00   1.340  0.19677   
## x8           2.958e-04  1.639e-04   1.804  0.08794 . 
## x10          1.735e-02  1.199e-02   1.446  0.16523   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 136 on 18 degrees of freedom
## Multiple R-squared:  0.7926, Adjusted R-squared:  0.7004 
## F-statistic: 8.597 on 8 and 18 DF,  p-value: 8.426e-05
model_ols1 <- lm(y ~ x3, data = data)
summary(model_ols1)
## 
## Call:
## lm(formula = y ~ x3, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -285.7 -173.2 -106.7  144.3  571.5 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    59.37     103.37   0.574   0.5709  
## x3            653.19     330.88   1.974   0.0595 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 235.6 on 25 degrees of freedom
## Multiple R-squared:  0.1349, Adjusted R-squared:  0.1003 
## F-statistic: 3.897 on 1 and 25 DF,  p-value: 0.05952

Uji Asumsi Heterogenitas OLS

##Heterogenitas##
bptest(model_ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 16.221, df = 8, p-value = 0.03932

Uji Asumsi Multikolinearitas OLS

##Multikolinearity##
vif(model_ols)
##       x1       x2       x3       x4       x6       x7       x8      x10 
## 1.693171 1.649316 1.643057 2.633562 2.241785 2.098027 3.527768 3.766920

Uji Asumsi Normalitas OLS

##Normalitas##
residual_ols <- residuals(model_ols)
shapiro.test(residual_ols)
## 
##  Shapiro-Wilk normality test
## 
## data:  residual_ols
## W = 0.95746, p-value = 0.3229

Cek Autokorelasi Residual Model OLS

##Autokorelasi##
coords <- cbind(data$longitude, data$latitude)
queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W") 

moran.test(residual_ols, queen_listw)
## 
##  Moran I test under randomisation
## 
## data:  residual_ols  
## weights: queen_listw    
## 
## Moran I statistic standard deviate = 0.51332, p-value = 0.3039
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.03355745       -0.03846154        0.01968442
# Hitung nilai Moran's I untuk KNN
nb_5 <- knn2nb(knearneigh(coords, k=5))
listw_5 <- nb2listw(nb_5, style="W")

nb_6 <- knn2nb(knearneigh(coords, k=6))
listw_6 <- nb2listw(nb_6, style="W")

queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W")

rook_nb <- poly2nb(jabar_sf, queen = FALSE)
rook_listw <- nb2listw(rook_nb, style = "W")

knn_moran_5 <- moran.test(residual_ols, listw_5)
print(knn_moran_5)
## 
##  Moran I test under randomisation
## 
## data:  residual_ols  
## weights: listw_5    
## 
## Moran I statistic standard deviate = 1.251, p-value = 0.1055
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.084773021      -0.038461538       0.009703921
knn_moran_6 <- moran.test(residual_ols, listw_6)
print(knn_moran_6)
## 
##  Moran I test under randomisation
## 
## data:  residual_ols  
## weights: listw_6    
## 
## Moran I statistic standard deviate = 1.8443, p-value = 0.03257
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.121219928      -0.038461538       0.007496251
knn_moran_7 <- moran.test(residual_ols, queen_listw)
print(knn_moran_7)
## 
##  Moran I test under randomisation
## 
## data:  residual_ols  
## weights: queen_listw    
## 
## Moran I statistic standard deviate = 0.51332, p-value = 0.3039
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.03355745       -0.03846154        0.01968442
knn_moran_8 <- moran.test(residual_ols, rook_listw)
print(knn_moran_8)
## 
##  Moran I test under randomisation
## 
## data:  residual_ols  
## weights: rook_listw    
## 
## Moran I statistic standard deviate = 0.51332, p-value = 0.3039
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.03355745       -0.03846154        0.01968442

Definisi Variabel pada Model GWR

######GWR#####
y <- as.numeric(data$Balita.Gizi.Buruk)
x1 <- as.numeric(data$PDRB.per.kapita)
x2 <- data$tingkat.pengangguran.terbuka
x3 <- data$Indeks.keparahan.kemiskinan
x4 <- as.numeric(data$kepadatan.penduduk)
x5 <- data$jumlah.puskesmas
x6 <- data$jumlah.dokter
x7 <- data$cakupan.imunisasi.dasar
x8 <- data$keluarga.dengan.akses.fasilitas.sanitasi
x9 <- data$jumlah.penduduk.dengan.akses.berkelanjutan
x10 <- data$balita.berat.badan.kurang

Penentuan Bandwtih Optimal

#optimal bandwith
bandwidth_optimal <- gwr.sel(
  y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8  +x10,
  data = data,
  coords = cbind(data$longitude, data$latitude),
  adapt = TRUE
)
## Adaptive q: 0.381966 CV score: 801523.8 
## Adaptive q: 0.618034 CV score: 736303.4 
## Adaptive q: 0.763932 CV score: 738355.4 
## Adaptive q: 0.6817313 CV score: 736847.3 
## Adaptive q: 0.586383 CV score: 740832.3 
## Adaptive q: 0.6180747 CV score: 736297.9 
## Adaptive q: 0.6423894 CV score: 734862.3 
## Adaptive q: 0.639664 CV score: 734853.9 
## Adaptive q: 0.6404948 CV score: 734855.9 
## Adaptive q: 0.6314176 CV score: 734867.8 
## Adaptive q: 0.6371811 CV score: 734851.6 
## Adaptive q: 0.6373114 CV score: 734851.6 
## Adaptive q: 0.6373521 CV score: 734851.6 
## Adaptive q: 0.6373928 CV score: 734851.6 
## Adaptive q: 0.6373521 CV score: 734851.6

Model GWR

#running GWR
model_gwr <- gwr(
  y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8  +x10,
  data = data,
  coords = cbind(data$longitude, data$latitude),
  adapt = bandwidth_optimal,
  hatmatrix = TRUE,
  se.fit = TRUE
)
print(model_gwr)
## Call:
## gwr(formula = y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 + x10, data = data, 
##     coords = cbind(data$longitude, data$latitude), adapt = bandwidth_optimal, 
##     hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Adaptive quantile: 0.6373521 (about 17 of 27 data points)
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -1.2595e+03 -1.2179e+03 -1.1656e+03 -1.0903e+03 -8.1918e+02
## x1            2.4099e-03  2.5537e-03  2.6297e-03  2.7340e-03  2.8809e-03
## x2           -3.9520e+01 -3.8227e+01 -3.5289e+01 -3.2945e+01 -3.2371e+01
## x3            6.5091e+02  8.4012e+02  8.7225e+02  1.0354e+03  1.0736e+03
## x4            1.0290e-02  1.3426e-02  1.6364e-02  2.0355e-02  2.1841e-02
## x6            1.0784e+00  1.1317e+00  1.3719e+00  1.5923e+00  1.6265e+00
## x7            6.6279e+00  7.3115e+00  8.2238e+00  9.0969e+00  9.6672e+00
## x8            3.0045e-04  3.0598e-04  3.1884e-04  3.4699e-04  3.6578e-04
## x10           1.4703e-02  1.6071e-02  1.7021e-02  1.7348e-02  1.8247e-02
##                 Global
## X.Intercept. -902.1530
## x1              0.0025
## x2            -34.8597
## x3            727.0945
## x4              0.0117
## x6              1.2290
## x7              6.9877
## x8              0.0003
## x10             0.0173
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 13.68578 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 13.31422 
## Sigma (residual: 2traceS - traceS'S): 134.4216 
## Effective number of parameters (model: traceS): 12.38756 
## Effective degrees of freedom (model: traceS): 14.61244 
## Sigma (model: traceS): 128.3115 
## Sigma (ML): 94.39412 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 379.5052 
## AIC (GWR p. 96, eq. 4.22): 334.5741 
## Residual sum of squares: 240576.8 
## Quasi-global R2: 0.8500062
#koefisien pada setiap lokasi
coef_gwr <- model_gwr$SDF@data[, c("x1",'x2','x3','x4','x6','x7','x8','x10')]

# Menampilkan ringkasan model
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
# Menampilkan ringkasan GWR Global (kuartil, median, dll.)
summary_gwr <- apply(coef_gwr, 2, summary)
print(summary_gwr)
##                  x1        x2        x3         x4       x6       x7
## Min.    0.002409909 -39.51971  650.9142 0.01029028 1.078432 6.627918
## 1st Qu. 0.002553685 -38.22654  840.1179 0.01342622 1.131732 7.311494
## Median  0.002629708 -35.28912  872.2458 0.01636355 1.371947 8.223832
## Mean    0.002628406 -35.81804  896.4846 0.01628267 1.346413 8.191521
## 3rd Qu. 0.002733953 -32.94530 1035.4412 0.02035543 1.592290 9.096944
## Max.    0.002880885 -32.37063 1073.6475 0.02184072 1.626486 9.667184
##                   x8        x10
## Min.    0.0003004484 0.01470271
## 1st Qu. 0.0003059780 0.01607122
## Median  0.0003188450 0.01702093
## Mean    0.0003257273 0.01680377
## 3rd Qu. 0.0003469948 0.01734795
## Max.    0.0003657821 0.01824698

Uji Asumsi Heterogenitas Model GWR

# Asumsi Heterogenitas Spasial
# Memasang paket lmtest
library(lmtest)

# Melakukan uji Breusch-Pagan pada model GWR
bp_test <- bptest(model_gwr$SDF$gwr.e ~ model_gwr$SDF$pred + 
                  I(model_gwr$SDF$pred^2))

# Menampilkan hasil uji
print(bp_test)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_gwr$SDF$gwr.e ~ model_gwr$SDF$pred + I(model_gwr$SDF$pred^2)
## BP = 3.9297, df = 2, p-value = 0.1402

Uji Asumsi Multikoleniaritas Model GWR

#Asumsi Multikoleniaritas
# Mengakses Condition Number dari model GWR
condition_number <- model_gwr$SDF@data$CN
print(condition_number)
## NULL
# Menyaring lokasi dengan potensi multikolinearitas tinggi
if (any(condition_number >= 30)) {
  print("Terdapat indikasi multikolinearitas di beberapa lokasi.")
} else {
  print("Tidak ada indikasi multikolinearitas.")
}
## [1] "Tidak ada indikasi multikolinearitas."

Uji Asumsi Normalitas Model GWR

##uji normality##
residual_gwr <- model_gwr$SDF$gwr.e
shapiro.test(residual_gwr)
## 
##  Shapiro-Wilk normality test
## 
## data:  residual_gwr
## W = 0.94385, p-value = 0.1516

Cek Autokorelasi Residual Model GWR

##Autokorelasi##
residual_gwr <- model_gwr$SDF$gwr.e
coords <- cbind(data$longitude, data$latitude)
queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W") 

moran.test(residual_gwr, queen_listw)
## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: queen_listw    
## 
## Moran I statistic standard deviate = 1.1675, p-value = 0.1215
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.12384698       -0.03846154        0.01932658
# Hitung nilai Moran's I untuk KNN
nb5 <- knn2nb(knearneigh(coords, k=5))
listw5 <- nb2listw(nb5, style="W")

nb6 <- knn2nb(knearneigh(coords, k=6))
listw6 <- nb2listw(nb6, style="W")

queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W")

rook_nb <- poly2nb(jabar_sf, queen = FALSE)
rook_listw <- nb2listw(rook_nb, style = "W")

knn_moran5 <- moran.test(residual_gwr, listw5)
print(knn_moran5)
## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: listw5    
## 
## Moran I statistic standard deviate = 1.6641, p-value = 0.04805
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.123973735      -0.038461538       0.009527975
knn_moran6 <- moran.test(residual_gwr, listw6)
print(knn_moran6)
## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: listw6    
## 
## Moran I statistic standard deviate = 2.2181, p-value = 0.01328
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.151848124      -0.038461538       0.007361712
knn_moran7 <- moran.test(residual_gwr, queen_listw)
print(knn_moran7)
## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: queen_listw    
## 
## Moran I statistic standard deviate = 1.1675, p-value = 0.1215
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.12384698       -0.03846154        0.01932658
knn_moran8 <- moran.test(residual_gwr, rook_listw)
print(knn_moran8)
## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: rook_listw    
## 
## Moran I statistic standard deviate = 1.1675, p-value = 0.1215
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.12384698       -0.03846154        0.01932658

Uji Kesesuaian Model GWR

# Menghitung residual untuk OLS dan GWR
residual_ols <- residuals(model_ols)
residual_gwr <- model_gwr$SDF$gwr.e

# Menghitung Sum of Squares (SS) untuk OLS dan GWR
SS_ols <- sum(residual_ols^2)  # Sum of Squares untuk model OLS
SS_gwr <- sum(residual_gwr^2)  # Sum of Squares untuk model GWR

# Degrees of Freedom (df) untuk OLS dan GWR
df_ols <- df.residual(model_ols)  # Derajat bebas model OLS
df_gwr <- sum(model_gwr$results$edf)  # Effective degrees of freedom untuk GWR

# Improvement dari GWR terhadap OLS
SS_improvement <- SS_ols - SS_gwr  # Perbaikan Sum of Squares
df_improvement <- df_ols - df_gwr  # Perbedaan derajat bebas

# Menghitung Mean Square (MS)
MS_ols <- SS_ols / df_ols  # Mean Square OLS
MS_improvement <- SS_improvement / df_improvement  # Mean Square improvement
MS_gwr <- SS_gwr / df_gwr  # Mean Square GWR residual

# Menghitung nilai F-statistik
F_value <- MS_improvement / MS_gwr  # F-statistik untuk membandingkan OLS dan GWR

# Menghitung p-value dari F-statistik
p_value <- pf(F_value, df_improvement, df_gwr, lower.tail = FALSE)

# Menampilkan tabel hasil perbandingan OLS dan GWR
tabel_hasil <- data.frame(
  "Keragaman Residual" = c("OLS Residual", "GWR Improvement", "GWR Residual"),
  "Df" = c(df_ols, df_improvement, df_gwr),
  "Sum Sq" = c(SS_ols, SS_improvement, SS_gwr),
  "Mean Sq" = c(MS_ols, MS_improvement, MS_gwr),
  "F" = c(NA, NA, F_value),  # Nilai F hanya muncul untuk GWR Improvement
  "p-value" = c(NA, NA, p_value)  # Menambahkan kolom p-value
)

# Menampilkan tabel
print(tabel_hasil)
##   Keragaman.Residual       Df    Sum.Sq  Mean.Sq        F   p.value
## 1       OLS Residual 18.00000 332689.22 18482.73       NA        NA
## 2    GWR Improvement  4.68578  92112.47 19657.87       NA        NA
## 3       GWR Residual 13.31422 240576.75 18069.16 1.087924 0.4084436

Pemilihan Model Terbaik

# Comparison of OLS and GWR models using AIC
aic_values <- data.frame(
  "MODEL" = c("GWR", "OLS"),
  "AIC" = c(model_gwr$results$AICh, AIC(model_ols))
)
aic_values
##   MODEL      AIC
## 1   GWR 334.5741
## 2   OLS 350.9391

Pemetaan Koefisien Model GWR

# Menampilkan koefisien dan p-value
coef_table <- as.data.frame(model_gwr$SDF@data)
print(coef_table[, c("x1", "x2", "x3", "x4", "x6", "x7", "x8", "x10")])
##             x1        x2        x3         x4       x6       x7           x8
## 1  0.002581543 -32.42267  842.1562 0.01345003 1.596316 8.985352 0.0003114371
## 2  0.002516878 -32.79590  842.4880 0.01400564 1.592530 9.256080 0.0003187514
## 3  0.002541646 -33.46475  834.4039 0.01412889 1.613452 9.483595 0.0003148354
## 4  0.002581912 -37.06315  916.7289 0.01767266 1.435713 9.128609 0.0003347314
## 5  0.002605647 -39.51971 1054.3479 0.02184072 1.130564 7.845026 0.0003617093
## 6  0.002620913 -38.79672 1073.6475 0.02179958 1.100690 7.398352 0.0003657821
## 7  0.002411897 -34.76915  651.4759 0.01031201 1.078434 6.627919 0.0003008658
## 8  0.002753446 -38.40670 1060.3660 0.02089069 1.132900 7.271602 0.0003546347
## 9  0.002791675 -38.26755 1049.1771 0.02043927 1.158007 7.323131 0.0003496963
## 10 0.002780088 -38.79380 1049.6114 0.02079337 1.151643 7.483515 0.0003501047
## 11 0.002761220 -39.50682 1024.1309 0.02070760 1.196582 8.034534 0.0003452199
## 12 0.002880885 -38.18553 1023.1510 0.01948506 1.216365 7.493316 0.0003375968
## 13 0.002845187 -38.59977  949.1596 0.01824999 1.371947 8.662774 0.0003225949
## 14 0.002690965 -35.28912  846.0544 0.01478916 1.604036 9.480896 0.0003059643
## 15 0.002409909 -34.75589  650.9142 0.01029028 1.080224 6.629720 0.0003004557
## 16 0.002637363 -33.09434  831.7497 0.01340241 1.626486 9.191968 0.0003037647
## 17 0.002632095 -37.19573  872.2458 0.01665451 1.533344 9.667184 0.0003197521
## 18 0.002636108 -37.83090 1069.2345 0.02111588 1.119273 7.188613 0.0003653089
## 19 0.002565723 -32.42859  844.1474 0.01357196 1.591087 9.001900 0.0003136732
## 20 0.002516569 -32.79627  842.7081 0.01401161 1.592050 9.253847 0.0003188450
## 21 0.002649549 -37.94518  928.5168 0.01815920 1.409497 9.065279 0.0003321409
## 22 0.002798780 -38.11182 1046.7515 0.02027158 1.164664 7.299856 0.0003487698
## 23 0.002629708 -32.62007  839.6238 0.01328203 1.606199 8.949836 0.0003059917
## 24 0.002590905 -32.37063  840.6120 0.01334146 1.599871 8.966540 0.0003099220
## 25 0.002714460 -32.53098  919.2182 0.01636355 1.492462 8.223832 0.0003007749
## 26 0.002411896 -34.76915  651.4753 0.01031201 1.078432 6.627918 0.0003008659
## 27 0.002409990 -34.75622  650.9873 0.01029099 1.080393 6.629882 0.0003004484
##           x10
## 1  0.01698135
## 2  0.01599085
## 3  0.01572963
## 4  0.01485624
## 5  0.01604837
## 6  0.01683812
## 7  0.01718762
## 8  0.01779161
## 9  0.01794844
## 10 0.01758856
## 11 0.01664293
## 12 0.01824698
## 13 0.01659912
## 14 0.01609408
## 15 0.01718457
## 16 0.01702093
## 17 0.01470271
## 18 0.01732378
## 19 0.01680949
## 20 0.01598902
## 21 0.01511680
## 22 0.01804406
## 23 0.01737211
## 24 0.01711993
## 25 0.01810237
## 26 0.01718762
## 27 0.01718450
##membuat plot-spplot##
coef_x1 <- model_gwr$SDF$x1
coef_x2 <- model_gwr$SDF$x2
coef_x3 <- model_gwr$SDF$x3
coef_x4 <- model_gwr$SDF$x4
coef_x5 <- model_gwr$SDF$x5
coef_x6 <- model_gwr$SDF$x6
coef_x7 <- model_gwr$SDF$x7
coef_x8 <- model_gwr$SDF$x8
coef_x9 <- model_gwr$SDF$x9
coef_x10 <- model_gwr$SDF$x10
##Membuat plot -gglot2##
gwr_df <- as.data.frame(model_gwr$SDF)
gwr_df
##       sum.w X.Intercept.          x1        x2        x3         x4       x6
## 1  17.37470   -1213.6388 0.002581543 -32.42267  842.1562 0.01345003 1.596316
## 2  17.24262   -1235.3255 0.002516878 -32.79590  842.4880 0.01400564 1.592530
## 3  17.28654   -1251.5650 0.002541646 -33.46475  834.4039 0.01412889 1.613452
## 4  16.57260   -1220.6626 0.002581912 -37.06315  916.7289 0.01767266 1.435713
## 5  16.04412   -1127.8756 0.002605647 -39.51971 1054.3479 0.02184072 1.130564
## 6  16.51345   -1098.2326 0.002620913 -38.79672 1073.6475 0.02179958 1.100690
## 7  17.95773    -819.2625 0.002411897 -34.76915  651.4759 0.01031201 1.078434
## 8  17.40602   -1088.2575 0.002753446 -38.40670 1060.3660 0.02089069 1.132900
## 9  17.51910   -1091.2367 0.002791675 -38.26755 1049.1771 0.02043927 1.158007
## 10 17.54177   -1102.5894 0.002780088 -38.79380 1049.6114 0.02079337 1.151643
## 11 16.89933   -1140.4684 0.002761220 -39.50682 1024.1309 0.02070760 1.196582
## 12 16.68149   -1101.4723 0.002880885 -38.18553 1023.1510 0.01948506 1.216365
## 13 16.07198   -1185.2137 0.002845187 -38.59977  949.1596 0.01824999 1.371947
## 14 16.46443   -1249.2509 0.002690965 -35.28912  846.0544 0.01478916 1.604036
## 15 17.95278    -819.1775 0.002409909 -34.75589  650.9142 0.01029028 1.080224
## 16 16.90314   -1228.8256 0.002637363 -33.09434  831.7497 0.01340241 1.626486
## 17 16.10434   -1259.4827 0.002632095 -37.19573  872.2458 0.01665451 1.533344
## 18 16.25037   -1084.1768 0.002636108 -37.83090 1069.2345 0.02111588 1.119273
## 19 17.45931   -1215.2295 0.002565723 -32.42859  844.1474 0.01357196 1.591087
## 20 17.24052   -1235.1637 0.002516569 -32.79627  842.7081 0.01401161 1.592050
## 21 16.52328   -1214.6486 0.002649549 -37.94518  928.5168 0.01815920 1.409497
## 22 17.44059   -1089.4467 0.002798780 -38.11182 1046.7515 0.02027158 1.164664
## 23 16.90283   -1210.1162 0.002629708 -32.62007  839.6238 0.01328203 1.606199
## 24 17.20623   -1212.0194 0.002590905 -32.37063  840.6120 0.01334146 1.599871
## 25 14.61569   -1165.6284 0.002714460 -32.53098  919.2182 0.01636355 1.492462
## 26 17.95767    -819.2620 0.002411896 -34.76915  651.4753 0.01031201 1.078432
## 27 17.95824    -819.2357 0.002409990 -34.75622  650.9873 0.01029099 1.080393
##          x7           x8        x10 X.Intercept._se       x1_se    x2_se
## 1  8.985352 0.0003114371 0.01698135        517.6940 0.001766681 17.64267
## 2  9.256080 0.0003187514 0.01599085        517.8146 0.001768092 17.59790
## 3  9.483595 0.0003148354 0.01572963        518.2396 0.001771664 17.77064
## 4  9.128609 0.0003347314 0.01485624        514.4076 0.001770659 17.60831
## 5  7.845026 0.0003617093 0.01604837        516.9336 0.001771186 17.42070
## 6  7.398352 0.0003657821 0.01683812        516.4236 0.001766353 17.36570
## 7  6.627919 0.0003008658 0.01718762        476.8166 0.001565843 16.12392
## 8  7.271602 0.0003546347 0.01779161        514.3254 0.001761514 17.36620
## 9  7.323131 0.0003496963 0.01794844        513.8364 0.001759948 17.37226
## 10 7.483515 0.0003501047 0.01758856        514.1727 0.001761698 17.40229
## 11 8.034534 0.0003452199 0.01664293        514.3090 0.001766845 17.53687
## 12 7.493316 0.0003375968 0.01824698        513.6876 0.001758481 17.47963
## 13 8.662774 0.0003225949 0.01659912        514.1146 0.001770960 18.15236
## 14 9.480896 0.0003059643 0.01609408        517.1864 0.001777348 18.43094
## 15 6.629720 0.0003004557 0.01718457        476.7919 0.001566241 16.12237
## 16 9.191968 0.0003037647 0.01702093        519.0886 0.001771919 18.05228
## 17 9.667184 0.0003197521 0.01470271        515.4180 0.001779836 18.11496
## 18 7.188613 0.0003653089 0.01732378        515.3387 0.001762553 17.33894
## 19 9.001900 0.0003136732 0.01680949        517.4553 0.001766254 17.59931
## 20 9.253847 0.0003188450 0.01598902        517.7968 0.001768050 17.59594
## 21 9.065279 0.0003321409 0.01511680        513.9146 0.001772232 17.73193
## 22 7.299856 0.0003487698 0.01804406        513.6803 0.001759406 17.36671
## 23 8.949836 0.0003059917 0.01737211        518.3167 0.001768187 17.81738
## 24 8.966540 0.0003099220 0.01711993        517.9354 0.001766975 17.67191
## 25 8.223832 0.0003007749 0.01810237        502.1377 0.001752437 16.92633
## 26 6.627918 0.0003008659 0.01718762        476.8166 0.001565843 16.12392
## 27 6.629882 0.0003004484 0.01718450        476.7852 0.001566194 16.12203
##       x3_se       x4_se     x6_se    x7_se        x8_se     x10_se       gwr.e
## 1  280.1248 0.011785217 1.0900669 5.246787 0.0002234409 0.01150636  175.563911
## 2  281.7249 0.011796945 1.0905776 5.244681 0.0002233612 0.01150831  224.985089
## 3  282.4750 0.011815502 1.0924069 5.256962 0.0002238421 0.01153569   10.531596
## 4  280.4130 0.011832863 1.0885846 5.202321 0.0002246551 0.01145189 -157.928870
## 5  280.1940 0.011990876 1.1097631 5.221019 0.0002278874 0.01148939 -136.410040
## 6  279.2165 0.011928970 1.1058535 5.216955 0.0002260751 0.01146668  183.868125
## 7  234.2089 0.008919521 0.7172558 4.934976 0.0001560973 0.01131969   51.455767
## 8  277.7656 0.011838650 1.0992553 5.194804 0.0002238297 0.01142352  -19.166125
## 9  277.7068 0.011809685 1.0961395 5.190393 0.0002232762 0.01140951   13.920179
## 10 277.8274 0.011841943 1.1007064 5.195290 0.0002239785 0.01142092  -18.162464
## 11 278.0422 0.011876261 1.1041305 5.202932 0.0002253272 0.01144142  -99.103542
## 12 278.8616 0.011774739 1.0931372 5.197286 0.0002228017 0.01140405   32.977222
## 13 279.7379 0.011798742 1.0989727 5.236319 0.0002249652 0.01148187  -77.443513
## 14 281.9936 0.011820438 1.0967432 5.276254 0.0002252315 0.01157563  -47.122907
## 15 234.2509 0.008920980 0.7169945 4.934879 0.0001560783 0.01131972   -3.426360
## 16 281.2947 0.011813487 1.0943764 5.279754 0.0002244971 0.01156964  -24.468625
## 17 282.5514 0.011864470 1.0958203 5.240003 0.0002257659 0.01152628  -41.515225
## 18 278.1688 0.011846946 1.0968161 5.206380 0.0002242518 0.01143639  -57.438057
## 19 280.1769 0.011781869 1.0896008 5.242031 0.0002232933 0.01149838  -30.103397
## 20 281.7158 0.011796694 1.0905404 5.244387 0.0002233540 0.01150779  102.358368
## 21 279.5843 0.011834147 1.0915713 5.204729 0.0002251158 0.01145524   -9.654120
## 22 277.6673 0.011798266 1.0945212 5.188585 0.0002230461 0.01140520   60.359725
## 23 280.2147 0.011796414 1.0916385 5.261465 0.0002239315 0.01153236  -56.977484
## 24 280.1342 0.011788656 1.0905096 5.250839 0.0002235630 0.01151303   -1.071166
## 25 268.5842 0.011665069 1.0540869 5.048940 0.0002135093 0.01133336   30.839991
## 26 234.2090 0.008919522 0.7172560 4.934976 0.0001560973 0.01131969  -86.364383
## 27 234.2456 0.008920811 0.7169707 4.934855 0.0001560762 0.01131972  163.070077
##           pred   pred.se   localR2 X.Intercept._se_EDF   x1_se_EDF x2_se_EDF
## 1   592.436089  86.56403 0.8470221            542.3462 0.001850808  18.48280
## 2   382.014911  57.04700 0.8392848            542.4725 0.001852288  18.43590
## 3   103.468404  53.38254 0.8384805            542.9178 0.001856030  18.61686
## 4   641.928870  75.16825 0.8362265            538.9033 0.001854976  18.44680
## 5   281.410040  81.03584 0.8488398            541.5496 0.001855529  18.25026
## 6   346.131875  72.46992 0.8575739            541.0153 0.001850466  18.19264
## 7   122.544233  63.01958 0.8390852            499.5222 0.001640407  16.89173
## 8    67.166125  70.57817 0.8698323            538.8172 0.001845396  18.19316
## 9   832.079821 121.01059 0.8718479            538.3049 0.001843755  18.19951
## 10  124.162464  59.67010 0.8687081            538.6571 0.001845588  18.23097
## 11  157.103542  44.93975 0.8591907            538.8000 0.001850980  18.37196
## 12  585.022778  88.90498 0.8758457            538.1490 0.001842219  18.31200
## 13  363.443513  88.69399 0.8596381            538.5963 0.001855291  19.01676
## 14  112.122907  67.87514 0.8467898            541.8144 0.001861984  19.30861
## 15  241.426360  83.02261 0.8390021            499.4963 0.001640825  16.89011
## 16  349.468625 112.19045 0.8490439            543.8071 0.001856296  18.91191
## 17  166.515225  43.93266 0.8361481            539.9618 0.001864591  18.97758
## 18   67.438057 105.07852 0.8623499            539.8788 0.001846484  18.16460
## 19   30.103397  78.70999 0.8456103            542.0961 0.001850362  18.43738
## 20  -57.358368  65.91062 0.8392705            542.4538 0.001852243  18.43385
## 21  445.654120 105.50254 0.8407535            538.3869 0.001856625  18.57631
## 22   -8.359725  86.77570 0.8726698            538.1414 0.001843188  18.19370
## 23  312.977484  94.61529 0.8507495            542.9985 0.001852386  18.66583
## 24   59.071166 106.11628 0.8480191            542.5991 0.001851117  18.51344
## 25  -12.839991  81.32986 0.8723448            526.0492 0.001835886  17.73235
## 26  204.364383  94.33970 0.8390851            499.5223 0.001640407  16.89173
## 27 -139.070077  67.63229 0.8390076            499.4893 0.001640775  16.88975
##    x3_se_EDF   x4_se_EDF x6_se_EDF x7_se_EDF    x8_se_EDF x10_se_EDF pred.se.1
## 1   293.4641 0.012346420 1.1419750  5.496635 0.0002340809 0.01205428  90.68614
## 2   295.1404 0.012358706 1.1425100  5.494428 0.0002339974 0.01205633  59.76353
## 3   295.9262 0.012378147 1.1444265  5.507294 0.0002345013 0.01208501  55.92458
## 4   293.7661 0.012396335 1.1404221  5.450051 0.0002353530 0.01199722  78.74770
## 5   293.5366 0.012561872 1.1626091  5.469640 0.0002387392 0.01203650  84.89471
## 6   292.5125 0.012497018 1.1585134  5.465382 0.0002368406 0.01201271  75.92088
## 7   245.3618 0.009344261 0.7514109  5.169975 0.0001635305 0.01185872  66.02052
## 8   290.9926 0.012402397 1.1516009  5.442176 0.0002344883 0.01196750  73.93905
## 9   290.9309 0.012372053 1.1483368  5.437555 0.0002339085 0.01195283 126.77302
## 10  291.0573 0.012405846 1.1531212  5.442686 0.0002346442 0.01196478  62.51155
## 11  291.2823 0.012441799 1.1567083  5.450691 0.0002360571 0.01198625  47.07975
## 12  292.1407 0.012335443 1.1451915  5.444777 0.0002334114 0.01194710  93.13857
## 13  293.0587 0.012360589 1.1513049  5.485668 0.0002356779 0.01202863  92.91753
## 14  295.4219 0.012383318 1.1489692  5.527505 0.0002359568 0.01212685  71.10730
## 15  245.4057 0.009345790 0.7511372  5.169874 0.0001635106 0.01185876  86.97608
## 16  294.6897 0.012376036 1.1464897  5.531171 0.0002351874 0.01212057 117.53287
## 17  296.0063 0.012429446 1.1480024  5.489527 0.0002365166 0.01207515  46.02470
## 18  291.4150 0.012411088 1.1490456  5.454304 0.0002349305 0.01198098 110.08227
## 19  293.5187 0.012342913 1.1414867  5.491652 0.0002339264 0.01204592  82.45809
## 20  295.1309 0.012358443 1.1424710  5.494121 0.0002339899 0.01205578  69.04923
## 21  292.8979 0.012397680 1.1435510  5.452574 0.0002358356 0.01200073 110.52648
## 22  290.8896 0.012360091 1.1466415  5.435662 0.0002336673 0.01194830  90.90789
## 23  293.5583 0.012358150 1.1436215  5.512012 0.0002345949 0.01208152  99.12079
## 24  293.4740 0.012350023 1.1424388  5.500880 0.0002342089 0.01206127 111.16946
## 25  281.3740 0.012220550 1.1042817  5.289366 0.0002236765 0.01187305  85.20272
## 26  245.3618 0.009344262 0.7514111  5.169976 0.0001635305 0.01185872  98.83208
## 27  245.4002 0.009345613 0.7511123  5.169849 0.0001635084 0.01185875  70.85289
##     coord.x   coord.y
## 1  -6479679 106824965
## 2  -6915727 106932576
## 3  -6822558 107139542
## 4  -7012851 107528627
## 5  -7202988 107885592
## 6  -7361212 108112488
## 7  -7325788   1083514
## 8  -6976233 108482982
## 9  -6764507 108478858
## 10 -6836531 108227876
## 11 -6861126 107920102
## 12 -6327234 108322903
## 13 -6571549 107762495
## 14 -6551701 107446541
## 15 -6301721  10730529
## 16 -6364614 107172509
## 17 -6840651 107512302
## 18 -7701397 108495155
## 19 -6594946 106794913
## 20 -6918366 106931496
## 21 -6910786 107609757
## 22 -6707076 108557818
## 23 -6236221 106994293
## 24 -6394473 106822692
## 25  -687121 107555486
## 26 -7316436   1081971
## 27 -7362487  10855887
jabar_merged1 <- jabar_merged %>%
  bind_cols(gwr_df %>%
              select(X.Intercept., x1, x2, x3, x4, x6, x7, x8, x10))
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = X.Intercept.), color = NA) +
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk intercept", fill = "Koefisien Beta0") +
  theme_minimal()

# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = x1), color = NA) +
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk x1", fill = "Koefisien Beta1") +
  theme_minimal()

# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = x2), color = NA) +
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk x2", fill = "Koefisien Beta2") +
  theme_minimal()

# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = x3), color = NA) +
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk x3", fill = "Koefisien Beta3") +
  theme_minimal()

# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = x4), color = NA) +
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk x4", fill = "Koefisien Beta4") +
  theme_minimal()

# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = x6), color = NA) +
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk x5", fill = "Koefisien Beta5") +
  theme_minimal()

# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = x7), color = NA) +
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk x6", fill = "Koefisien Beta6") +
  theme_minimal()

# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = x8), color = NA) + # Gunakan 'fill' untuk mengisi warna berdasarkan x1
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk x7", fill = "Koefisien Beta7") +
  theme_minimal()

# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
  geom_sf(aes(fill = x10), color = NA) +
  scale_fill_viridis_c() +
  geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
  labs(title = "Koefisien GWR untuk x8", fill = "Koefisien Beta8") +
  theme_minimal()

## Penentuan Kelompok Variabel Signifikan Tiap Kabupaten/Kota

# Ambil hasil koefisien dan standar error
hasil_gwr <- as.data.frame(model_gwr$SDF)

# Contoh kolom koefisien dan p-value yang dihasilkan
# Kolom biasanya bernama seperti "x1", "x1_se", "x2", "x2_se", dll.
# Periksa struktur hasil GWR
str(hasil_gwr)
## 'data.frame':    27 obs. of  35 variables:
##  $ sum.w              : num  17.4 17.2 17.3 16.6 16 ...
##  $ X.Intercept.       : num  -1214 -1235 -1252 -1221 -1128 ...
##  $ x1                 : num  0.00258 0.00252 0.00254 0.00258 0.00261 ...
##  $ x2                 : num  -32.4 -32.8 -33.5 -37.1 -39.5 ...
##  $ x3                 : num  842 842 834 917 1054 ...
##  $ x4                 : num  0.0135 0.014 0.0141 0.0177 0.0218 ...
##  $ x6                 : num  1.6 1.59 1.61 1.44 1.13 ...
##  $ x7                 : num  8.99 9.26 9.48 9.13 7.85 ...
##  $ x8                 : num  0.000311 0.000319 0.000315 0.000335 0.000362 ...
##  $ x10                : num  0.017 0.016 0.0157 0.0149 0.016 ...
##  $ X.Intercept._se    : num  518 518 518 514 517 ...
##  $ x1_se              : num  0.00177 0.00177 0.00177 0.00177 0.00177 ...
##  $ x2_se              : num  17.6 17.6 17.8 17.6 17.4 ...
##  $ x3_se              : num  280 282 282 280 280 ...
##  $ x4_se              : num  0.0118 0.0118 0.0118 0.0118 0.012 ...
##  $ x6_se              : num  1.09 1.09 1.09 1.09 1.11 ...
##  $ x7_se              : num  5.25 5.24 5.26 5.2 5.22 ...
##  $ x8_se              : num  0.000223 0.000223 0.000224 0.000225 0.000228 ...
##  $ x10_se             : num  0.0115 0.0115 0.0115 0.0115 0.0115 ...
##  $ gwr.e              : num  175.6 225 10.5 -157.9 -136.4 ...
##  $ pred               : num  592 382 103 642 281 ...
##  $ pred.se            : num  86.6 57 53.4 75.2 81 ...
##  $ localR2            : num  0.847 0.839 0.838 0.836 0.849 ...
##  $ X.Intercept._se_EDF: num  542 542 543 539 542 ...
##  $ x1_se_EDF          : num  0.00185 0.00185 0.00186 0.00185 0.00186 ...
##  $ x2_se_EDF          : num  18.5 18.4 18.6 18.4 18.3 ...
##  $ x3_se_EDF          : num  293 295 296 294 294 ...
##  $ x4_se_EDF          : num  0.0123 0.0124 0.0124 0.0124 0.0126 ...
##  $ x6_se_EDF          : num  1.14 1.14 1.14 1.14 1.16 ...
##  $ x7_se_EDF          : num  5.5 5.49 5.51 5.45 5.47 ...
##  $ x8_se_EDF          : num  0.000234 0.000234 0.000235 0.000235 0.000239 ...
##  $ x10_se_EDF         : num  0.0121 0.0121 0.0121 0.012 0.012 ...
##  $ pred.se.1          : num  90.7 59.8 55.9 78.7 84.9 ...
##  $ coord.x            : num  -6479679 -6915727 -6822558 -7012851 -7202988 ...
##  $ coord.y            : num  1.07e+08 1.07e+08 1.07e+08 1.08e+08 1.08e+08 ...
library(dplyr)
# Tambahkan kolom Kabupaten.kota ke hasil_gwr sebelum melakukan analisis
hasil_gwr <- hasil_gwr %>%
  mutate(Kabupaten.kota = data$Kabupaten.kota)  # Pastikan 'data' adalah nama dataframe awal yang berisi Kabupaten.kota
# Hitung p-value manual (z = koefisien / se)
hasil_gwr <- hasil_gwr %>%
  mutate(
    p_value_x1 = 2 * (1 - pnorm(abs(coef_x1 / x1_se))),
    p_value_x2 = 2 * (1 - pnorm(abs(coef_x2 / x2_se))),
    p_value_x3 = 2 * (1 - pnorm(abs(coef_x3 / x3_se))),
    p_value_x4 = 2 * (1 - pnorm(abs(coef_x4 / x4_se))),
    p_value_x6 = 2 * (1 - pnorm(abs(coef_x6 / x6_se))),
    p_value_x7 = 2 * (1 - pnorm(abs(coef_x7 / x7_se))),
    p_value_x8 = 2 * (1 - pnorm(abs(coef_x8 / x8_se))),
    p_value_x10 = 2 * (1 - pnorm(abs(coef_x10 / x10_se)))
  ) %>%
  # Tentukan variabel mana yang signifikan (p < 0.05)
  mutate(
    signifikan_x1 = ifelse(p_value_x1 < 0.05, "x1",NA),
    signifikan_x2 = ifelse(p_value_x2 < 0.05, "x2",NA),
    signifikan_x3 = ifelse(p_value_x3 < 0.05, "x3",NA),
    signifikan_x4 = ifelse(p_value_x4 < 0.05, "x4",NA),
    signifikan_x6 = ifelse(p_value_x6 < 0.05, "x6",NA),
    signifikan_x7 = ifelse(p_value_x7 < 0.05, "x7",NA),
    signifikan_x8 = ifelse(p_value_x8 < 0.05, "x8",NA),
    signifikan_x10 = ifelse(p_value_x10 < 0.05, "x10",NA)
  )
# Tambahkan kolom Kabupaten.kota ke hasil_gwr sebelum melakukan analisis
hasil_gwr <- hasil_gwr %>%
  mutate(Kabupaten.kota = data$Kabupaten.kota)  # Pastikan 'data' adalah nama dataframe awal yang berisi Kabupaten.kota
# Menggabungkan variabel signifikan menjadi satu kolom
hasil_gwr <- hasil_gwr %>%
  mutate(
    variabel_signifikan = paste(
      signifikan_x1, signifikan_x2, signifikan_x3, signifikan_x4,
      signifikan_x6, signifikan_x7, signifikan_x8, signifikan_x10, sep = ", "
    )
  ) %>%
  # Hapus NA dalam hasil variabel_signifikan
  mutate(variabel_signifikan = gsub("NA, |, NA", "", variabel_signifikan)) %>%
  mutate(variabel_signifikan = ifelse(variabel_signifikan == "NA", "Tidak ada", variabel_signifikan))

# Membuat tabel akhir
tabel_hasil <- hasil_gwr %>%
  dplyr::select(Kabupaten.kota, variabel_signifikan)
# Tampilkan tabel hasil
print(tabel_hasil)
##      Kabupaten.kota variabel_signifikan
## 1             bogor                  x3
## 2          sukabumi                  x3
## 3           cianjur                  x3
## 4           bandung              x2, x3
## 5             garut              x2, x3
## 6       tasikmalaya              x2, x3
## 7            ciamis              x2, x3
## 8          kuningan              x2, x3
## 9           cirebon              x2, x3
## 10       majalengka              x2, x3
## 11         sumedang              x2, x3
## 12        indramayu              x2, x3
## 13           subang              x2, x3
## 14       purwakarta                  x3
## 15         karawang              x2, x3
## 16           bekasi                  x3
## 17    bandung barat              x2, x3
## 18      pangandaran              x2, x3
## 19       kota bogor                  x3
## 20    kota sukabumi                  x3
## 21     kota bandung              x2, x3
## 22     kota cirebon              x2, x3
## 23      kota bekasi                  x3
## 24       kota depok                  x3
## 25      kota cimahi                  x3
## 26 kota tasikmalaya              x2, x3
## 27      kota banjar              x2, x3
view(tabel_hasil)