LOAD LIBRARY

Memuat seluruh package R yang dibutuhkan untuk analisis

library(readr)
library(tmap)
library(sf)

## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

library(spdep)

## 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')`

library(spatialreg)

## 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(dplyr)

## 
## 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(ggplot2)
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

library(nortest)
library(DT)

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:/KULIAH/Semester 6/Stat Spasial/PROJEK/ntt.shp")

## Reading layer `ntt' from data source 
##   `D:\KULIAH\Semester 6\Stat 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:/Downloads/NTT.STUNTING.csv", show_col_types = FALSE)
data_ntt

## # A tibble: 23 × 18
##    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
## # ℹ 13 more rows
## # ℹ 10 more variables: Miskin <dbl>, Gizi.Buruk <dbl>, Wasting <dbl>,
## #   Desa.Posbindu <dbl>, Posyandu <dbl>, RLS <dbl>, Hunian.Layak <dbl>,
## #   Air.Minum.Layak <dbl>, Gini.Ratio <dbl>, Pengeluaran <dbl>

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, Gini.Ratio, Pengeluaran, geometry

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

AUTOKORELASI SPASIAL

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)

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

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).

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)

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. 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.

Model Spasial

1. Model OLS

Variabel respons: Y = presentase balita stunting (%) Variabel prediktor: X1 = Indeks Pembangunan Manusia X2 = persentase penduduk miskin (%) X3 = Rumah tangga dengan sanitasi layak (%) X4 = Pemberian ASI pada bayi umur 0-23 bulan (Baduta) (%)

ols_model <- lm(Stunting ~  IPM + Miskin + Sanitasi.Layak + ASI , data = ntt_merge)
summary(ols_model)

## 
## Call:
## lm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + ASI, 
##     data = ntt_merge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.963  -5.646  -1.247   4.045  16.300 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)    139.6464   137.0098   1.019    0.322
## IPM             -0.4042     0.7284  -0.555    0.586
## Miskin           1.2639     1.8066   0.700    0.494
## Sanitasi.Layak  -0.1254     0.1768  -0.709    0.488
## ASI             -0.9372     0.9317  -1.006    0.329
## 
## Residual standard error: 8.313 on 17 degrees of freedom
## Multiple R-squared:  0.312,  Adjusted R-squared:  0.1501 
## F-statistic: 1.927 on 4 and 17 DF,  p-value: 0.1521

AIC(ols_model)

## [1] 161.944

Uji Asumsi OLS

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.10543, p-value = 0.946
## 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.946 > 0.05 sehingga H₀ tidak ditolak. Dengan demikian Residual berdistribusi normal.

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 = 0.89434, df = 4, p-value = 0.9254

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

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.3092, p-value = 0.0311
## alternative hypothesis: true autocorrelation is greater than 0

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

Uji Multikolinearitas (VIF)

Nilai VIF > 10 mengindikasikan multikolinearitas.

vif_result <- vif(ols_model)
print(round(vif_result, 3))

##            IPM         Miskin Sanitasi.Layak            ASI 
##          2.685          2.888          1.407          1.658

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

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.

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 = 2.6673, p-value = 0.003823
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.35339801       -0.04761905        0.02260324

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

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 + Miskin + Sanitasi.Layak + ASI,
## data = ntt_merge)
## test weights: lw_knn
## 
## RSerr = 4.6104, df = 1, p-value = 0.03178
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + ASI,
## data = ntt_merge)
## test weights: lw_knn
## 
## RSlag = 7.6387, df = 1, p-value = 0.005713
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + ASI,
## data = ntt_merge)
## test weights: lw_knn
## 
## adjRSerr = 1.4088, df = 1, p-value = 0.2352
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + ASI,
## data = ntt_merge)
## test weights: lw_knn
## 
## adjRSlag = 4.4371, df = 1, p-value = 0.03517
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + ASI,
## data = ntt_merge)
## test weights: lw_knn
## 
## SARMA = 9.0475, df = 2, p-value = 0.01085

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)	4.610	0.032	Ya
LM-Lag (RSlag)	7.639	0.006	Ya
Robust LM-Error	1.409	0.235	Tidak
Robust LM-Lag	4.437	0.035	Ya
SARMA	9.048	0.011	Ya

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 + Miskin + Sanitasi.Layak + ASI, data = ntt_merge,listw = lw_knn)

summary(sar_model, Nagelkerke = TRUE)

