Nama: Fatih Zahrani
NPM: 140610230014
Dosen Pengampu: I Gede Nyoman Mindra Jaya, M.Si., Ph.D.

Mata Kuliah: Statistik Spasial

S1 Statistika FMIPA UNPAD

Dataset

https://bit.ly/4stLbr9

Syntax Analisis

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(sf)
library(ggplot2)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(openxlsx)
## Warning: package 'openxlsx' was built under R version 4.3.3
# -------------------------- MASUKKAN PLOT SPASIAL -----------------------------
Indo <- st_read("D:/Kuliah/Sem 5/Spasial/RBI_50K_2023_Jawa Barat.x26272/RBI_50K_2023_Jawa Barat.shp")
## Reading layer `RBI_50K_2023_Jawa Barat' from data source 
##   `D:\Kuliah\Sem 5\Spasial\RBI_50K_2023_Jawa Barat.x26272\RBI_50K_2023_Jawa Barat.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 25 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY, XYZ
## Bounding box:  xmin: 106.3703 ymin: -7.82099 xmax: 108.8468 ymax: -5.806538
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
names(Indo)
##  [1] "NAMOBJ"     "FCODE"      "REMARK"     "METADATA"   "SRS_ID"    
##  [6] "KDBBPS"     "KDCBPS"     "KDCPUM"     "KDEBPS"     "KDEPUM"    
## [11] "KDPBPS"     "KDPKAB"     "KDPPUM"     "LUASWH"     "TIPADM"    
## [16] "WADMKC"     "WADMKD"     "WADMKK"     "WADMPR"     "WIADKC"    
## [21] "WIADKK"     "WIADPR"     "WIADKD"     "SHAPE_Leng" "SHAPE_Area"
## [26] "geometry"
Indo$WADMKK
##  [1] "Bandung"          "Bandung Barat"    "Bekasi"           "Bogor"           
##  [5] "Ciamis"           "Cianjur"          "Cirebon"          "Garut"           
##  [9] "Indramayu"        "Karawang"         "Kota Bandung"     "Kota Banjar"     
## [13] "Kota Bekasi"      "Kota Bogor"       "Kota Cimahi"      "Kota Cirebon"    
## [17] "Kota Depok"       "Kota Sukabumi"    "Kota Tasikmalaya" "Kuningan"        
## [21] "Majalengka"       "Pangandaran"      "Purwakarta"       "Subang"          
## [25] "Sukabumi"         "Sumedang"         "Tasikmalaya"
Indo <- st_zm(Indo, drop = TRUE, what = "ZM")

ggplot(Indo) +
  geom_sf(fill = "lightblue", color = "black") +
  labs(title = "Peta Kabupaten/Kota di Provinsi Jawa Barat") +
  theme_minimal()

# ------------------------------- MASUKKAN DATA --------------------------------
data <- read.xlsx("D:/Kuliah/Sem 5/Spasial/UAS Spasial/Laporan/TBC fix.xlsx")
head(data)
##        kab_kota    JK      JP HIV   ST    DB Sanitasi Faskes Miskin   KP
## 1       Bandung 14012 3753120 598 8.87 44200    96.79     79  239.9 2156
## 2 Bandung Barat  4569 1884190 153 4.95 15911    94.91     38  179.7 1468
## 3        Bekasi 15458 3273870 893 1.70 38666    97.00    102  204.5 2617
## 4         Bogor 29110 5682300 831 1.91 65751    88.77    129  446.8 1899
## 5        Ciamis  3477 1259230 110 2.52 14906    95.16     42   90.8  789
## 6       Cianjur 12197 2584990 293 2.51 15838    69.99     54  239.3  712
##    Rate_1k
## 1 3.733427
## 2 2.424915
## 3 4.721629
## 4 5.122926
## 5 2.761211
## 6 4.718393
data$`kab_kota`
##  [1] "Bandung"          "Bandung Barat"    "Bekasi"           "Bogor"           
##  [5] "Ciamis"           "Cianjur"          "Cirebon"          "Garut"           
##  [9] "Indramayu"        "Karawang"         "Kota Bandung"     "Kota Banjar"     
## [13] "Kota Bekasi"      "Kota Bogor"       "Kota Cimahi"      "Kota Cirebon"    
## [17] "Kota Depok"       "Kota Sukabumi"    "Kota Tasikmalaya" "Kuningan"        
## [21] "Majalengka"       "Pangandaran"      "Purwakarta"       "Subang"          
## [25] "Sukabumi"         "Sumedang"         "Tasikmalaya"
setdiff(data$`kab_kota`, Indo$WADMKK)
## character(0)
# Merged data
jabar_merged <- Indo %>%
  left_join(data, by = c("WADMKK" = "kab_kota"))

