1 LOAD LIBRARY

Memuat seluruh package R yang dibutuhkan untuk analisis spasial dan pemodelan.

library(readr)
library(tmap)
library(sf)
library(spdep)
library(spatialreg)
library(dplyr)
library(ggplot2)
library(lmtest)
library(car)
library(nortest)
library(DT)

2 LOAD DATA

Dua sumber data digunakan: 1. Shapefile (ntt_shp): Data geometri batas wilayah kabupaten/kota Provinsi NTT. 2. Data atribut 22 wilayah dengan 18 variabel, meliputi: Variabel respons: Stunting (prevalensi stunting, %) Variabel prediktor: Jumlah Faskes, ASI, PDRB, IPM, Jumlah Tenaga Kesehatan, Sanitasi Layak, Kemiskinan (Miskin), Gizi Buruk, Wasting, Desa Posbindu, Posyandu, RLS, Hunian Layak, Air Minum Layak, Gini Ratio, Pengeluaran

# Data Shapefile
ntt_shp <- st_read("D:/stat/6th Semester/SPASIAL/PROJEK/ntt.shp")
## Reading layer `ntt' from data source `D:\stat\6th Semester\SPASIAL\PROJEK\ntt.shp' using driver `ESRI Shapefile'
## Simple feature collection with 22 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 118.927 ymin: -11.00762 xmax: 125.1933 ymax: -7.778658
## Geodetic CRS:  WGS 84
# Data stunting
data_ntt <- read_csv("D:/stat/6th Semester/SPASIAL/PROJEK/NTT.STUNTING.csv", show_col_types = FALSE)
data_ntt
## # A tibble: 22 × 16
##    WADMKK  Stunting Jumlah.Faskes   ASI  PDRB   IPM Jumlah.Tenaga Sanitasi.Layak
##    <chr>      <dbl>         <dbl> <dbl> <dbl> <dbl>         <dbl>          <dbl>
##  1 Sumba …     23.1            13 100   10769  68.5          1028           62.8
##  2 Sumba …     17.7            29  98.6 17518  71.0          1539           70.1
##  3 Kupang      15.0            29  97.5 14732  69.0          1552           77.1
##  4 Timor …     41.4            40  96.3 11467  66.9          1786           71.1
##  5 Timor …     30.0            30 100   12052  67.9          1931           84.3
##  6 Belu        29.2            21  98.4 14710  68.9          1297           86.9
##  7 Alor        22.6            31  95.0 10471  68.4          1850           88.4
##  8 Lembata     10.5            15  96.6  9232  69.7          1082           92.6
##  9 Flores…     18.7            24 100   13436  70.6          1836           95.1
## 10 Sikka       12.6            28  97.5 11456  70.6          2162           87.1
## # ℹ 12 more rows
## # ℹ 8 more variables: Miskin <dbl>, Gizi.Buruk <dbl>, Wasting <dbl>,
## #   Desa.Posbindu <dbl>, Posyandu <dbl>, RLS <dbl>, Hunian.Layak <dbl>,
## #   Air.Minum.Layak <dbl>

3 JOIN SHAPEFILE DAN DATA KASUS

Menggabungkan shapefile (ntt_shp) dengan data stunting (data_ntt) berdasarkan kolom kunci WADMKK (nama kabupaten/kota)

ntt_merge <- ntt_shp %>%
  left_join(data_ntt, by = "WADMKK") %>%
  st_make_valid() %>%
  filter(!is.na(Stunting))

cat("Jumlah kabupaten/kota:", nrow(ntt_merge), "\n")
## Jumlah kabupaten/kota: 22
cat("Variabel:", paste(names(ntt_merge), collapse = ", "), "\n")
## Variabel: KDPKAB, KDPPUM, WADMKK, WADMPR, METADATA, UPDATED, Stunting, Jumlah.Faskes, ASI, PDRB, IPM, Jumlah.Tenaga, Sanitasi.Layak, Miskin, Gizi.Buruk, Wasting, Desa.Posbindu, Posyandu, RLS, Hunian.Layak, Air.Minum.Layak, geometry

4 VISUALISASI PETA SHP

tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
q_breaks <- quantile(ntt_merge$Stunting, probs = seq(0, 1, 0.2), na.rm = TRUE)