## 
## Call:lagsarlm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + 
##     ASI, data = ntt_merge, listw = lw_knn)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -9.07806 -5.04435 -0.95296  2.25506 15.78062 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)    118.194691 105.084094  1.1248   0.2607
## IPM             -0.404009   0.562059 -0.7188   0.4723
## Miskin           1.265153   1.386721  0.9123   0.3616
## Sanitasi.Layak  -0.050416   0.141404 -0.3565   0.7214
## ASI             -0.862518   0.716000 -1.2046   0.2283
## 
## Rho: 0.40065, LR test value: 4.8932, p-value: 0.026962
## Asymptotic standard error: 0.19124
##     z-value: 2.095, p-value: 0.03617
## Wald statistic: 4.3891, p-value: 0.03617
## 
## Log likelihood: -72.5254 for lag model
## ML residual variance (sigma squared): 40.627, (sigma: 6.374)
## Nagelkerke pseudo-R-squared: 0.44919 
## Number of observations: 22 
## Number of parameters estimated: 7 
## AIC: 159.05, (AIC for lm: 161.94)
## LM test for residual autocorrelation
## test value: 4.6049, p-value: 0.031881

Rho (ρ) = 0.401 ≠ 0 dengan p-value 0.031881 < 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.

4. SEM

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

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

summary(sem_model,Nagelkerke = TRUE)

## 
## Call:errorsarlm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + 
##     ASI, data = ntt_merge, listw = lw_knn)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -9.23531 -4.47257 -0.47335  2.69611 16.37723 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)    118.36101  109.71034  1.0789   0.2807
## IPM             -0.24156    0.59406 -0.4066   0.6843
## Miskin           1.04591    1.47527  0.7090   0.4783
## Sanitasi.Layak  -0.13073    0.16080 -0.8130   0.4162
## ASI             -0.81802    0.71875 -1.1381   0.2551
## 
## Lambda: 0.37771, LR test value: 3.3194, p-value: 0.068468
## Asymptotic standard error: 0.20867
##     z-value: 1.8101, p-value: 0.070282
## Wald statistic: 3.2764, p-value: 0.070282
## 
## Log likelihood: -73.31233 for error model
## ML residual variance (sigma squared): 43.916, (sigma: 6.6269)
## Nagelkerke pseudo-R-squared: 0.40834 
## Number of observations: 22 
## Number of parameters estimated: 7 
## AIC: 160.62, (AIC for lm: 161.94)

parameter λ = 0.38 dengan p-value 0.0702 > 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.

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 + Miskin + Sanitasi.Layak + ASI, data = ntt_merge,listw = lw_knn, type = "mixed")

summary(sdm_model,Nagelkerke = TRUE)

## 
## Call:lagsarlm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + 
##     ASI, data = ntt_merge, listw = lw_knn, type = "mixed")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -8.26062 -3.93772  0.78134  3.65922  9.92190 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)        491.762383 343.425894  1.4319  0.15216
## IPM                 -1.060954   0.520381 -2.0388  0.04147
## Miskin               1.777826   1.242588  1.4307  0.15250
## Sanitasi.Layak       0.091392   0.151013  0.6052  0.54505
## ASI                 -1.927740   0.856153 -2.2516  0.02435
## lag.IPM             -2.047199   1.410836 -1.4511  0.14676
## lag.Miskin           4.556906   2.924068  1.5584  0.11914
## lag.Sanitasi.Layak   0.503696   0.215114  2.3415  0.01920
## lag.ASI             -1.683834   2.238156 -0.7523  0.45185
## 
## Rho: 0.07436, LR test value: 0.14061, p-value: 0.70767
## Asymptotic standard error: 0.21079
##     z-value: 0.35277, p-value: 0.72426
## Wald statistic: 0.12444, p-value: 0.72426
## 
## Log likelihood: -67.24717 for mixed model
## ML residual variance (sigma squared): 26.417, (sigma: 5.1398)
## Nagelkerke pseudo-R-squared: 0.65911 
## Number of observations: 22 
## Number of parameters estimated: 11 
## AIC: 156.49, (AIC for lm: 154.63)
## LM test for residual autocorrelation
## test value: 1.5219, p-value: 0.21734