# --------------------------- BUAT PLOT SPASIAL --------------------------------
# JK
ggplot(jabar_merged) +
  geom_sf(aes(fill = JK)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Jumlah Kasus TBC", fill = "Jumlah Kasus (orang)") +
  theme_minimal()

# JP
ggplot(jabar_merged) +
  geom_sf(aes(fill = JP)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Jumlah Penduduk di Jawa Barat", fill = "Jumlah Penduduk (Jiwa)") +
  theme_minimal()

# HIV
ggplot(jabar_merged) +
  geom_sf(aes(fill = HIV)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Jumlah Kasus HIV di Jawa Barat", fill = "Jumlah Kasus (orang)") +
  theme_minimal()

# ST
ggplot(jabar_merged) +
  geom_sf(aes(fill = ST)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Kasus Stunting di Jawa Barat", fill = "Persentase Kasus (%)") +
  theme_minimal()

# DB
ggplot(jabar_merged) +
  geom_sf(aes(fill = DB)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Jumlah Penderita Diabetes di Jawa Barat", fill = "Jumlah Penderita (orang)") +
  theme_minimal()

# Sanitasi
ggplot(jabar_merged) +
  geom_sf(aes(fill = Sanitasi)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Keluarga dengan Sanitasi Layak Jawa Barat", fill = "Persentase Keluarga (%)") +
  theme_minimal()

# Faskes
ggplot(jabar_merged) +
  geom_sf(aes(fill = Faskes)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Jumlah Fasilitas Kesehatan di Jawa Barat", fill = "Fasilitas Kesehatan (unit)") +
  theme_minimal()

# Miskin
ggplot(jabar_merged) +
  geom_sf(aes(fill = Miskin)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Penduduk Miskin di Jawa Barat", fill = "Jumlah Penduduk Miskin (ribu orang)") +
  theme_minimal()

# KP
ggplot(jabar_merged) +
  geom_sf(aes(fill = KP)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Kepadatan Penduduk di Jawa Barat", fill = "Kepadatan Penduduk (jiwa/km2)") +
  theme_minimal()

# Rate_1k
ggplot(jabar_merged) +
  geom_sf(aes(fill = Rate_1k)) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Sebaran Prevalensi Kasus TBC di Jawa Barat", fill = "Prevalensi Kasus (per 1000 penduduk)") +
  theme_minimal()

# ------------------------- PLOT JARINGAN TETANGGA -----------------------------
# Ubah ke format Spatial untuk fungsi spdep
jabar_sp <- as_Spatial(jabar_merged)
row.names(jabar_sp) <- jabar_sp$WADMKK
# Buat daftar ketetanggaan (queen contiguity)
W <- poly2nb(jabar_sp, row.names = row.names(jabar_sp), queen = TRUE)
WL <- nb2listw(W, style = "W", zero.policy = TRUE)
summary(W)
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 106 
## Percentage nonzero weights: 14.54047 
## Average number of links: 3.925926 
## Link number distribution:
## 
## 1 2 3 4 5 6 7 8 
## 4 3 6 4 1 7 1 1 
## 4 least connected regions:
## Kota Banjar Kota Bogor Kota Cirebon Kota Sukabumi with 1 link
## 1 most connected region:
## Bogor with 8 links
# Buat Plot
CoordJ <- coordinates(jabar_sp)
plot(jabar_sp, axes = TRUE, col = "gray90",
     main = "Jaringan Ketetanggaan Kabupaten/Kota di Jawa Barat")
text(CoordJ[,1], CoordJ[,2], labels = row.names(jabar_sp),
     col = "black", cex = 0.6, pos = 1.5)
points(CoordJ[,1], CoordJ[,2], pch = 19, cex = 0.7, col = "blue")
plot.nb(W, CoordJ, add = TRUE, col = "red", lwd = 1.5)

# ------------------------- CEK AUTOKORELASI SPASIAL ---------------------------

################################
# Variabel Y 
################################

# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$JK,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$JK  
## weights: WL    
## 
## Moran I statistic standard deviate = 3.2168, p-value = 0.0006482
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.38398688       -0.03846154        0.01724651
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$JK,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$JK 
## weights: WL   
## 
## Geary C statistic standard deviate = 0.15773, p-value = 0.4373
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.96891958        1.00000000        0.03882643
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$JK,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$JK, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$JK 
## weights: WL   
## 
## standard deviate = 3.7021, p-value = 0.0001069
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       5.552318e-02       3.846154e-02       2.123951e-05
# ----------- Peta Getis-Ord Local G-----------

# Hitung Getis-Ord Local G*
Local_Getis <- localG(jabar_merged$JK, WL, zero.policy = TRUE)
                      
jabar_merged$GiZ <- as.numeric(Local_Getis)

jabar_merged$Gi_cluster <- cut(jabar_merged$GiZ,
                               breaks = c(-Inf, -1.96, 1.96, Inf),
                               labels = c("Coldspot", "Tidak signifikan", "Hotspot")
)

# 3. Now run your ggplot
library(ggplot2)
ggplot(jabar_merged) +
  geom_sf(aes(fill = Gi_cluster), color = "white", linewidth = 0.2) +
  scale_fill_manual(
    values = c(
      "Hotspot" = "red",
      "Coldspot" = "blue",
      "Tidak signifikan" = "grey80"
    )
  ) +
  labs(
    title = "Peta Hotspot dan Coldspot Getis–Ord (Gi*)",
    subtitle = "Variabel: Jumlah Kasus TBC (JK) di Provinsi Jawa Barat",
    fill = "Kategori"
  ) +
  theme_minimal()

# -----------LISA---------------
# Local Moran’s I
Local_Moran <- localmoran(jabar_merged$JK, WL, zero.policy = TRUE)

# Tambahkan hasil LISA ke data spasial
colnames(Local_Moran) <- c("Ii", "E.Ii", "Var.Ii", "Z.Ii", "P.value")
jabar_merged$Ii <- Local_Moran[, "Ii"]
jabar_merged$Z.Ii <- Local_Moran[, "Z.Ii"]
jabar_merged$P.value <- Local_Moran[, "P.value"]

# Klasifikasi Klaster LISA
mean_var <- mean(jabar_merged$JK, na.rm = TRUE)
mean_lag <- lag.listw(WL, jabar_merged$JK)

jabar_merged$cluster <- NA
jabar_merged$cluster[jabar_merged$JK >= mean_var & mean_lag >= mean(mean_lag)] <- "High-High"
jabar_merged$cluster[jabar_merged$JK <= mean_var & mean_lag <= mean(mean_lag)] <- "Low-Low"
jabar_merged$cluster[jabar_merged$JK >= mean_var & mean_lag <= mean(mean_lag)] <- "High-Low"
jabar_merged$cluster[jabar_merged$JK <= mean_var & mean_lag >= mean(mean_lag)] <- "Low-High"
jabar_merged$cluster[jabar_merged$P.value > 0.05] <- "Non-signifikan"

# Visualisasi LISA
ggplot(jabar_merged) +
  geom_sf(aes(fill = cluster)) +
  scale_fill_manual(values = c(
    "High-High" = "red",
    "Low-Low" = "blue",
    "High-Low" = "orange",
    "Low-High" = "green",
    "Non-signifikan" = "grey80"
  )) +
  labs(
    title = "Peta Klaster Autokorelasi Lokal (LISA)",
    subtitle = "Variabel: Jumlah Kasus TBC di Jawa Barat",
    fill = "Tipe Klaster"
  ) +
  theme_minimal()

# ---------------------------------- MODEL OLS ---------------------------------
ols_model <- lm(JK ~ 
                  JP +
                  HIV +
                  ST +
                  DB +
                  Sanitasi +
                  Faskes +
                  Miskin +
                  KP,
                data = data
)
summary(ols_model)
## 
## Call:
## lm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + 
##     KP, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3332.1 -1479.8  -421.3  1143.3  4128.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  2.368e+03  4.129e+03   0.574   0.5734  
## JP           5.013e-04  6.266e-04   0.800   0.4341  
## HIV          2.195e+00  3.800e+00   0.578   0.5706  
## ST          -1.819e+02  1.332e+02  -1.365   0.1890  
## DB          -8.767e-02  5.001e-02  -1.753   0.0966 .
## Sanitasi    -4.894e+01  3.921e+01  -1.248   0.2279  
## Faskes       1.399e+02  5.463e+01   2.561   0.0196 *
## Miskin       2.198e+01  1.283e+01   1.714   0.1038  
## KP           2.768e-01  1.538e-01   1.800   0.0887 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2217 on 18 degrees of freedom
## Multiple R-squared:  0.9127, Adjusted R-squared:  0.8738 
## F-statistic: 23.51 on 8 and 18 DF,  p-value: 5.095e-08
# Uji Multikolinearitas
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
vif(ols_model)
##        JP       HIV        ST        DB  Sanitasi    Faskes    Miskin        KP 
##  3.132431  8.096913  1.536464  4.041859  1.804266 13.726707  8.199660  2.726582
ols_model <- lm(JK ~ 
                  JP +
                  HIV +
                  ST +
                  DB +
                  Sanitasi +
                  Faskes +
                  Miskin +
                  KP,
                data = data
)
summary(ols_model)
## 
## Call:
## lm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + 
##     KP, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3332.1 -1479.8  -421.3  1143.3  4128.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  2.368e+03  4.129e+03   0.574   0.5734  
## JP           5.013e-04  6.266e-04   0.800   0.4341  
## HIV          2.195e+00  3.800e+00   0.578   0.5706  
## ST          -1.819e+02  1.332e+02  -1.365   0.1890  
## DB          -8.767e-02  5.001e-02  -1.753   0.0966 .
## Sanitasi    -4.894e+01  3.921e+01  -1.248   0.2279  
## Faskes       1.399e+02  5.463e+01   2.561   0.0196 *
## Miskin       2.198e+01  1.283e+01   1.714   0.1038  
## KP           2.768e-01  1.538e-01   1.800   0.0887 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2217 on 18 degrees of freedom
## Multiple R-squared:  0.9127, Adjusted R-squared:  0.8738 
## F-statistic: 23.51 on 8 and 18 DF,  p-value: 5.095e-08
# Uji Normalitas Residual
shapiro.test(residuals(ols_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(ols_model)
## W = 0.94793, p-value = 0.1909
# Uji Homoskedastisitas
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(ols_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols_model
## BP = 13.307, df = 8, p-value = 0.1017
# Uji Linearitas
library(lmtest)
resettest(ols_model)
## 
##  RESET test
## 
## data:  ols_model
## RESET = 3.0695, df1 = 2, df2 = 16, p-value = 0.07442
# Uji Autokorelasi Spasial Residual
jabar_merged$residual_ols <- residuals(ols_model)
moran_res <- moran.test(
  jabar_merged$residual_ols,
  WL,
  zero.policy = TRUE
)
moran_res
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$residual_ols  
## weights: WL    
## 
## Moran I statistic standard deviate = 1.5723, p-value = 0.05794
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.18141415       -0.03846154        0.01955598
# ------------------------- CEK AUTOKORELASI SPASIAL DENGAN SEMUA X --------------------------

# ============================================================
# JP
# ============================================================
# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$JP,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$JP  
## weights: WL    
## 
## Moran I statistic standard deviate = 0.22309, p-value = 0.4117
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.008562506      -0.038461538       0.017962121
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$JP,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$JP 
## weights: WL   
## 
## Geary C statistic standard deviate = -1.6217, p-value = 0.9476
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        1.30372571        1.00000000        0.03507526
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$JP,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$JP, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$JP 
## weights: WL   
## 
## standard deviate = 2.0697, p-value = 0.01924
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       4.685361e-02       3.846154e-02       1.644126e-05
# ============================================================
# HIV
# ============================================================
# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$HIV,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$HIV  
## weights: WL    
## 
## Moran I statistic standard deviate = 2.3904, p-value = 0.008415
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.28189530       -0.03846154        0.01796086
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$HIV,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$HIV 
## weights: WL   
## 
## Geary C statistic standard deviate = 1.772, p-value = 0.03819
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.66809558        1.00000000        0.03508188
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$HIV,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$HIV, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$HIV 
## weights: WL   
## 
## standard deviate = 1.9444, p-value = 0.02593
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       4.911834e-02       3.846154e-02       3.004006e-05
# ============================================================
# ST
# ============================================================
# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$ST,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$ST  
## weights: WL    
## 
## Moran I statistic standard deviate = 1.6448, p-value = 0.05001
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.18586762       -0.03846154        0.01860240
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$ST,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$ST 
## weights: WL   
## 
## Geary C statistic standard deviate = 1.9452, p-value = 0.02588
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.65356903        1.00000000        0.03171895
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$ST,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$ST, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$ST 
## weights: WL   
## 
## standard deviate = 0.084846, p-value = 0.4662
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       3.886751e-02       3.846154e-02       2.289465e-05
# ============================================================
# DB
# ============================================================
# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$DB,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$DB  
## weights: WL    
## 
## Moran I statistic standard deviate = 2.0205, p-value = 0.02166
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.24113652       -0.03846154        0.01914872
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$DB,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$DB 
## weights: WL   
## 
## Geary C statistic standard deviate = 0.88116, p-value = 0.1891
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.85031867        1.00000000        0.02885517
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$DB,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$DB, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$DB 
## weights: WL   
## 
## standard deviate = 2.5117, p-value = 0.006007
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       4.974998e-02       3.846154e-02       2.019871e-05
# ============================================================
# Sanitasi
# ============================================================
# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$Sanitasi,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$Sanitasi  
## weights: WL    
## 
## Moran I statistic standard deviate = 0.97, p-value = 0.166
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.08120644       -0.03846154        0.01522002
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$Sanitasi,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$Sanitasi 
## weights: WL   
## 
## Geary C statistic standard deviate = 1.3936, p-value = 0.08172
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.69010149        1.00000000        0.04944916
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$Sanitasi,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$Sanitasi, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$Sanitasi 
## weights: WL   
## 
## standard deviate = 0.84938, p-value = 0.1978
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       3.922151e-02       3.846154e-02       8.005431e-07
# ============================================================
# Faskes
# ============================================================
# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$Faskes,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$Faskes  
## weights: WL    
## 
## Moran I statistic standard deviate = 1.565, p-value = 0.05879
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.17916284       -0.03846154        0.01933754
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$Faskes,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$Faskes 
## weights: WL   
## 
## Geary C statistic standard deviate = 0.4475, p-value = 0.3273
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.92529976        1.00000000        0.02786542
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$Faskes,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$Faskes, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$Faskes 
## weights: WL   
## 
## standard deviate = 2.6856, p-value = 0.00362
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       4.702572e-02       3.846154e-02       1.016933e-05
# ============================================================
# Miskin
# ============================================================
# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$Miskin,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$Miskin  
## weights: WL    
## 
## Moran I statistic standard deviate = -0.48851, p-value = 0.6874
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.10389258       -0.03846154        0.01794002
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$Miskin,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$Miskin 
## weights: WL   
## 
## Geary C statistic standard deviate = -1.9553, p-value = 0.9747
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        1.36680787        1.00000000        0.03519112
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$Miskin,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$Miskin, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$Miskin 
## weights: WL   
## 
## standard deviate = 2.7156, p-value = 0.003308
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       4.996325e-02       3.846154e-02       1.793915e-05
# ============================================================
# KP 
# ============================================================
# -----------Morans'I-----------
Global_Moran <- moran.test(
  jabar_merged$KP,
  WL,
  zero.policy = TRUE
)
Global_Moran
## 
##  Moran I test under randomisation
## 
## data:  jabar_merged$KP  
## weights: WL    
## 
## Moran I statistic standard deviate = 2.1464, p-value = 0.01592
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.25844583       -0.03846154        0.01913399
# -----------Geary's C-----------
Global_Geary <- geary.test(
  jabar_merged$KP,
  WL,
  zero.policy = TRUE
)
Global_Geary
## 
##  Geary C test under randomisation
## 
## data:  jabar_merged$KP 
## weights: WL   
## 
## Geary C statistic standard deviate = 2.3887, p-value = 0.008454
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##         0.5936931         1.0000000         0.0289324
# -----------Getis-Ord General G-----------
Global_Getis <- globalG.test(
  jabar_merged$KP,
  WL,
  zero.policy = TRUE
)
## Warning in globalG.test(jabar_merged$KP, WL, zero.policy = TRUE): Binary
## weights recommended (especially for distance bands)
Global_Getis
## 
##  Getis-Ord global G statistic
## 
## data:  jabar_merged$KP 
## weights: WL   
## 
## standard deviate = 0.38583, p-value = 0.3498
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       0.0422751479       0.0384615385       0.0000976989
################### UJI LM TEST #################

library(spdep)

LM_tests <- lm.LMtests(
  ols_model,
  WL,
  test = c(
    "LMerr",     # LM Error
    "LMlag",     # LM Lag
    "RLMerr",    # Robust LM Error
    "RLMlag",    # Robust LM Lag
    "SARMA"      # LM SARMA
  ),
  zero.policy = TRUE
)
## Please update scripts to use lm.RStests in place of lm.LMtests
LM_tests
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes +
## Miskin + KP, data = data)
## test weights: listw
## 
## RSerr = 1.4336, df = 1, p-value = 0.2312
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes +
## Miskin + KP, data = data)
## test weights: listw
## 
## RSlag = 6.0566, df = 1, p-value = 0.01385
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes +
## Miskin + KP, data = data)
## test weights: listw
## 
## adjRSerr = 0.13417, df = 1, p-value = 0.7141
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes +
## Miskin + KP, data = data)
## test weights: listw
## 
## adjRSlag = 4.7572, df = 1, p-value = 0.02918
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes +
## Miskin + KP, data = data)
## test weights: listw
## 
## SARMA = 6.1908, df = 2, p-value = 0.04526
lm_df <- data.frame(
  Test = names(LM_tests),
  Statistic = sapply(LM_tests, function(x) x$statistic),
  P_value   = sapply(LM_tests, function(x) x$p.value)
)

lm_df
##                       Test Statistic    P_value
## RSerr.RSerr          RSerr 1.4336091 0.23117634
## RSlag.RSlag          RSlag 6.0566457 0.01385406
## adjRSerr.adjRSerr adjRSerr 0.1341717 0.71414527
## adjRSlag.adjRSlag adjRSlag 4.7572084 0.02917583
## SARMA.SARMA          SARMA 6.1908174 0.04525651
# =============================================================
# ========  MODEL SPASIAL EKONOMETRIK LANJUTAN  ===============
# =============================================================

library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.3.3
## 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
library(lmtest)
library(car)

# Pastikan data & listw sudah benar
# Variabel dependen: Rata-Rata.Lama.Sekolah

# =============================================================
# 1. SPATIAL LAG MODEL (SLM)
# =============================================================
library(spdep)

slm_model <- lagsarlm(
  JK ~ 
    JP +
    HIV +
    ST +
    DB +
    Sanitasi +
    Faskes +
    Miskin +
    KP,
  data   = data,
  listw  = WL,
  method = "LU"
)
## Warning in lagsarlm(JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
##   reciprocal condition number = 3.69182e-20 - using numerical Hessian.
summary(slm_model)
## 
## Call:lagsarlm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + 
##     Miskin + KP, data = data, listw = WL, method = "LU")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2597.59 -1015.88  -123.34   676.86  4562.44 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -1.4247e+02  3.3308e+03 -0.0428 0.965882
## JP           5.1626e-04  4.6526e-04  1.1096 0.267163
## HIV          3.6152e+00  2.8959e+00  1.2484 0.211883
## ST          -6.7170e+01  1.1140e+02 -0.6030 0.546526
## DB          -8.6924e-02  3.7137e-02 -2.3406 0.019252
## Sanitasi    -3.7648e+01  3.0169e+01 -1.2479 0.212068
## Faskes       1.1251e+02  4.2412e+01  2.6528 0.007984
## Miskin       2.5133e+01  9.6586e+00  2.6021 0.009264
## KP           1.9326e-01  1.1953e-01  1.6168 0.105924
## 
## Rho: 0.15945, LR test value: 5.0785, p-value: 0.024224
## Approximate (numerical Hessian) standard error: 0.067493
##     z-value: 2.3624, p-value: 0.018155
## Wald statistic: 5.5812, p-value: 0.018155
## 
## Log likelihood: -238.3018 for lag model
## ML residual variance (sigma squared): 2697100, (sigma: 1642.3)
## Number of observations: 27 
## Number of parameters estimated: 11 
## AIC: 498.6, (AIC for lm: 501.68)
# 1. Multikolinearitas
vif(lm(JK ~ 
         JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
       data = data))
##        JP       HIV        ST        DB  Sanitasi    Faskes    Miskin        KP 
##  3.132431  8.096913  1.536464  4.041859  1.804266 13.726707  8.199660  2.726582
# 2. Normalitas residual
shapiro.test(residuals(slm_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(slm_model)
## W = 0.95478, p-value = 0.2793
# 3. Homoskedastisitas
bptest.Sarlm(slm_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 12.567, df = 8, p-value = 0.1276
# 4. Linearitas
resettest(lm(JK ~ 
               JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
             data = data))
## 
##  RESET test
## 
## data:  lm(JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,     data = data)
## RESET = 3.0695, df1 = 2, df2 = 16, p-value = 0.07442
# 5. Autokorelasi spasial residual
moran.test(residuals(slm_model), WL, zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  residuals(slm_model)  
## weights: WL    
## 
## Moran I statistic standard deviate = 0.30568, p-value = 0.3799
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.00346073       -0.03846154        0.01880910
# =============================================================
# 2. SPATIAL ERROR MODEL (SEM)
# =============================================================
sem_model <- errorsarlm(
  JK ~ 
    JP +
    HIV +
    ST +
    DB +
    Sanitasi +
    Faskes +
    Miskin +
    KP,
  data   = data,
  listw  = WL,
  method = "LU"
)
## Warning in errorsarlm(JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
##   reciprocal condition number = 7.04439e-20 - using numerical Hessian.
summary(sem_model)
## 
## Call:errorsarlm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + 
##     Miskin + KP, data = data, listw = WL, method = "LU")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2603.51 -1127.59  -266.65  1054.53  4938.86 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -2.3581e+02  3.5466e+03 -0.0665  0.94699
## JP           2.2607e-04  3.8531e-04  0.5867  0.55739
## HIV          4.4902e+00  3.1436e+00  1.4284  0.15318
## ST          -8.4348e+00  1.0697e+02 -0.0789  0.93715
## DB          -6.9746e-02  3.7816e-02 -1.8443  0.06513
## Sanitasi    -2.2612e+01  3.3025e+01 -0.6847  0.49354
## Faskes       9.6836e+01  4.4653e+01  2.1686  0.03011
## Miskin       2.7883e+01  1.0971e+01  2.5416  0.01104
## KP           2.4926e-01  1.3536e-01  1.8415  0.06555
## 
## Lambda: 0.45538, LR test value: 3.0422, p-value: 0.081127
## Approximate (numerical Hessian) standard error: 0.21188
##     z-value: 2.1492, p-value: 0.031618
## Wald statistic: 4.6191, p-value: 0.031618
## 
## Log likelihood: -239.32 for error model
## ML residual variance (sigma squared): 2762800, (sigma: 1662.2)
## Number of observations: 27 
## Number of parameters estimated: 11 
## AIC: 500.64, (AIC for lm: 501.68)
# 1. Multikolinearitas
vif(lm(JK ~ 
         JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
       data = data))
##        JP       HIV        ST        DB  Sanitasi    Faskes    Miskin        KP 
##  3.132431  8.096913  1.536464  4.041859  1.804266 13.726707  8.199660  2.726582
# 2. Normalitas residual
shapiro.test(residuals(sem_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(sem_model)
## W = 0.95329, p-value = 0.2574
# 3. Homoskedastisitas
bptest.Sarlm(sem_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 15.61, df = 8, p-value = 0.04831
# 4. Linearitas
resettest(lm(JK ~ 
               JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
             data = data))
## 
##  RESET test
## 
## data:  lm(JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,     data = data)
## RESET = 3.0695, df1 = 2, df2 = 16, p-value = 0.07442
# 5. Autokorelasi spasial residual
moran.test(residuals(sem_model), WL, zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  residuals(sem_model)  
## weights: WL    
## 
## Moran I statistic standard deviate = 0.50856, p-value = 0.3055
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.03080089       -0.03846154        0.01854863
# =============================================================
# 3. SPATIAL AUTOREGRESSIVE COMBINED MODEL (SAC)
# =============================================================
library(spdep)

sac_model <- sacsarlm(
  JK ~ 
    JP +
    HIV +
    ST +
    DB +
    Sanitasi +
    Faskes +
    Miskin +
    KP,
  data   = data,
  listw  = WL,
  method = "LU"
)

summary(sac_model)
## 
## Call:sacsarlm(formula = JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + 
##     Miskin + KP, data = data, listw = WL, method = "LU")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2601.66 -1013.28  -118.78   674.48  4576.85 
## 
## Type: sac 
## Coefficients: (numerical Hessian approximate standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -2.2069e+02  3.7125e+03 -0.0594  0.95260
## JP           5.0825e-04  4.8316e-04  1.0519  0.29284
## HIV          3.6807e+00  3.1749e+00  1.1593  0.24633
## ST          -6.3302e+01  1.3324e+02 -0.4751  0.63472
## DB          -8.6222e-02  3.9165e-02 -2.2015  0.02770
## Sanitasi    -3.6774e+01  3.4845e+01 -1.0553  0.29127
## Faskes       1.1151e+02  4.6730e+01  2.3862  0.01702
## Miskin       2.5321e+01  1.0405e+01  2.4335  0.01496
## KP           1.9451e-01  1.2251e-01  1.5877  0.11236
## 
## Rho: 0.15659
## Approximate (numerical Hessian) standard error: 0.085716
##     z-value: 1.8269, p-value: 0.067715
## Lambda: 0.01953
## Approximate (numerical Hessian) standard error: 0.35697
##     z-value: 0.054711, p-value: 0.95637
## 
## LR test value: 5.0803, p-value: 0.078853
## 
## Log likelihood: -238.3009 for sac model
## ML residual variance (sigma squared): 2697300, (sigma: 1642.4)
## Number of observations: 27 
## Number of parameters estimated: 12 
## AIC: 500.6, (AIC for lm: 501.68)
# ------------------ Uji Asumsi SDM ------------------

# 1. Multikolinearitas
library(car)

vif(lm(
  JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
  data = data
))
##        JP       HIV        ST        DB  Sanitasi    Faskes    Miskin        KP 
##  3.132431  8.096913  1.536464  4.041859  1.804266 13.726707  8.199660  2.726582
# 2. Normalitas residual
shapiro.test(residuals(sac_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(sac_model)
## W = 0.95474, p-value = 0.2787
# 3. Homoskedastisitas
library(spdep)

bptest.Sarlm(sac_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 12.658, df = 8, p-value = 0.1242
# 4. Linearitas
resettest(lm(JK ~ 
               JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,
             data = data))
## 
##  RESET test
## 
## data:  lm(JK ~ JP + HIV + ST + DB + Sanitasi + Faskes + Miskin + KP,     data = data)
## RESET = 3.0695, df1 = 2, df2 = 16, p-value = 0.07442
# 5. Autokorelasi spasial residual
moran.test(
  residuals(sac_model),
  WL,
  zero.policy = TRUE
)
## 
##  Moran I test under randomisation
## 
## data:  residuals(sac_model)  
## weights: WL    
## 
## Moran I statistic standard deviate = 0.29116, p-value = 0.3855
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.00145773       -0.03846154        0.01879702
# =============================================================
# 4. PERBANDINGAN MODEL
# =============================================================
data.frame(
  Model = c("OLS", "SLM", "SEM", "SAC"),
  LogLik = c(
    as.numeric(logLik(ols_model)),
    as.numeric(logLik(slm_model)),
    as.numeric(logLik(sem_model)),
    as.numeric(logLik(sac_model))
  ),
  AIC = c(
    AIC(ols_model),
    AIC(slm_model),
    AIC(sem_model),
    AIC(sac_model)
  )
)
##   Model    LogLik      AIC
## 1   OLS -240.8411 501.6821
## 2   SLM -238.3018 498.6036
## 3   SEM -239.3200 500.6399
## 4   SAC -238.3009 500.6018
# =============================================================
# ========  MODEL SPASIAL INTERPOLASI           ===============
# =============================================================

# 1. Transformasi Data Area → Titik (Centroid)
library(sf)

# Ambil centroid tiap kabupaten/kota
centroid_jabar <- st_centroid(jabar_merged)
## Warning: st_centroid assumes attributes are constant over geometries
# Cek hasil
plot(st_geometry(jabar_merged))
plot(st_geometry(centroid_jabar), add = TRUE, col = "red", pch = 19)

# 2. Persiapan Grid Prediksi (Permukaan Kontinu)
library(sp)

# Ubah ke sp object (dibutuhkan gstat)
centroid_sp <- as(centroid_jabar, "Spatial")

# Buat grid prediksi
grid <- st_make_grid(jabar_merged, cellsize = 0.05, what = "centers")
grid_sf <- st_sf(geometry = grid)
grid_sp <- as(grid_sf, "Spatial")

# 3. Interpolasi IDW 
library(sf)
library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.29
## 
## Attaching package: 'terra'
## The following object is masked from 'package:zoo':
## 
##     time<-
library(automap)
library(gstat)
## Warning: package 'gstat' was built under R version 4.3.3
idw_model <- gstat(
  formula = Rate_1k ~ 1,
  data = centroid_sp,
  set = list(idp = 2)
)

idw_pred <- predict(idw_model, grid_sp)
## [inverse distance weighted interpolation]
# Visualisasi IDW
spplot(
  idw_pred["var1.pred"],
  main = "Interpolasi Prevalensi TBC - IDW"
)

# Interpolasi Trend Surface (Polynomial Regression)
# Ambil data atribut
df_trend <- centroid_sp@data

head(data)
##        kab_kota    JK      JP HIV   ST    DB Sanitasi Faskes Miskin   KP
## 1       Bandung 14012 3753120 598 8.87 44200    96.79     79  239.9 2156
## 2 Bandung Barat  4569 1884190 153 4.95 15911    94.91     38  179.7 1468
## 3        Bekasi 15458 3273870 893 1.70 38666    97.00    102  204.5 2617
## 4         Bogor 29110 5682300 831 1.91 65751    88.77    129  446.8 1899
## 5        Ciamis  3477 1259230 110 2.52 14906    95.16     42   90.8  789
## 6       Cianjur 12197 2584990 293 2.51 15838    69.99     54  239.3  712
##    Rate_1k
## 1 3.733427
## 2 2.424915
## 3 4.721629
## 4 5.122926
## 5 2.761211
## 6 4.718393
# Tambahkan koordinat
coords <- coordinates(centroid_sp)
df_trend$x <- coords[,1]
df_trend$y <- coords[,2]

# Pastikan nama variabel BENAR
names(df_trend)
##  [1] "NAMOBJ"       "FCODE"        "REMARK"       "METADATA"     "SRS_ID"      
##  [6] "KDBBPS"       "KDCBPS"       "KDCPUM"       "KDEBPS"       "KDEPUM"      
## [11] "KDPBPS"       "KDPKAB"       "KDPPUM"       "LUASWH"       "TIPADM"      
## [16] "WADMKC"       "WADMKD"       "WADMKK"       "WADMPR"       "WIADKC"      
## [21] "WIADKK"       "WIADPR"       "WIADKD"       "SHAPE_Leng"   "SHAPE_Area"  
## [26] "JK"           "JP"           "HIV"          "ST"           "DB"          
## [31] "Sanitasi"     "Faskes"       "Miskin"       "KP"           "Rate_1k"     
## [36] "GiZ"          "Gi_cluster"   "Ii"           "Z.Ii"         "P.value"     
## [41] "cluster"      "residual_ols" "x"            "y"
trend_model <- lm(
  Rate_1k ~ poly(x, 2) + poly(y, 2),
  data = df_trend
)


grid_coords <- coordinates(grid_sp)

trend_pred <- predict(
  trend_model,
  newdata = data.frame(
    x = grid_coords[,1],
    y = grid_coords[,2]
  )
)

grid_sp$trend <- trend_pred


# Visualisasi Trend Surface 
spplot(
  grid_sp["trend"],
  main = "Interpolasi Prevalensi TBC - Trend Surface"
)

# 5. Interpolasi Ordinary Kriging

# a. Variogram Empiris
vgm_emp <- variogram(Rate_1k ~ 1, centroid_sp)
plot(vgm_emp)

# b. Fitting Variogram
vgm_fit <- fit.variogram(vgm_emp, model = vgm("Sph"))
## Warning in fit.variogram(vgm_emp, model = vgm("Sph")): No convergence after 200
## iterations: try different initial values?
plot(vgm_emp, vgm_fit)

# c. Kriging
kriging_model <- gstat(
  formula = Rate_1k ~ 1,
  data = centroid_sp,
  model = vgm_fit
)

kriging_pred <- predict(kriging_model, grid_sp)
## [using ordinary kriging]
# Visualisasi Kriging
spplot(
  kriging_pred["var1.pred"],
  main = "Interpolasi Prevalensi TBC - Ordinary Kriging"
)

# =============================================================
# ======== EVALUASI MODEL SPASIAL INTERPOLASI   ===============
# =============================================================

library(gstat)
library(sp)

# ----------------------- Cross-Validation IDW ----------------------- #

# Cross-validation IDW
library(gstat)

idw_model <- gstat(
  formula = Rate_1k ~ 1,
  data = centroid_sp,
  set = list(idp = 2)
)

cv_idw <- gstat.cv(
  object = idw_model,
  nfold = nrow(centroid_sp)
)
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
# Error metrics IDW
ME_idw   <- mean(cv_idw$residual)
MAE_idw  <- mean(abs(cv_idw$residual))
RMSE_idw <- sqrt(mean(cv_idw$residual^2))

ME_idw
## [1] -0.009213774
MAE_idw
## [1] 3.997565
RMSE_idw
## [1] 6.895552
# ----------------------- Cross-Validation Trend Surface ----------------------- #
# Data untuk trend surface
df_trend <- centroid_sp@data
coords <- coordinates(centroid_sp)
df_trend$x <- coords[,1]
df_trend$y <- coords[,2]

# Leave-One-Out CV
pred_trend <- numeric(nrow(df_trend))

for(i in 1:nrow(df_trend)){
  train <- df_trend[-i, ]
  test  <- df_trend[i, ]
  
  model <- lm(
    Rate_1k ~ poly(x, 2) + poly(y, 2),
    data = train
  )
  
  pred_trend[i] <- predict(model, newdata = test)
}

# Residual
res_trend <- df_trend$Rate_1k - pred_trend

# Error metrics Trend Surface
ME_trend   <- mean(res_trend)
MAE_trend  <- mean(abs(res_trend))
RMSE_trend <- sqrt(mean(res_trend^2))

ME_trend
## [1] 0.1994893
MAE_trend
## [1] 4.284163
RMSE_trend
## [1] 6.929764
# ----------------------- Cross-Validation Ordinary Kriging ----------------------- #
# Variogram
vgm_emp <- variogram(Rate_1k ~ 1, centroid_sp)
vgm_fit <- fit.variogram(vgm_emp, vgm("Sph"))
## Warning in fit.variogram(vgm_emp, vgm("Sph")): No convergence after 200
## iterations: try different initial values?
# Cross-validation Kriging
cv_ok <- krige.cv(
  formula = Rate_1k ~ 1,
  locations = centroid_sp,
  model = vgm_fit,
  nfold = nrow(centroid_sp)
)

# Error metrics Kriging
ME_ok   <- mean(cv_ok$residual)
MAE_ok  <- mean(abs(cv_ok$residual))
RMSE_ok <- sqrt(mean(cv_ok$residual^2))

ME_ok
## [1] 0.007087704
MAE_ok
## [1] 3.466216
RMSE_ok
## [1] 6.280008
########## Ringkasan

cv_result <- data.frame(
  Model = c("IDW", "Trend Surface", "Ordinary Kriging"),
  ME    = c(ME_idw, ME_trend, ME_ok),
  MAE   = c(MAE_idw, MAE_trend, MAE_ok),
  RMSE  = c(RMSE_idw, RMSE_trend, RMSE_ok)
)

cv_result
##              Model           ME      MAE     RMSE
## 1              IDW -0.009213774 3.997565 6.895552
## 2    Trend Surface  0.199489306 4.284163 6.929764
## 3 Ordinary Kriging  0.007087704 3.466216 6.280008

```