tm_shape(ntt_merge) +
  tm_polygons(
    col          = "Stunting",
    style        = "fixed",
    breaks       = q_breaks,
    palette      = "YlOrRd",
    colorNA      = "gray80",
    border.col   = "black",
    border.alpha = 0.8,
    title        = "Stunting (%)"
  ) +
  tm_text("WADMKK", size = 0.45, col = "black") +
  tm_compass(position = c("right", "top")) +
  tm_scalebar(position = c("left", "bottom")) +
  tm_layout(
    main.title     = "Peta Sebaran Stunting Provinsi NTT",
    main.title.size = 1,
    legend.outside = TRUE,
    frame          = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values'),
##   'colorNA' (rename to 'value.na') to 'tm_scale_intervals(<HERE>)'
## For small multiples, specify a 'tm_scale_' for each multiple, and put them in a
## list: 'fill'.scale = list(<scale1>, <scale2>, ...)'[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: use `col_alpha` instead of `border.alpha`.[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
## "brewer.yl_or_rd"

Peta menunjukkan distribusi prevalensi stunting di 22 kabupaten/kota Provinsi Nusa Tenggara Timur (NTT) yang dibagi menjadi 5 kelas berdasarkan kuantil.Prevalensi stunting sangat tinggi berada pada wilayah Sumba Barat Daya, Belu, Timor Tengah Utara, Malaka, Timor Tengah Selatan

5 AUTOKORELASI SPASIAL

5.1 1. Membentuk Matriks Spasial

Menghitung koordinat sentroid tiap kabupaten/kota sebagai titik acuan. Mencoba nilai k = 1 hingga 10 untuk metode k-Nearest Neighbor (kNN), menghitung Moran’s I pada tiap nilai k untuk menemukan k yang menghasilkan autokorelasi spasial tertinggi.

coords <- st_coordinates(st_centroid(ntt_merge))
## Warning: st_centroid assumes attributes are constant over geometries
# Pemilihan k optimal
k_values <- 1:10

moran_I <- sapply(k_values, function(k) {
  nb_k <- knn2nb(knearneigh(coords, k = k))
  
  # Skip jika ada sub-graph (wilayah terisolasi)
  if (n.comp.nb(nb_k)$nc > 1) return(NA)
  
  lw_k <- nb2listw(nb_k, style = "W")
  tryCatch(
    moran.test(ntt_merge$Stunting, lw_k)$estimate[1],
    warning = function(w) NA,
    error   = function(e) NA
  )
})
## Warning in knn2nb(knearneigh(coords, k = k)): neighbour object has 6 sub-graphs
## Warning in knearneigh(coords, k = k): k greater than one-third of the number of
## data points
## Warning in knearneigh(coords, k = k): k greater than one-third of the number of
## data points
## Warning in knearneigh(coords, k = k): k greater than one-third of the number of
## data points
k_df <- data.frame(k = k_values, Moran_I = round(moran_I, 4))
cat("Moran's I per nilai k:\n")
## Moran's I per nilai k:
print(k_df)
##     k Moran_I
## 1   1      NA
## 2   2  0.3827
## 3   3  0.4357
## 4   4  0.3258
## 5   5  0.3043
## 6   6  0.2173
## 7   7  0.1672
## 8   8  0.1307
## 9   9  0.0919
## 10 10  0.0672
k_optimal <- k_values[which.max(moran_I)]
cat("\nk Optimal:", k_optimal, "\n")
## 
## k Optimal: 3

k = 3 dipilih sebagai k optimal karena menghasilkan Moran’s I global tertinggi (0.4357), Nilai NA pada k = 1 terjadi karena menghasilkan sub-graph (wilayah terisolasi), sehingga matriks bobot tidak valid.

# Bobot spasial kNN dengan k optimal
nb_knn  <- knn2nb(knearneigh(coords, k = k_optimal))
lw_knn  <- nb2listw(nb_knn, style = "W", zero.policy = TRUE)

5.2 Visualisasi Jaringan Tetangga KNN

Menggambar peta batas wilayah NTT beserta garis koneksi kNN (k=3), di mana setiap kabupaten/kota dihubungkan ke 3 tetangga terdekatnya.

plot(st_geometry(ntt_merge), border = "gray70", main = paste0("Jaringan Tetangga KNN (k=", k_optimal, ")"))
plot(nb_knn, coords, add = TRUE, col = "steelblue", lwd = 1.2)
points(coords, pch = 16, col = "red", cex = 0.8)

Peta jaringan menunjukkan titik-titik merah (sentroid kabupaten) yang terhubung oleh garis biru yang enggambarkan struktur ketetanggaan spasial

5.3 2. Moran’s I Global

Pengujian autokorelasi spasial menggunakan Indeks Moran’s I.

  • H₀: ρ = 0 (tidak ada autokorelasi spasial)
  • H₁: ρ ≠ 0 (ada autokorelasi spasial)
moran_global <- moran.test(ntt_merge$Stunting, lw_knn, zero.policy = TRUE)
print(moran_global)
## 
##  Moran I test under randomisation
## 
## data:  ntt_merge$Stunting  
## weights: lw_knn    
## 
## Moran I statistic standard deviate = 3.277, p-value = 0.0005245
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.43572115       -0.04761905        0.02175405

Kesimpulan: p-value < 0.05 maka H ₀ditolak . Terdapat autokorelasi spasial positif yang signifikan pada prevalensi stunting di NTT. Artinya, wilayah dengan prevalensi stunting tinggi cenderung bertetangga dengan wilayah yang juga memiliki prevalensi stunting tinggi (dan sebaliknya).

5.4 3. Moran’s I Lokal (LISA)

Menghitung indeks Moran lokal (LISA) untuk mengidentifikasi kluster spasial di tingkat kabupaten/kota.

local_moran <- localmoran(ntt_merge$Stunting, lw_knn, zero.policy = TRUE)

# Tambahkan hasil LISA ke data
ntt_merge$Ii   <- local_moran[, 1]   # Indeks Moran lokal
ntt_merge$P_Ii <- local_moran[, 5]   # p-value

# Visualisasi LISA
tm_shape(ntt_merge) +
  tm_polygons(
    col      = "Ii",
    palette  = "RdBu",
    midpoint = 0,
    title    = "Indeks Moran Lokal (Ii)"
  ) +
  tm_layout(
    main.title     = "LISA Stunting NTT",
    legend.outside = TRUE,
    frame          = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'midpoint', 'palette' (rename to 'values') to
## fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdBu" is named
## "brewer.rd_bu"
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.

Nilai Ii positif: wilayah tersebut dan tetangganya sama-sama tinggi (High-High) atau sama-sama rendah (Low-Low) Nilai Ii negatif: wilayah berbeda dengan tetangganya (High-Low atau Low-High)

5.5 4. Geary’s C

Uji Geary’s C sebagai pelengkap Moran’s I untuk mengkonfirmasi autokorelasi spasial.

geary.test(ntt_merge$Stunting, lw_knn)
## 
##  Geary C test under randomisation
## 
## data:  ntt_merge$Stunting 
## weights: lw_knn   
## 
## Geary C statistic standard deviate = 2.2927, p-value = 0.01093
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.64396873        1.00000000        0.02411413

Interpretasi: Nilai Geary’s C = 0.64 < 1 menunjukkan autokorelasi spasial positif (nilai mendekati 0 = autokorelasi positif kuat; mendekati 2 = negatif). p-value < 0.05 → signifikan. Hasil ini konsisten dengan Moran’s I global.

5.6 5. Getis Ord Gi*

Menghitung statistik Gi* untuk mengidentifikasi hot spot (konsentrasi nilai tinggi) dan cold spot (konsentrasi nilai rendah) prevalensi stunting. Wilayah diklasifikasikan berdasarkan Z-score

gi_star <- localG(ntt_merge$Stunting, lw_knn, zero.policy = TRUE)
ntt_merge$GiZ <- as.numeric(gi_star)

ntt_merge$Gi_cluster <- cut(
  ntt_merge$GiZ,
  breaks = c(-Inf, -2.58, -1.96, 1.96, 2.58, Inf),
  labels = c("Cold Spot 99%","Cold Spot 95%","Not Significant",
             "Hot Spot 95%","Hot Spot 99%")
)

tm_shape(ntt_merge) +
  tm_polygons(
    col     = "Gi_cluster",
    palette = c("blue4","lightskyblue","gray85","orange","red3"),
    title   = "Getis-Ord Gi*"
  ) +
  tm_text("WADMKK", size = 0.4) +
  tm_compass(position = c("right","top")) +
  tm_scalebar(position = c("left","bottom")) +
  tm_layout(
    main.title     = "Peta Hotspot Stunting NTT (Getis-Ord Gi*)",
    legend.outside = TRUE,
    frame          = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

Wilayah Belu, Timor Tengah Utara (TTU), dan Malaka menjadi hotspot paling kuat (99%), yang artinya ketiga kabupaten ini memiliki prevalensi stunting yang tinggi sekaligus saling berdekatan, sehingga membentuk kluster.

6 Model Spasial

6.1 1. Model OLS

Variabel respons: Y = presentase balita stunting (%) Variabel prediktor: X1 = Indeks Pembangunan Manusia X2 = Pemberian ASI pada bayi umur 0-23 bulan (Baduta) (%) X3 = persentase penduduk miskin (%) X4 = Jumlah Desa yang memiliki Posbindu tiap Kabupaten/Kota X5 = Persentase balita Wasting (%)

ols_model <- lm(Stunting ~  IPM+ASI+Miskin + Desa.Posbindu+ Wasting,  data = ntt_merge)
summary(ols_model)
## 
## Call:
## lm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu + 
##     Wasting, data = ntt_merge)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0226  -4.7669  -0.5528   5.1405  11.0871 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   132.99514  125.50310   1.060   0.3050  
## IPM            -0.58065    0.66332  -0.875   0.3943  
## ASI            -0.98635    0.86843  -1.136   0.2728  
## Miskin          1.16080    1.58513   0.732   0.4746  
## Desa.Posbindu   0.02801    0.02896   0.967   0.3479  
## Wasting         1.42265    0.61654   2.307   0.0347 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.504 on 16 degrees of freedom
## Multiple R-squared:  0.4723, Adjusted R-squared:  0.3074 
## F-statistic: 2.864 on 5 and 16 DF,  p-value: 0.04935
AIC(ols_model)
## [1] 158.1067

6.1.1 Uji Asumsi OLS

6.1.1.1 Uji Normalitas Residual (Kolmogorov-Smirnov)

  • H₀: Residual berdistribusi normal
  • H₁: Residual tidak berdistribusi normal
ks_test <- ks.test(residuals(ols_model), "pnorm",
                   mean = mean(residuals(ols_model)),
                   sd   = sd(residuals(ols_model)))
print(ks_test)
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  residuals(ols_model)
## D = 0.13365, p-value = 0.7789
## alternative hypothesis: two-sided
# QQ Plot
qqnorm(residuals(ols_model), main = "Q-Q Plot Residual OLS")
qqline(residuals(ols_model), col = "red", lwd = 2)

p-value = 0.7789 > 0.05 sehingga H₀ tidak ditolak. Dengan demikian Residual berdistribusi normal.

6.1.1.2 Uji Homoskedastisitas / Residual Identik (Breusch-Pagan)

  • H₀: σ₁² = σ₂² = … = σₙ² (residual homogen/identik)
  • H₁: minimal ada satu σᵢ² ≠ σ² (heteroskedastisitas)
bp_test <- bptest(ols_model)
print(bp_test)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols_model
## BP = 7.0471, df = 5, p-value = 0.2172

p-value = 0.2172 > 0.05 sehingga H₀ tidak ditolak. Dengan demikian Variansi residual homogen

6.1.1.3 Uji Bebas Autokorelasi Residual (Durbin-Watson)

  • H₀: ρ = 0 (tidak ada autokorelasi)
  • H₁: ρ ≠ 0 (ada autokorelasi)
library(lmtest)
dw_test <- dwtest(ols_model)
print(dw_test)
## 
##  Durbin-Watson test
## 
## data:  ols_model
## DW = 1.2282, p-value = 0.009828
## alternative hypothesis: true autocorrelation is greater than 0

p-value = 0.009828 < 0.05 sehingga H₀ ditolak artinya terdapat autokorelasi positif pada residual OLS

6.1.1.4 Uji Multikolinearitas (VIF)

Nilai VIF > 10 mengindikasikan multikolinearitas.

vif_result <- vif(ols_model)
print(round(vif_result, 3))
##           IPM           ASI        Miskin Desa.Posbindu       Wasting 
##         2.732         1.767         2.729         1.173         1.093

Semua nilai VIF < 10 sehingga tidak ada multikolinearitas. Dengan demikian Variabel prediktor tidak saling berkorelasi

6.1.2 Peta Residual OLS

ntt_merge$resid_ols <- residuals(ols_model)

tm_shape(ntt_merge) +
  tm_polygons(
    col      = "resid_ols",
    palette  = "RdBu",
    midpoint = 0,
    title    = "Residual OLS"
  ) +
  tm_text("WADMKK", size = 0.4) +
  tm_layout(
    main.title     = "Peta Residual Model OLS",
    legend.outside = TRUE,
    frame          = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'midpoint', 'palette' (rename to 'values') to
## fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdBu" is named
## "brewer.rd_bu"
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.

Peta choropleth yang memetakan nilai residual OLS per kabupaten. Pola spasial yang terlihat pada peta residual memperkuat dugaan adanya autokorelasi spasial yang belum tertangkap oleh OLS.

6.1.3 Moran’s I pada Residual OLS

moran_resid <- moran.test(residuals(ols_model), lw_knn, zero.policy = TRUE)
print(moran_resid)
## 
##  Moran I test under randomisation
## 
## data:  residuals(ols_model)  
## weights: lw_knn    
## 
## Moran I statistic standard deviate = 1.8177, p-value = 0.03455
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.23196339       -0.04761905        0.02365744

Moran I = 0.23196 > E(I) = -0.0476 serta p-value = 0.003455 < 0.05 yang berarti Residual OLS secara signifikan mengandung autokorelasi spasial

6.2 Uji Efek Spasial (Lagrange Multiplier)

Berdasarkan skema analisis Anselin (2005), pengujian LM digunakan untuk menentukan model spasial yang tepat.

  • LM-Lag: Menguji dependensi spasial lag (→ SAR)
  • LM-Error: Menguji dependensi spasial error (→ SEM)
lm_tests <- lm.RStests(
  model       = ols_model,
  listw       = lw_knn,
  test        = "all"
)
print(lm_tests)
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu +
## Wasting, data = ntt_merge)
## test weights: lw_knn
## 
## RSerr = 1.9863, df = 1, p-value = 0.1587
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu +
## Wasting, data = ntt_merge)
## test weights: lw_knn
## 
## RSlag = 7.9702, df = 1, p-value = 0.004755
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu +
## Wasting, data = ntt_merge)
## test weights: lw_knn
## 
## adjRSerr = 3.0381, df = 1, p-value = 0.08133
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu +
## Wasting, data = ntt_merge)
## test weights: lw_knn
## 
## adjRSlag = 9.022, df = 1, p-value = 0.002668
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu +
## Wasting, data = ntt_merge)
## test weights: lw_knn
## 
## SARMA = 11.008, df = 2, p-value = 0.00407

Keputusan: Jika LM-Lag signifikan → pilih SAR/SDM. Jika LM-Error signifikan → pilih SEM/SDEM. Jika keduanya signifikan → lihat Robust LM, pilih yang lebih signifikan.

Uji Statistik p-value Signifikansi
LM-Error (RSerr) 1.9863 0.1587 Tidak
LM-Lag (RSlag) 7.9702 0.0048 Ya
Robust LM-Error 3.0381 0.0813 Tidak
Robust LM-Lag 9.0220 0.0027 Ya
SARMA 11.008 0.0041 Ya

6.3 2. SAR

SAR digunakan jika ρ ≠ 0 dan λ = 0 nilai Y di suatu wilayah dipengaruhi oleh nilai Y di wilayah tetangga (efek lag spasial pada variabel respons) Persamaan: Y = ρWY + Xβ + ε .

sar_model <- lagsarlm(Stunting ~ IPM+ASI+Miskin + Desa.Posbindu+ Wasting , data = ntt_merge,listw = lw_knn)

summary(sar_model, Nagelkerke = TRUE)
## 
## Call:lagsarlm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu + 
##     Wasting, data = ntt_merge, listw = lw_knn)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -10.35296  -4.02871  -0.84964   4.85421  10.50671 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)   140.389024  91.598428  1.5327  0.12536
## IPM            -0.593940   0.483882 -1.2274  0.21965
## ASI            -1.103955   0.635570 -1.7370  0.08240
## Miskin          0.912847   1.157077  0.7889  0.43016
## Desa.Posbindu   0.040558   0.021164  1.9163  0.05532
## Wasting         1.121314   0.478383  2.3440  0.01908
## 
## Rho: 0.40639, LR test value: 5.7381, p-value: 0.016601
## Asymptotic standard error: 0.16805
##     z-value: 2.4184, p-value: 0.01559
## Wald statistic: 5.8485, p-value: 0.01559
## 
## Log likelihood: -69.18429 for lag model
## ML residual variance (sigma squared): 29.935, (sigma: 5.4713)
## Nagelkerke pseudo-R-squared: 0.59347 
## Number of observations: 22 
## Number of parameters estimated: 8 
## AIC: 154.37, (AIC for lm: 158.11)
## LM test for residual autocorrelation
## test value: 3.6852, p-value: 0.054896

Rho (ρ) = 0.4064 ≠ 0 dengan p-value 0.01559 < 0.05 artinya secara signifikan prevalensi stunting di suatu kabupaten memang dipengaruhi oleh kabupaten-kabupaten di sekitarnya.

ntt_merge$resid_sar <- residuals(sar_model)

tm_shape(ntt_merge) +
  tm_polygons(
    col      = "resid_sar",
    palette  = "RdBu",
    midpoint = 0,
    title    = "Residual SAR"
  ) +
  tm_layout(
    main.title = "Peta Residual Model SAR",
    legend.outside = TRUE, frame = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'midpoint', 'palette' (rename to 'values') to
## fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdBu" is named
## "brewer.rd_bu"
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.

6.4 4. SEM

SEM digunakan jika λ ≠ 0 dan ρ = 0. Persamaan: Y = Xβ + u, di mana u = λWu + ε

#SEM
sem_model <- errorsarlm(Stunting ~  IPM+ASI+Miskin + Desa.Posbindu+ Wasting, data = ntt_merge,listw = lw_knn)

summary(sem_model,Nagelkerke = TRUE)
## 
## Call:errorsarlm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu + 
##     Wasting, data = ntt_merge, listw = lw_knn)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.5634 -4.1305 -1.2062  3.3311 12.1200 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)   134.856545  97.271391  1.3864  0.16563
## IPM            -0.506155   0.517395 -0.9783  0.32794
## ASI            -1.010068   0.641620 -1.5742  0.11543
## Miskin          0.242134   1.326829  0.1825  0.85520
## Desa.Posbindu   0.041909   0.020885  2.0067  0.04479
## Wasting         1.396634   0.616166  2.2667  0.02341
## 
## Lambda: 0.44825, LR test value: 3.0217, p-value: 0.082159
## Asymptotic standard error: 0.19152
##     z-value: 2.3405, p-value: 0.019258
## Wald statistic: 5.478, p-value: 0.019258
## 
## Log likelihood: -70.5425 for error model
## ML residual variance (sigma squared): 33.426, (sigma: 5.7815)
## Nagelkerke pseudo-R-squared: 0.54004 
## Number of observations: 22 
## Number of parameters estimated: 8 
## AIC: 157.09, (AIC for lm: 158.11)

parameter λ = 0.4483, p-value (Wald) = 0.0193 < 0.05, artinya efek error spasial tidak signifikan secara statistik sehingga model SEM tidak cocok digunakan.

ntt_merge$resid_sem <- residuals(sem_model)

tm_shape(ntt_merge) +
  tm_polygons(
    col      = "resid_sem",
    palette  = "RdBu",
    midpoint = 0,
    title    = "Residual SEM"
  ) +
  tm_layout(
    main.title = "Peta Residual Model SEM",
    legend.outside = TRUE, frame = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'midpoint', 'palette' (rename to 'values') to
## fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdBu" is named
## "brewer.rd_bu"
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.

6.5 5. SDM

SDM adalah model paling umum — menggabungkan lag Y dan lag X. Persamaan: Y = ρWY + Xβ + WXθ + ε

SDM dapat direduksi menjadi SAR (jika θ=0) atau SEM (jika θ+ρβ=0).

#SDM
sdm_model = lagsarlm(Stunting ~  IPM+ASI+Miskin + Desa.Posbindu+ Wasting, data = ntt_merge,listw = lw_knn, type = "mixed")

summary(sdm_model,Nagelkerke = TRUE)
## 
## Call:lagsarlm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu + 
##     Wasting, data = ntt_merge, listw = lw_knn, type = "mixed")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -7.20810 -2.33912  0.49488  2.82512  5.22901 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                      Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)       769.9306769 311.5973073  2.4709 0.0134768
## IPM                -1.3192828   0.3839771 -3.4358 0.0005907
## ASI                -2.4962294   0.7456593 -3.3477 0.0008149
## Miskin              0.1631982   0.9074183  0.1798 0.8572712
## Desa.Posbindu       0.0430371   0.0174241  2.4700 0.0135124
## Wasting             1.6499568   0.4823370  3.4208 0.0006245
## lag.IPM            -2.0283779   0.9806188 -2.0685 0.0385961
## lag.ASI            -3.2413130   2.0143697 -1.6091 0.1075955
## lag.Miskin          2.7887980   1.9404981  1.4372 0.1506737
## lag.Desa.Posbindu   0.0037255   0.0375991  0.0991 0.9210719
## lag.Wasting         0.8251666   0.6916406  1.1931 0.2328471
## 
## Rho: -0.29532, LR test value: 1.9059, p-value: 0.16742
## Asymptotic standard error: 0.21133
##     z-value: -1.3974, p-value: 0.16229
## Wald statistic: 1.9527, p-value: 0.16229
## 
## Log likelihood: -59.4379 for mixed model
## ML residual variance (sigma squared): 12.737, (sigma: 3.5689)
## Nagelkerke pseudo-R-squared: 0.83239 
## Number of observations: 22 
## Number of parameters estimated: 13 
## AIC: 144.88, (AIC for lm: 144.78)
## LM test for residual autocorrelation
## test value: 2.6923, p-value: 0.10083
Variabel Estimasi p-value Signifikan?
(Intercept) 769.931 0.013 Ya
IPM −1.319 0.001 Ya
ASI −2.496 0.001 Ya
Miskin 0.163 0.857 Tidak
Desa.Posbindu 0.043 0.014 Ya
Wasting 1.650 0.001 Ya
lag.IPM −2.028 0.039 Ya
lag.ASI −3.241 0.108 Tidak
lag.Miskin 2.789 0.151 Tidak
lag.Desa.Posbindu 0.004 0.921 Tidak
lag.Wasting 0.825 0.233 Tidak
Rho (ρ) −0.295 0.162 Tidak

Rho (ρ) tidak signifikan artinya stunting di suatu kabupaten/kota tidak terbukti dipengaruhi oleh stunting di kabupaten/kota tetangganya.

ntt_merge$resid_sdm <- residuals(sdm_model)

tm_shape(ntt_merge) +
  tm_polygons(
    col      = "resid_sdm",
    palette  = "RdBu",
    midpoint = 0,
    title    = "Residual SDM"
  ) +
  tm_layout(
    main.title = "Peta Residual Model SDM",
    legend.outside = TRUE, frame = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'midpoint', 'palette' (rename to 'values') to
## fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdBu" is named
## "brewer.rd_bu"
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.

## 6. SDEM

SDEM menggabungkan error spasial dan lag X. Persamaan: Y = Xβ + WXθ + u, di mana u = λWu + ε

#SDEM
sdem_model = errorsarlm(Stunting ~ IPM+ASI+Miskin + Desa.Posbindu+ Wasting,data   = ntt_merge,listw  = lw_knn,etype  = "emixed"
)
summary(sdem_model,Nagelkerke = TRUE)
## 
## Call:errorsarlm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu + 
##     Wasting, data = ntt_merge, listw = lw_knn, etype = "emixed")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.65757 -1.72022 -0.04561  1.47743  4.60735 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                     Estimate Std. Error z value  Pr(>|z|)
## (Intercept)       984.525900 196.341570  5.0144 5.321e-07
## IPM                -1.289156   0.405665 -3.1779 0.0014835
## ASI                -2.578740   0.635416 -4.0584 4.942e-05
## Miskin              0.076708   1.015903  0.0755 0.9398109
## Desa.Posbindu       0.042084   0.018645  2.2572 0.0239964
## Wasting             1.406862   0.535961  2.6249 0.0086666
## lag.IPM            -2.268566   0.833935 -2.7203 0.0065220
## lag.ASI            -5.124372   1.338236 -3.8292 0.0001286
## lag.Miskin          1.131869   1.962902  0.5766 0.5641891
## lag.Desa.Posbindu   0.043806   0.022896  1.9132 0.0557201
## lag.Wasting         0.425853   0.592654  0.7186 0.4724167
## 
## Lambda: -1.4512, LR test value: 9.9019, p-value: 0.001651
## Asymptotic standard error: 0.20996
##     z-value: -6.9119, p-value: 4.7806e-12
## Wald statistic: 47.775, p-value: 4.7807e-12
## 
## Log likelihood: -55.43988 for error model
## ML residual variance (sigma squared): 5.7559, (sigma: 2.3991)
## Nagelkerke pseudo-R-squared: 0.88347 
## Number of observations: 22 
## Number of parameters estimated: 13 
## AIC: 136.88, (AIC for lm: 144.78)
Variabel Estimasi p-value Signifikan
(Intercept) 984.526 < 0.001 Ya
IPM −1.289 0.001 Ya
ASI −2.579 < 0.001 Ya
Miskin 0.077 0.940 Tidak
Desa.Posbindu 0.042 0.024 Ya
Wasting 1.407 0.009 Ya
lag.IPM −2.269 0.007 Ya
lag.ASI −5.124 < 0.001 Ya
lag.Miskin 1.132 0.564 Tidak
lag.Desa.Posbindu 0.044 0.056 Tidak
lag.Wasting 0.426 0.472 Tidak
Lambda (λ) −1.451 < 0.001 Ya 
sarma_model <-sacsarlm(Stunting ~ IPM+ASI+Miskin + Desa.Posbindu+ Wasting , data   = ntt_merge,listw  = lw_knn)
summary(sarma_model, Nagelkerke = TRUE)
## 
## Call:sacsarlm(formula = Stunting ~ IPM + ASI + Miskin + Desa.Posbindu + 
##     Wasting, data = ntt_merge, listw = lw_knn)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -4.909492 -2.609758  0.018563  2.483887  4.943520 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value  Pr(>|z|)
## (Intercept)   129.293069  65.670599  1.9688 0.0489747
## IPM            -0.644585   0.273634 -2.3557 0.0184903
## ASI            -0.997532   0.535058 -1.8643 0.0622736
## Miskin          1.488989   0.497060  2.9956 0.0027391
## Desa.Posbindu   0.025314   0.013118  1.9296 0.0536512
## Wasting         0.687332   0.191123  3.5963 0.0003228
## 
## Rho: 0.66851
## Asymptotic standard error: 0.074374
##     z-value: 8.9886, p-value: < 2.22e-16
## Lambda: -1.3441
## Asymptotic standard error: 0.25237
##     z-value: -5.3258, p-value: 1.0049e-07
## 
## LR test value: 17.696, p-value: 0.00014364
## 
## Log likelihood: -63.20515 for sac model
## ML residual variance (sigma squared): 10.408, (sigma: 3.2261)
## Nagelkerke pseudo-R-squared: 0.76394 
## Number of observations: 22 
## Number of parameters estimated: 9 
## AIC: 144.41, (AIC for lm: 158.11)
Variabel Estimasi Std. Error z-value p-value Keterangan
Intercept 129.2931 65.6706 1.9688 0.0490 Signifikan (α=5%)
IPM -0.6446 0.2736 -2.3557 0.0185 Signifikan (α=5%)
ASI -0.9975 0.5351 -1.8643 0.0623 Signifikan (α=10%)
Miskin 1.4890 0.4971 2.9956 0.0027 Signifikan (α=5%)
Desa Posbindu 0.0253 0.0131 1.9296 0.0537 Signifikan (α=10%)
Wasting 0.6873 0.1911 3.5963 0.0003 Signifikan (α=5%)
Rho (ρ) 0.6685 0.0744 8.9886 <0.0001 Signifikan
Lambda (λ) -1.3441 0.2524 -5.3258 <0.0001 Signifikan

6.5.1 Statistik Model SARMA

Statistik Nilai
Log Likelihood -63.2052
Varians Residual (σ²) 10.408
Sigma (σ) 3.2261
Nagelkerke Pseudo-R² 0.7639
AIC SARMA 144.41
AIC OLS 158.11
LR test 17.696
p-value LR test 0.00014
Jumlah observasi 22
Jumlah parameter 9

7 Perbandingan Kinerja Model

Pemilihan model terbaik dapat dilakukan dengan kriteria AIC terkecil.

aic_table <- data.frame(
  Model  = c("OLS", "SAR", "SEM", "SDM", "SDEM","SARMA"),
  AIC    = c(
    AIC(ols_model),
    AIC(sar_model),
    AIC(sem_model),
    AIC(sdm_model),
    AIC(sdem_model),
    AIC(sarma_model)
  )
)
aic_table$AIC    <- round(aic_table$AIC, 4)
aic_table        <- aic_table[order(aic_table$AIC), ]
aic_table$Rank   <- 1:nrow(aic_table)

print(aic_table)
##   Model      AIC Rank
## 5  SDEM 136.8798    1
## 6 SARMA 144.4103    2
## 4   SDM 144.8758    3
## 2   SAR 154.3686    4
## 3   SEM 157.0850    5
## 1   OLS 158.1067    6

8 EKSTRAKSI ERROR SARMA: epsilon (ε̂) dan u (û)

8.1 Matriks bobot spasial (W) dalam bentuk matrix

W <- listw2mat(lw_knn)
n <- nrow(ntt_merge)
I <- diag(n)

# Parameter estimasi
rho    <- sarma_model$rho
lambda <- sarma_model$lambda

# ε̂: pure error 
# Dari u = λWu + ε → ε = (I - λW)u
u_hat     <- residuals(sarma_model)           # û = residual dari sacsarlm
eps_hat   <- (I - lambda * W) %*% u_hat      # ε̂ = (I - λW)û

# Tabel hasil per kabupaten/kota
error_table <- data.frame(
  Kabupaten_Kota = ntt_merge$WADMKK,
  u_hat          = round(as.numeric(u_hat), 4),
  eps_hat        = round(as.numeric(eps_hat), 4)
)

8.2 Peta SARMA

library(DT)
datatable(
  error_table,
  colnames = c("Kabupaten/Kota", "û (Spasial Error)", "ε̂ (Pure Error)"),
  caption  = "Nilai Error Model SARMA per Kabupaten/Kota",
  options  = list(pageLength = 22, dom = "t")
) %>%
  formatRound(columns = c("u_hat", "eps_hat"), digits = 4)
ntt_merge$u_hat   <- as.numeric(u_hat)
ntt_merge$eps_hat <- as.numeric(eps_hat)

8.2.1 Peta û

tm_shape(ntt_merge) +
  tm_polygons(
    col      = "u_hat",
    palette  = "RdBu",
    midpoint = 0,
    title    = "û (Spasial Error)"
  ) +
  tm_text("WADMKK", size = 0.4) +
  tm_layout(
    main.title     = "Peta û — Error Spasial Model SARMA",
    legend.outside = TRUE, frame = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'midpoint', 'palette' (rename to 'values') to
## fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdBu" is named
## "brewer.rd_bu"
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.

8.2.2 Peta ε̂

tm_shape(ntt_merge) +
  tm_polygons(
    col      = "eps_hat",
    palette  = "RdBu",
    midpoint = 0,
    title    = "ε̂ (Pure Error)"
  ) +
  tm_text("WADMKK", size = 0.4) +
  tm_layout(
    main.title     = "Peta ε̂ — Pure Error Model SARMA",
    legend.outside = TRUE, frame = FALSE
  )
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'midpoint', 'palette' (rename to 'values') to
## fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdBu" is named
## "brewer.rd_bu"
## Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.