Variabel	Estimasi	p-value	Signifikan?
(Intercept)	491.762	0.152	Tidak
IPM	−1.061	0.041	Ya
Miskin	1.778	0.153	Tidak
Sanitasi.Layak	0.091	0.545	Tidak
ASI	−1.928	0.024	Ya
lag.IPM	−2.047	0.147	Tidak
lag.Miskin	4.557	0.119	Tidak
lag.Sanitasi.Layak	0.504	0.019	Ya
lag.ASI	−1.684	0.452	Tidak
Rho (ρ)	0.074	0.724	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 + Miskin + Sanitasi.Layak + ASI,data   = ntt_merge,listw  = lw_knn,etype  = "emixed"
)
summary(sdem_model,Nagelkerke = TRUE)

## 
## Call:errorsarlm(formula = Stunting ~ IPM + Miskin + Sanitasi.Layak + 
##     ASI, data = ntt_merge, listw = lw_knn, etype = "emixed")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.0418 -2.6809  0.8716  2.4973  5.9122 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                       Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)        1477.638512  199.975220  7.3891 1.479e-13
## IPM                  -1.277939    0.438910 -2.9116 0.0035956
## Miskin                1.421327    1.111607  1.2786 0.2010298
## Sanitasi.Layak       -0.096042    0.136441 -0.7039 0.4814885
## ASI                  -3.382138    0.655150 -5.1624 2.438e-07
## lag.IPM              -4.196187    1.087262 -3.8594 0.0001137
## lag.Miskin            3.354564    2.656191  1.2629 0.2066170
## lag.Sanitasi.Layak    0.990479    0.161483  6.1336 8.589e-10
## lag.ASI              -8.783675    1.599069 -5.4930 3.952e-08
## 
## Lambda: -1.4608, LR test value: 3.5551, p-value: 0.059362
## Asymptotic standard error: 0.20727
##     z-value: -7.048, p-value: 1.8145e-12
## Wald statistic: 49.675, p-value: 1.8145e-12
## 
## Log likelihood: -65.53993 for error model
## ML residual variance (sigma squared): 14.331, (sigma: 3.7856)
## Nagelkerke pseudo-R-squared: 0.70812 
## Number of observations: 22 
## Number of parameters estimated: 11 
## AIC: 153.08, (AIC for lm: 154.63)

Variabel	Estimasi	p-value	Signifikan
(Intercept)	1477.639	< 0.001	Ya
IPM	−1.278	0.004	Ya
Miskin	1.421	0.201	Tidak
Sanitasi.Layak	−0.096	0.481	Tidak
ASI	−3.382	< 0.001	Ya
lag.IPM	−4.196	< 0.001	Ya
lag.Miskin	3.355	0.207	Tidak
lag.Sanitasi.Layak	0.990	< 0.001	Ya
lag.ASI	−8.784	< 0.001	Ya
Lambda (λ)	−1.461	< 0.001	Ya

Model memenuhi karakteristik Spatial Durbin Error Model (SDEM), yaitu terdapat pengaruh langsung dari variabel independen (IPM dan ASI) serta pengaruh tidak langsung atau efek spasial (spillover effect) dari wilayah sekitar (lag.IPM, lag.Sanitasi.Layak, dan lag.ASI). Selain itu, parameter Lambda (λ) yang signifikan menunjukkan adanya ketergantungan spasial pada komponen error.

Nilai Nagelkerke pseudo-Rsquare sebesar 0,70812 yang berarti 70,81% variasi prevalensi stunting di Provinsi NTT dapat dijelaskan oleh variabel IPM, kemiskinan, sanitasi layak, ASI, serta efek spasial antarwilayah.

Perbandingan Kinerja Model

Pemilihan model terbaik dapat dilakukan dengan kriteria AIC terkecil.

aic_table <- data.frame(
  Model  = c("OLS", "SAR", "SEM", "SDM", "SDEM"),
  AIC    = c(
    AIC(ols_model),
    AIC(sar_model),
    AIC(sem_model),
    AIC(sdm_model),
    AIC(sdem_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 153.0799    1
## 4   SDM 156.4943    2
## 2   SAR 159.0508    3
## 3   SEM 160.6247    4
## 1   OLS 161.9440    5

Model Spatial Durbin Error Model (SDEM) memiliki nilai AIC terkecil yaitu 153,0799, sehingga dipilih sebagai model terbaik dalam menjelaskan prevalensi stunting di Provinsi NTT dibandingkan model OLS, SAR, SEM, dan SDM.

Analisis Spasial Prevalensi Stunting di NTT

2026-05-14