library(terra)
## Warning: package 'terra' was built under R version 4.5.1
## terra 1.8.70
library (dplyr)
## Warning: package 'dplyr' was built under R version 4.5.1
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
##
## intersect, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library (readxl)
library (sf)
## Warning: package 'sf' was built under R version 4.5.1
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library (ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.1
library (stringr)
library (plotly)
## Warning: package 'plotly' was built under R version 4.5.1
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library (spatstat.geom)
## Warning: package 'spatstat.geom' was built under R version 4.5.1
## Loading required package: spatstat.data
## Warning: package 'spatstat.data' was built under R version 4.5.1
## Loading required package: spatstat.univar
## Warning: package 'spatstat.univar' was built under R version 4.5.1
## spatstat.univar 3.1-4
## spatstat.geom 3.5-0
##
## Attaching package: 'spatstat.geom'
## The following objects are masked from 'package:terra':
##
## area, delaunay, is.empty, rescale, rotate, shift, where.max,
## where.min
library (patchwork)
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:spatstat.geom':
##
## area
## The following object is masked from 'package:terra':
##
## area
library (spatstat)
## Warning: package 'spatstat' was built under R version 4.5.1
## Loading required package: spatstat.random
## Warning: package 'spatstat.random' was built under R version 4.5.1
## spatstat.random 3.4-1
## Loading required package: spatstat.explore
## Warning: package 'spatstat.explore' was built under R version 4.5.1
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## spatstat.explore 3.5-2
## Loading required package: spatstat.model
## Warning: package 'spatstat.model' was built under R version 4.5.1
## Loading required package: rpart
## spatstat.model 3.4-0
## Loading required package: spatstat.linnet
## Warning: package 'spatstat.linnet' was built under R version 4.5.1
## spatstat.linnet 3.3-1
##
## spatstat 3.4-0
## For an introduction to spatstat, type 'beginner'
library (purrr)
library (psych)
##
## Attaching package: 'psych'
## The following objects are masked from 'package:spatstat.geom':
##
## reflect, rescale
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
## The following objects are masked from 'package:terra':
##
## describe, distance, rescale
library (spdep)
## Warning: package 'spdep' was built under R version 4.5.1
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.5.1
## 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)
## Warning: package 'spatialreg' was built under R version 4.5.1
## 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 (raster)
## Warning: package 'raster' was built under R version 4.5.1
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.5.1
##
## Attaching package: 'raster'
## The following object is masked from 'package:nlme':
##
## getData
## The following object is masked from 'package:plotly':
##
## select
## The following object is masked from 'package:dplyr':
##
## select
library (lmtest)
## Warning: package 'lmtest' was built under R version 4.5.1
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following object is masked from 'package:terra':
##
## time<-
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(tmap)
## Warning: package 'tmap' was built under R version 4.5.1
data<-read_xlsx("C:/Users/Fathoni Sabri/OneDrive/Documents/Semester 5/RegSpas/Data/Data Baru(Jawa Barat).xlsx")
data
## # A tibble: 5,311 × 13
## `KODE PROV` `NAMA PROVINSI` `KODE KAB` `NAMA KABUPATEN` `KODE KEC`
## <dbl> <chr> <dbl> <chr> <dbl>
## 1 32 JAWA BARAT 3201 BOGOR 320102
## 2 32 JAWA BARAT 3201 BOGOR 320102
## 3 32 JAWA BARAT 3201 BOGOR 320102
## 4 32 JAWA BARAT 3201 BOGOR 320102
## 5 32 JAWA BARAT 3201 BOGOR 320102
## 6 32 JAWA BARAT 3201 BOGOR 320102
## 7 32 JAWA BARAT 3201 BOGOR 320102
## 8 32 JAWA BARAT 3201 BOGOR 320102
## 9 32 JAWA BARAT 3201 BOGOR 320102
## 10 32 JAWA BARAT 3201 BOGOR 320102
## # ℹ 5,301 more rows
## # ℹ 8 more variables: `NAMA KECAMATAN` <chr>, `KODE DESA` <dbl>,
## # `NAMA DESA` <chr>, `Jarak sarkes terdekat` <dbl>, `Jumlah Bidan` <dbl>,
## # `Jumlah dokter` <dbl>, `Jumlah posyandu` <dbl>,
## # `Jumlah bayi stunting` <dbl>
# Agregasi berdasarkan nama kabupaten
data_jabar_kab <- data %>%
group_by(`NAMA KABUPATEN`) %>%
summarise(
total_stunting = sum(`Jumlah bayi stunting`, na.rm = TRUE)
) %>%
arrange(`NAMA KABUPATEN`) # urutkan nama kabupaten
# Lihat hasil agregasi
data_jabar_kab
## # A tibble: 19 × 2
## `NAMA KABUPATEN` total_stunting
## <chr> <dbl>
## 1 BANDUNG 3132
## 2 BANDUNG BARAT 1765
## 3 BEKASI 429
## 4 BOGOR 1521
## 5 CIAMIS 924
## 6 CIANJUR 1459
## 7 CIREBON 2566
## 8 GARUT 5498
## 9 INDRAMAYU 437
## 10 KARAWANG 436
## 11 KOTA BANJAR 104
## 12 KUNINGAN 1346
## 13 MAJALENGKA 529
## 14 PANGANDARAN 124
## 15 PURWAKARTA 400
## 16 SUBANG 250
## 17 SUKABUMI 6316
## 18 SUMEDANG 4291
## 19 TASIKMALAYA 5365
desa <- st_read("C:/Users/Fathoni Sabri/OneDrive/Documents/Semester 5/RegSpas/idn_admbnda_adm4_bps_20200401.shp")
## Reading layer `idn_admbnda_adm4_bps_20200401' from data source
## `C:\Users\Fathoni Sabri\OneDrive\Documents\Semester 5\RegSpas\idn_admbnda_adm4_bps_20200401.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 81912 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS: WGS 84
kecamatan <- st_read("C:/Users/Fathoni Sabri/OneDrive/Documents/Semester 5/RegSpas/idn_admbnda_adm3_bps_20200401.shp")
## Reading layer `idn_admbnda_adm3_bps_20200401' from data source
## `C:\Users\Fathoni Sabri\OneDrive\Documents\Semester 5\RegSpas\idn_admbnda_adm3_bps_20200401.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7069 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS: WGS 84
kabupaten <- st_read("C:/Users/Fathoni Sabri/OneDrive/Documents/Semester 5/RegSpas/idn_admbnda_adm2_bps_20200401.shp")
## Reading layer `idn_admbnda_adm2_bps_20200401' from data source
## `C:\Users\Fathoni Sabri\OneDrive\Documents\Semester 5\RegSpas\idn_admbnda_adm2_bps_20200401.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 522 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS: WGS 84
provinsi <- st_read("C:/Users/Fathoni Sabri/OneDrive/Documents/Semester 5/RegSpas/idn_admbnda_adm1_bps_20200401.shp")
## Reading layer `idn_admbnda_adm1_bps_20200401' from data source
## `C:\Users\Fathoni Sabri\OneDrive\Documents\Semester 5\RegSpas\idn_admbnda_adm1_bps_20200401.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS: WGS 84
shp_prov <- provinsi %>%
mutate(ADM1_PCODE_clean = sub("^ID", "", ADM1_PCODE))
shp_kab <- kabupaten %>%
mutate(ADM2_PCODE_clean = sub("^ID", "", ADM2_PCODE))
shp_jabar <- shp_kab %>%
filter(ADM1_EN == "Jawa Barat")
shp_desa <- desa %>%
mutate(ADM4_PCODE_clean= sub("^ID", "", ADM4_PCODE))
shp_desa <- shp_desa %>%
filter(ADM1_EN == "Jawa Barat")
shp_desa_clean <- shp_desa %>%
group_by(ADM4_EN) %>%
summarise(geometry = st_union(geometry),
.groups = "drop")
# 1. Pastikan kolom kode kabupaten bertipe karakter
data_jabar_kab$`NAMA KABUPATEN` <- as.character(data_jabar_kab$`NAMA KABUPATEN`)
# 2. Buat data agregasi jumlah bidan per kabupaten
stunting_kab <- data_jabar_kab %>%
group_by(`NAMA KABUPATEN`) %>%
summarise(total_stunting = sum(total_stunting, na.rm = TRUE))
# 3. Siapkan shapefile kabupaten Jawa Barat
kab_jabar <- kabupaten %>%
filter(ADM1_EN == "Jawa Barat") %>%
mutate(NAMA_KABUPATEN = toupper(ADM2_EN))
# 4. Gabungkan shapefile dengan data total bidan
shp_gab_stunting <- kab_jabar %>%
left_join(stunting_kab, by = c("NAMA_KABUPATEN" = "NAMA KABUPATEN"))
# 5. Plot peta total bidan
ggplot(shp_gab_stunting) +
geom_sf(aes(fill = total_stunting), color = "grey40", size = 0.2) +
scale_fill_gradient(
low = "#81C784", # hijau muda (baik)
high = "#E53935", # merah tua (buruk)
na.value = "white"
) +
theme_minimal() +
labs(
title = "Total Jumlah Bayi Stunting per Kabupaten di Jawa Barat",
fill = "Total Stunting"
)
Peta choropleth yang menampilkan Total Jumlah Bayi Stunting per Kabupaten di Jawa Barat menunjukkan adanya variasi dan pola spasial yang signifikan dalam masalah ini. Peta ini menggunakan gradasi warna (dari hijau muda ke merah tua) untuk memvisualisasikan total kasus stunting, di mana warna merah tua mengindikasikan total kasus stunting yang sangat tinggi (lebih dari 6.000) dan warna hijau muda mengindikasikan total kasus yang rendah (kurang dari 2.000).
# Daftar nama kabupaten yang ingin dipilih
kabupaten_keep <- c(
"BANDUNG", "BANDUNG BARAT", "BEKASI", "BOGOR",
"CIAMIS", "CIANJUR", "CIREBON", "GARUT",
"INDRAMAYU", "KARAWANG","KOTA BANJAR", "KUNINGAN", "MAJALENGKA",
"PANGANDARAN", "PURWAKARTA", "SUBANG", "SUKABUMI",
"SUMEDANG", "TASIKMALAYA"
)
# Filter shapefile agar hanya kabupaten yang dipilih
shp_kab_stunting <- shp_gab_stunting %>%
filter(NAMA_KABUPATEN %in% kabupaten_keep)
# Lihat hasil nama kabupaten
shp_kab_stunting %>%
st_drop_geometry() %>%
dplyr::select(NAMA_KABUPATEN)
## NAMA_KABUPATEN
## 1 BANDUNG
## 2 BANDUNG BARAT
## 3 BEKASI
## 4 BOGOR
## 5 CIAMIS
## 6 CIANJUR
## 7 CIREBON
## 8 GARUT
## 9 INDRAMAYU
## 10 KARAWANG
## 11 KOTA BANJAR
## 12 KUNINGAN
## 13 MAJALENGKA
## 14 PANGANDARAN
## 15 PURWAKARTA
## 16 SUBANG
## 17 SUKABUMI
## 18 SUMEDANG
## 19 TASIKMALAYA
# Ambil centroid setiap kabupaten
coords <- st_coordinates(st_centroid(st_geometry(shp_kab_stunting)))
# Hitung matriks jarak Euclidean
D <- as.matrix(dist(coords, method = "euclidean"))
# Tentukan ambang batas jarak (misalnya 1/3 dari maksimum)
threshold <- max(D) / 3
# Buat matriks pembobot
n <- nrow(D)
W.eu <- matrix(0, n, n)
for (i in 1:n) {
for (j in 1:n) {
if (D[i, j] <= threshold && i != j) {
W.eu[i, j] <- 1
}
}
}
# Tambahkan nama kabupaten sebagai label
rownames(W.eu) <- shp_kab_stunting$NAMA_KABUPATEN
colnames(W.eu) <- shp_kab_stunting$NAMA_KABUPATEN
W.eu
## BANDUNG BANDUNG BARAT BEKASI BOGOR CIAMIS CIANJUR CIREBON GARUT
## BANDUNG 0 1 0 0 0 1 0 1
## BANDUNG BARAT 1 0 0 0 0 1 0 1
## BEKASI 0 0 0 1 0 0 0 0
## BOGOR 0 0 1 0 0 0 0 0
## CIAMIS 0 0 0 0 0 0 1 1
## CIANJUR 1 1 0 0 0 0 0 1
## CIREBON 0 0 0 0 1 0 0 0
## GARUT 1 1 0 0 1 1 0 0
## INDRAMAYU 0 0 0 0 0 0 1 0
## KARAWANG 0 1 1 1 0 0 0 0
## KOTA BANJAR 0 0 0 0 1 0 1 0
## KUNINGAN 0 0 0 0 1 0 1 0
## MAJALENGKA 0 0 0 0 1 0 1 0
## PANGANDARAN 0 0 0 0 1 0 0 0
## PURWAKARTA 1 1 1 1 0 1 0 0
## SUBANG 1 1 1 0 0 0 0 0
## SUKABUMI 0 0 0 1 0 1 0 0
## SUMEDANG 1 1 0 0 1 0 1 1
## TASIKMALAYA 1 0 0 0 1 0 0 1
## INDRAMAYU KARAWANG KOTA BANJAR KUNINGAN MAJALENGKA PANGANDARAN
## BANDUNG 0 0 0 0 0 0
## BANDUNG BARAT 0 1 0 0 0 0
## BEKASI 0 1 0 0 0 0
## BOGOR 0 1 0 0 0 0
## CIAMIS 0 0 1 1 1 1
## CIANJUR 0 0 0 0 0 0
## CIREBON 1 0 1 1 1 0
## GARUT 0 0 0 0 0 0
## INDRAMAYU 0 0 0 1 1 0
## KARAWANG 0 0 0 0 0 0
## KOTA BANJAR 0 0 0 1 1 1
## KUNINGAN 1 0 1 0 1 1
## MAJALENGKA 1 0 1 1 0 0
## PANGANDARAN 0 0 1 1 0 0
## PURWAKARTA 0 1 0 0 0 0
## SUBANG 1 1 0 0 1 0
## SUKABUMI 0 0 0 0 0 0
## SUMEDANG 1 0 0 1 1 0
## TASIKMALAYA 0 0 1 1 0 1
## PURWAKARTA SUBANG SUKABUMI SUMEDANG TASIKMALAYA
## BANDUNG 1 1 0 1 1
## BANDUNG BARAT 1 1 0 1 0
## BEKASI 1 1 0 0 0
## BOGOR 1 0 1 0 0
## CIAMIS 0 0 0 1 1
## CIANJUR 1 0 1 0 0
## CIREBON 0 0 0 1 0
## GARUT 0 0 0 1 1
## INDRAMAYU 0 1 0 1 0
## KARAWANG 1 1 0 0 0
## KOTA BANJAR 0 0 0 0 1
## KUNINGAN 0 0 0 1 1
## MAJALENGKA 0 1 0 1 0
## PANGANDARAN 0 0 0 0 1
## PURWAKARTA 0 1 0 1 0
## SUBANG 1 0 0 1 0
## SUKABUMI 0 0 0 0 0
## SUMEDANG 1 1 0 0 0
## TASIKMALAYA 0 0 0 0 0
Matriks biner ini menunjukkan struktur ketetanggaan spasial berdasarkan jarak: \(Nilai 1 → ada hubungan spasial (wilayah berdekatan).\) \(Nilai 0 → tidak ada hubungan spasial.\)
#inverse weight matrix
alpha1=1
w.idw=1/(D^alpha1)
# inverse weight matrix - row-normalized
diag(w.idw)<-0
rtot<-rowSums(w.idw, na.rm =T)
w.idw.sd<-w.idw/rtot
rowSums(w.idw.sd, na.rm=T)
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
w.idw.s <- mat2listw(w.idw.sd,style='W') #untuk melihat matriks W
summary(w.idw.s)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 19
## Number of nonzero links: 342
## Percentage nonzero weights: 94.73684
## Average number of links: 18
## Link number distribution:
##
## 18
## 19
## 19 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 with 18 links
## 19 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 with 18 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 19 361 19 2.832473 76.35231
#invers distance weight - plot
plot(st_geometry(shp_kab_stunting), border = "green", col = "gray90",
main = "Inverse Distance Weight (alpha = 1)")
plot(w.idw.s, coords, add = TRUE, col = "blue")
alpha2=2
w.idw2<-1/(D^alpha2)
diag(w.idw2)<-0
rtot<-rowSums(w.idw2,na.rm=TRUE)
w.idw2.sd<-w.idw2/rtot #row-normalized
w.idw2.s = mat2listw(w.idw2.sd,style='W')
w.idw2.s <- mat2listw(w.idw2.sd,style='W')
summary(w.idw2.s)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 19
## Number of nonzero links: 342
## Percentage nonzero weights: 94.73684
## Average number of links: 18
## Link number distribution:
##
## 18
## 19
## 19 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 with 18 links
## 19 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 with 18 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 19 361 19 5.565555 77.19852
#invers distance weight - plot
plot(st_geometry(shp_kab_stunting), border = "green", col = "gray90",
main = "Inverse Distance Weight (alpha = 2)")
plot(w.idw2.s, coords, add = TRUE, col = "blue")
Dari kedua peta Inverse Distance Weight tersebut, terlihat bahwa ketika nilai \(alpha = 1\), hubungan spasial antar wilayah (garis biru) tampak lebih merata karena bobot kedekatan antar kabupaten menurun secara perlahan terhadap jarak. Artinya, wilayah yang agak jauh masih memiliki pengaruh spasial meskipun kecil. Namun, ketika \(alpha = 2\), pengaruh jarak menjadi jauh lebih kuat — bobot antar wilayah menurun sangat cepat sehingga hanya wilayah yang benar-benar berdekatan yang memiliki pengaruh berarti. Dengan kata lain, semakin besar nilai alpha, semakin lokal pengaruh spasial yang ditangkap oleh model. Jadi, pemilihan alpha menentukan sejauh mana suatu wilayah “merasakan” efek dari wilayah lain di sekitarnya: alpha kecil menggambarkan keterkaitan yang luas, sedangkan alpha besar menggambarkan keterkaitan yang lebih terbatas di sekitar tetangga terdekat.
# Ambil koordinat centroid kabupaten (jika belum ada)
coords <- st_coordinates(st_centroid(st_geometry(shp_kab_stunting)))
# 1️⃣ K = 4
w.knn.stunt <- knn2nb(knearneigh(coords, k = 4, longlat = TRUE))
w.knn.stunt.s <- nb2listw(w.knn.stunt, style = 'W')
summary(w.knn.stunt.s)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 19
## Number of nonzero links: 76
## Percentage nonzero weights: 21.05263
## Average number of links: 4
## Non-symmetric neighbours list
## Link number distribution:
##
## 4
## 19
## 19 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 with 4 links
## 19 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 with 4 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 19 361 19 8.75 78.25
# Plot hasil KNN (k=4)
plot(st_geometry(shp_kab_stunting), border = "darkgreen", col = "gray90",
main = "KNN Distance Weight (K = 4) - Data Stunting")
plot(w.knn.stunt.s, coords, add = TRUE, col = "red")
# 2️⃣ K = 3
w.knn2.stunt <- knn2nb(knearneigh(coords, k = 3, longlat = TRUE))
w.knn2.stunt.s <- nb2listw(w.knn2.stunt, style = 'W')
summary(w.knn2.stunt.s)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 19
## Number of nonzero links: 57
## Percentage nonzero weights: 15.78947
## Average number of links: 3
## Non-symmetric neighbours list
## Link number distribution:
##
## 3
## 19
## 19 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 with 3 links
## 19 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 with 3 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 19 361 19 11.66667 77.55556
# Plot hasil KNN (k=3)
plot(st_geometry(shp_kab_stunting), border = "darkgreen", col = "gray90",
main = "KNN Distance Weight (K = 3) - Data Stunting")
plot(w.knn2.stunt.s, coords, add = TRUE, col = "red")
Berdasarkan hasil visualisasi dan ringkasan bobot spasial, perbedaan antara K = 3 dan K = 4 pada metode K-Nearest Neighbors (KNN) menunjukkan variasi tingkat keterhubungan antar wilayah dalam data stunting. Ketika nilai K meningkat dari 3 menjadi 4, setiap wilayah memiliki satu tetangga tambahan, sehingga jumlah koneksi meningkat dari 57 menjadi 76 (rata-rata dari 3 menjadi 4 koneksi per wilayah). Hal ini tercermin pada peta, di mana pola garis merah tampak lebih rapat pada K = 4. Peningkatan jumlah tetangga menyebabkan struktur spasial menjadi lebih padat dan hubungan antarwilayah semakin kuat, yang dapat meningkatkan stabilitas model spasial tetapi berpotensi mengurangi kejelasan pengaruh lokal. Nilai konstanta S0, S1, dan S2 juga menunjukkan adanya perubahan distribusi bobot, di mana S1 menurun dan S2 meningkat pada K = 4, menandakan hubungan spasial yang lebih seimbang namun kurang tajam dalam membedakan wilayah yang sangat berdekatan.
# Membuat matriks keterhubungan spasial dengan metode ROOK untuk data stunting
rook_stunting <- poly2nb(as(shp_kab_stunting, "Spatial"), queen = FALSE)
# Melihat ringkasan struktur keterhubungan
summary(rook_stunting)
## Neighbour list object:
## Number of regions: 19
## Number of nonzero links: 78
## Percentage nonzero weights: 21.60665
## Average number of links: 4.105263
## Link number distribution:
##
## 1 2 3 4 5 6
## 1 3 2 4 5 4
## 1 least connected region:
## 11 with 1 link
## 4 most connected regions:
## 6 13 16 18 with 6 links
# Mengubah daftar tetangga menjadi matriks pembobot spasial
w.rook_stunting <- nb2listw(rook_stunting, style = 'W')
# Melihat sebagian isi matriks bobot
head(w.rook_stunting, 3)
## $style
## [1] "W"
##
## $neighbours
## Neighbour list object:
## Number of regions: 19
## Number of nonzero links: 78
## Percentage nonzero weights: 21.60665
## Average number of links: 4.105263
##
## $weights
## $weights[[1]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## $weights[[2]]
## [1] 0.25 0.25 0.25 0.25
##
## $weights[[3]]
## [1] 0.5 0.5
##
## $weights[[4]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## $weights[[5]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## $weights[[6]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## $weights[[7]]
## [1] 0.3333333 0.3333333 0.3333333
##
## $weights[[8]]
## [1] 0.25 0.25 0.25 0.25
##
## $weights[[9]]
## [1] 0.25 0.25 0.25 0.25
##
## $weights[[10]]
## [1] 0.25 0.25 0.25 0.25
##
## $weights[[11]]
## [1] 1
##
## $weights[[12]]
## [1] 0.3333333 0.3333333 0.3333333
##
## $weights[[13]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## $weights[[14]]
## [1] 0.5 0.5
##
## $weights[[15]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## $weights[[16]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## $weights[[17]]
## [1] 0.5 0.5
##
## $weights[[18]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## $weights[[19]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## attr(,"mode")
## [1] "binary"
## attr(,"W")
## [1] TRUE
## attr(,"comp")
## attr(,"comp")$d
## [1] 5 4 2 5 5 6 3 4 4 4 1 3 6 2 5 6 2 6 5
# Plot hubungan spasial ROOK untuk data stunting
plot(st_geometry(shp_kab_stunting), border = "blue", col = "gray",
main = "Rook Contiguity Weight - Stunting")
plot(rook_stunting, coords, add = TRUE, col = "red")
Hasil dari pembentukan matriks ketetanggaan menggunakan metode Rook Contiguity menunjukkan bahwa masalah stunting di Jawa Barat perlu didekati secara spasial karena adanya keterhubungan geografis antar-kabupaten. Secara rata-rata, setiap kabupaten berinteraksi langsung (berbagi batas darat) dengan sekitar \(4\) tetangga (\(4.105\) links), namun terdapat variasi besar; ada kabupaten yang relatif terisolasi (hanya 1 tetangga), dan ada pula yang sangat terhubung (hingga 6 tetangga). Karena kita menggunakan normalisasi baris (Style ‘W’), bobot pada matriks memastikan bahwa pengaruh total dari semua tetangga selalu bernilai \(1\). Ini berarti setiap kabupaten yang bertetangga akan membagi rata pengaruhnya, di mana bobot \(w_{ij}\) dihitung sebagai \(w_{ij} = \frac{1}{N_i}\) (\(N_i\) adalah jumlah tetangga kabupaten \(i\)). Struktur \(\mathbf{W}\) yang dihasilkan ini adalah fondasi utama untuk menguji dan memodelkan apakah tingkat stunting pada suatu kabupaten dipengaruhi secara signifikan oleh tingkat stunting (atau faktor penentu stunting) di kabupaten tetangganya.
# Membuat matriks keterhubungan spasial dengan metode QUEEN untuk data stunting
queen_stunting <- poly2nb(as(shp_kab_stunting, "Spatial"), queen = TRUE)
# Menampilkan ringkasan struktur keterhubungan
summary(queen_stunting)
## Neighbour list object:
## Number of regions: 19
## Number of nonzero links: 78
## Percentage nonzero weights: 21.60665
## Average number of links: 4.105263
## Link number distribution:
##
## 1 2 3 4 5 6
## 1 3 2 4 5 4
## 1 least connected region:
## 11 with 1 link
## 4 most connected regions:
## 6 13 16 18 with 6 links
# Mengubah daftar tetangga menjadi matriks pembobot spasial
w.queen_stunting <- nb2listw(queen_stunting, style = "W")
# Melihat sebagian isi matriks bobot
head(w.queen_stunting, 3)
## $style
## [1] "W"
##
## $neighbours
## Neighbour list object:
## Number of regions: 19
## Number of nonzero links: 78
## Percentage nonzero weights: 21.60665
## Average number of links: 4.105263
##
## $weights
## $weights[[1]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## $weights[[2]]
## [1] 0.25 0.25 0.25 0.25
##
## $weights[[3]]
## [1] 0.5 0.5
##
## $weights[[4]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## $weights[[5]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## $weights[[6]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## $weights[[7]]
## [1] 0.3333333 0.3333333 0.3333333
##
## $weights[[8]]
## [1] 0.25 0.25 0.25 0.25
##
## $weights[[9]]
## [1] 0.25 0.25 0.25 0.25
##
## $weights[[10]]
## [1] 0.25 0.25 0.25 0.25
##
## $weights[[11]]
## [1] 1
##
## $weights[[12]]
## [1] 0.3333333 0.3333333 0.3333333
##
## $weights[[13]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## $weights[[14]]
## [1] 0.5 0.5
##
## $weights[[15]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## $weights[[16]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## $weights[[17]]
## [1] 0.5 0.5
##
## $weights[[18]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## $weights[[19]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## attr(,"mode")
## [1] "binary"
## attr(,"W")
## [1] TRUE
## attr(,"comp")
## attr(,"comp")$d
## [1] 5 4 2 5 5 6 3 4 4 4 1 3 6 2 5 6 2 6 5
# Plot hubungan spasial QUEEN untuk data stunting
plot(st_geometry(shp_kab_stunting), border = "blue", col = "gray",
main = "Queen Contiguity Weight - Data Stunting")
plot(queen_stunting, coords, add = TRUE, col = "red")
Hasil pembentukan matriks bobot spasial menggunakan metode Queen Contiguity untuk data stunting di Jawa Barat menunjukkan bahwa, secara rata-rata, setiap kabupaten memiliki sekitar 4 tetangga (\(4.105\) links) yang berinteraksi dengannya, baik dengan berbagi batas maupun hanya berbagi titik sudut. Menariknya, hasil statistiknya (Number of nonzero links dan Average number of links) ternyata identik dengan metode Rook, mengindikasikan bahwa kedua definisi ketetanggaan menghasilkan struktur keterhubungan yang sama untuk dataset ini. Karena matriks telah dinormalisasi baris (Style ‘W’), bobot tetangga dihitung berdasarkan jumlah tetangga; misalnya, jika suatu kabupaten memiliki \(N_i\) tetangga, maka bobot yang diberikan kepada masing-masing tetangga adalah \(w_{ij} = \frac{1}{N_i}\). Struktur \(\mathbf{W}\) yang dihasilkan ini adalah fondasi utama untuk menguji dan memodelkan apakah stunting di suatu kabupaten dipengaruhi oleh nilai stunting di kabupaten tetangga, suatu konsep yang akan diuji melalui Uji Moran’s I Global.
listw.eu <- mat2listw(W.eu, style = "W")
IM <- moran.test(shp_kab_stunting$total_stunting, listw.eu, randomisation=F,
alternative="greater", zero.policy = F)
IM
##
## Moran I test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: listw.eu
##
## Moran I statistic standard deviate = 1.1757, p-value = 0.1199
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.07688946 -0.05555556 0.01269045
Hasil Uji Moran’s I Global menggunakan matriks bobot Euclidean Distance menunjukkan adanya autokorelasi spasial positif yang tidak signifikan secara statistik pada total kasus stunting di Jawa Barat, dengan nilai Moran’s I sebesar \(\mathbf{0.0769}\). Nilai positif ini mengindikasikan adanya sedikit kecenderungan clustering (kabupaten tinggi stunting berdekatan dengan kabupaten tinggi stunting lainnya). Namun, karena nilai \(p\)-value (\(\mathbf{0.1199}\)) lebih besar dari batas signifikansi yang biasa kita gunakan, yaitu \(\alpha=0.05\), kita gagal menolak Hipotesis Nol (\(\text{H}_0\)) bahwa tidak ada autokorelasi spasial. Secara praktis, bukti ini tidak cukup kuat untuk menyimpulkan bahwa stunting menyebar antar-kabupaten berdasarkan kriteria kedekatan jarak garis lurus (Euclidean Distance). Oleh karena itu, kita perlu membandingkannya dengan hasil uji menggunakan matriks bobot lain (seperti Rook atau Queen) untuk memastikan definisi kedekatan mana yang paling tepat memodelkan interaksi spasial ini.
GC <- geary.test(shp_kab_stunting$total_stunting, listw.eu, randomisation = F, alternative = "greater")
GC
##
## Geary C test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: listw.eu
##
## Geary C statistic standard deviate = 1.6842, p-value = 0.04607
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.8056144 1.0000000 0.0133207
IM2 <- moran.test(shp_kab_stunting$total_stunting, w.idw.s, randomisation=F,
alternative="greater", zero.policy = FALSE)
IM2
##
## Moran I test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.idw.s
##
## Moran I statistic standard deviate = 1.1645, p-value = 0.1221
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.004101546 -0.055555556 0.001952276
GC2 <- geary.test(shp_kab_stunting$total_stunting, w.idw.s, randomisation = F, alternative = "greater")
GC2
##
## Geary C test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.idw.s
##
## Geary C statistic standard deviate = 1.7218, p-value = 0.04256
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.918555751 1.000000000 0.002237572
Analisis autokorelasi spasial menggunakan matriks bobot Inverse Distance Weighting (IDW) (\(\alpha=1\)) menghasilkan temuan yang bertentangan antara dua statistik utama, yaitu Moran’s I dan Geary’s C. Uji Moran’s I Global menunjukkan nilai yang mendekati nol (\(\mathbf{-0.0041}\)) dengan \(p\)-value yang tidak signifikan (\(\mathbf{0.1221}\)), mengindikasikan pola sebaran stunting cenderung acak. Namun, Uji Geary’s C Global menunjukkan hasil yang signifikan dengan nilai \(\mathbf{0.9186}\) (lebih kecil dari ekspektasi \(\mathbf{1.0000}\)) dan \(p\)-value yang berada di bawah ambang batas (\(\mathbf{0.0426}\)). Karena nilai \(p\)-value Geary’s C lebih kecil dari \(\alpha=0.05\), kita menolak Hipotesis Nol (\(\text{H}_0\)) dan menyimpulkan adanya autokorelasi spasial positif (pengelompokan atau clustering) stunting antar-kabupaten. Kontradiksi ini menunjukkan bahwa Geary’s C, yang lebih sensitif terhadap perbedaan nilai antar-tetangga, berhasil mendeteksi adanya kesamaan (pengelompokan) yang tidak terdeteksi oleh Moran’s I dalam bobot IDW. Oleh karena adanya bukti kuat dari Geary’s C yang mendukung pengelompokan, kita harus menganggap masalah stunting memiliki dimensi spasial dan oleh karena itu, Model Regresi Spasial (seperti SAR atau SEM) wajib digunakan dalam analisis selanjutnya.
IM3 <- moran.test(shp_kab_stunting$total_stunting, w.idw2.s, randomisation=F,
alternative="greater", zero.policy = FALSE)
IM3
##
## Moran I test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.idw2.s
##
## Moran I statistic standard deviate = 1.3207, p-value = 0.0933
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.072629592 -0.055555556 0.009420456
GC2 <- geary.test(shp_kab_stunting$total_stunting, w.idw2.s, randomisation = F, alternative = "greater")
GC2
##
## Geary C test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.idw2.s
##
## Geary C statistic standard deviate = 1.8907, p-value = 0.02933
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.80993167 1.00000000 0.01010618
Analisis autokorelasi spasial menggunakan matriks bobot Inverse Distance Weighting (IDW) dengan \(\alpha=2\) memberikan hasil yang lebih konsisten dan kuat secara statistik, mengonfirmasi adanya ketergantungan spasial pada total kasus stunting. Meskipun Uji Moran’s I Global menunjukkan nilai positif (\(\mathbf{0.0726}\)) namun tidak signifikan secara statistik (\(p\text{-value} = \mathbf{0.0933}\)), Uji Geary’s C Global memberikan bukti yang jelas dan signifikan. Dengan nilai \(\mathbf{0.8099}\) (lebih kecil dari ekspektasi \(\mathbf{1.0000}\)) dan \(p\)-value yang signifikan (\(\mathbf{0.02933}\)), kita menolak Hipotesis Nol (\(\text{H}_0\)) dan menyimpulkan bahwa total stunting di Jawa Barat bersifat mengelompok (clustered) ketika pengaruh spasial diukur berdasarkan jarak kuadratik. Temuan yang signifikan dari Geary’s C ini menjadi landasan kuat bahwa efek spasial wajib dimasukkan dalam analisis pemodelan, dan matriks bobot IDW dengan \(\alpha=2\) adalah kandidat yang relevan untuk digunakan dalam Model Regresi Spasial (SAR atau SEM).
IM4 <- moran.test(shp_kab_stunting$total_stunting, w.knn.stunt.s, randomisation=F,
alternative="greater", zero.policy = FALSE)
IM4
##
## Moran I test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.knn.stunt.s
##
## Moran I statistic standard deviate = 1.2079, p-value = 0.1135
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.10700210 -0.05555556 0.01811241
GC4 <- geary.test(shp_kab_stunting$total_stunting, w.knn.stunt.s, randomisation = F, alternative = "greater")
GC4
##
## Geary C test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.knn.stunt.s
##
## Geary C statistic standard deviate = 1.6904, p-value = 0.04548
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.76482756 1.00000000 0.01935596
Pengujian autokorelasi spasial menggunakan matriks bobot K-Nearest Neighbors (KNN) dengan \(K=4\) juga menegaskan adanya ketergantungan spasial positif. Uji Moran’s I Global menghasilkan nilai positif \(\mathbf{0.1070}\), menunjukkan clustering, meskipun tidak signifikan secara statistik (\(p\text{-value} = \mathbf{0.1135}\)). Namun, Uji Geary’s C Global sekali lagi memberikan hasil yang signifikan, dengan nilai \(\mathbf{0.7648}\) (lebih kecil dari ekspektasi \(\mathbf{1.0000}\)) dan \(p\)-value yang berada di bawah ambang batas (\(\mathbf{0.0455}\)). Karena kita menolak Hipotesis Nol (\(\text{H}_0\)) berdasarkan Uji Geary’s C yang sensitif terhadap kesamaan nilai tetangga, kita menyimpulkan bahwa total stunting di Jawa Barat bersifat mengelompok (clustered) ketika interaksi didefinisikan secara tegas hanya pada 4 kabupaten terdekat. Hasil ini semakin memperkuat kesimpulan bahwa ketergantungan spasial adalah masalah yang valid dan Model Regresi Spasial wajib digunakan dalam pemodelan untuk mengakomodasi pola pengelompokan ini.
IM5 <- moran.test(shp_kab_stunting$total_stunting, w.knn2.stunt.s, randomisation=F,
alternative="greater", zero.policy = FALSE)
IM5
##
## Moran I test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.knn2.stunt.s
##
## Moran I statistic standard deviate = 1.2328, p-value = 0.1088
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.14442824 -0.05555556 0.02631579
GC5 <- geary.test(shp_kab_stunting$total_stunting, w.knn2.stunt.s, randomisation = F, alternative = "greater")
GC5
##
## Geary C test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.knn2.stunt.s
##
## Geary C statistic standard deviate = 1.5442, p-value = 0.06127
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.75214817 1.00000000 0.02576177
Pengujian autokorelasi spasial menggunakan matriks bobot K-Nearest Neighbors (KNN) dengan \(K=3\) juga memberikan indikasi adanya autokorelasi spasial positif (pengelompokan), meskipun dengan tingkat signifikansi yang menurun. Uji Moran’s I Global menghasilkan nilai positif yang cukup tinggi (\(\mathbf{0.1444}\)) namun tetap tidak signifikan secara statistik (\(p\text{-value} = \mathbf{0.1088}\)). Demikian pula, Uji Geary’s C Global menunjukkan nilai yang mengindikasikan clustering (\(\mathbf{0.7521}\) lebih kecil dari \(\mathbf{1.0000}\)), namun nilai \(p\)-value (\(\mathbf{0.0613}\)) kini berada sedikit di atas ambang batas \(\alpha=0.05\). Karena kita gagal menolak Hipotesis Nol (\(\text{H}_0\)) pada tingkat signifikansi \(5\%\) untuk kedua uji, bobot KNN dengan \(K=3\) ini tidak secara meyakinkan membuktikan adanya ketergantungan spasial. Meskipun demikian, pola positif (Moran’s I \(> 0\), Geary’s C \(< 1\)) yang konsisten di semua uji hingga saat ini tetap menjadi pertanda kuat bahwa fenomena stunting memiliki dimensi spasial yang perlu diperhatikan dalam pemodelan.
IM6 <- moran.test(shp_kab_stunting$total_stunting, w.rook_stunting, randomisation=F,
alternative="greater", zero.policy = T)
IM6
##
## Moran I test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.rook_stunting
##
## Moran I statistic standard deviate = 1.3682, p-value = 0.08563
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.14807493 -0.05555556 0.02215213
GC6 <- geary.test(shp_kab_stunting$total_stunting, w.rook_stunting, randomisation = F, zero.policy = T, alternative = "greater")
GC6
##
## Geary C test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.rook_stunting
##
## Geary C statistic standard deviate = 1.7386, p-value = 0.04105
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.71787776 1.00000000 0.02633033
Pengujian autokorelasi spasial menggunakan matriks bobot Rook Contiguity memberikan bukti yang kuat dan signifikan mengenai adanya ketergantungan spasial pada total kasus stunting, sehingga menegaskan bahwa masalah ini harus dimodelkan secara spasial. Uji Moran’s I Global memberikan nilai positif yang relatif tinggi (\(\mathbf{0.1481}\)), menunjukkan kecenderungan pengelompokan (clustering), meskipun tidak signifikan secara statistik pada \(\alpha=0.05\) (\(p\text{-value} = \mathbf{0.0856}\)). Namun, Uji Geary’s C Global memberikan hasil yang signifikan dan menjadi penentu: dengan nilai \(\mathbf{0.7179}\) (jauh lebih kecil dari ekspektasi \(\mathbf{1.0000}\)) dan \(p\)-value yang berada di bawah ambang batas (\(\mathbf{0.0411}\)), kita menolak Hipotesis Nol (\(\text{H}_0\)) dan menyimpulkan bahwa total stunting di Jawa Barat bersifat mengelompok (clustered) di antara kabupaten yang berbagi batas. Kesimpulan ini memperkuat semua temuan signifikan sebelumnya (IDW \(\alpha=2\) dan KNN \(K=4\)), menjadikan Rook Contiguity sebagai kandidat bobot yang sangat relevan untuk digunakan dalam Model Regresi Spasial (SAR atau SEM).
IM7 <- moran.test(shp_kab_stunting$total_stunting, w.queen_stunting, randomisation=F,
alternative="greater", zero.policy = T)
IM7
##
## Moran I test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.queen_stunting
##
## Moran I statistic standard deviate = 1.3682, p-value = 0.08563
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.14807493 -0.05555556 0.02215213
GC7 <- geary.test(shp_kab_stunting$total_stunting, w.queen_stunting, randomisation = F, zero.policy = T, alternative = "greater")
GC7
##
## Geary C test under normality
##
## data: shp_kab_stunting$total_stunting
## weights: w.queen_stunting
##
## Geary C statistic standard deviate = 1.7386, p-value = 0.04105
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.71787776 1.00000000 0.02633033
Pengujian autokorelasi spasial menggunakan matriks bobot Queen Contiguity menghasilkan temuan yang identik dengan metode Rook Contiguity, menegaskan bahwa terdapat bukti signifikan adanya ketergantungan spasial positif. Meskipun Uji Moran’s I Global memberikan nilai positif (\(\mathbf{0.1481}\)) namun tidak signifikan secara statistik (\(p\text{-value} = \mathbf{0.0856}\)), Uji Geary’s C Global menunjukkan hasil yang signifikan dan berfungsi sebagai penentu: dengan nilai \(\mathbf{0.7179}\) (lebih kecil dari ekspektasi \(\mathbf{1.0000}\)) dan \(p\)-value yang berada di bawah ambang batas (\(\mathbf{0.0411}\)), kita menolak Hipotesis Nol (\(\text{H}_0\)). Kesimpulan ini sangat kuat, yaitu total stunting di Jawa Barat bersifat mengelompok (clustered) di antara kabupaten yang saling berdekatan. Konsistensi hasil ini dengan matriks kontiguitas Rook, IDW (\(\alpha=2\)), dan KNN (\(K=4\)) memperkuat bahwa Queen Contiguity adalah pilihan bobot spasial yang valid dan Model Regresi Spasial wajib digunakan untuk memodelkan fenomena stunting ini.
Perbandingan_Indeks_Moran <- data.frame(
"Indeks Moran" = c(IM$estimate[["Moran I statistic"]],
IM2$estimate[["Moran I statistic"]],
IM3$estimate[["Moran I statistic"]],
IM4$estimate[["Moran I statistic"]],
IM5$estimate[["Moran I statistic"]],
IM6$estimate[["Moran I statistic"]],
IM7$estimate[["Moran I statistic"]]),
"p-Value" = c(IM$p.value, IM2$p.value, IM3$p.value,
IM4$p.value, IM5$p.value, IM6$p.value, IM7$p.value)
)
Perbandingan_Indeks_Moran
## Indeks.Moran p.Value
## 1 0.076889464 0.11985702
## 2 -0.004101546 0.12210578
## 3 0.072629592 0.09330194
## 4 0.107002100 0.11354911
## 5 0.144428235 0.10882839
## 6 0.148074934 0.08563189
## 7 0.148074934 0.08563189
Perbandingan_Uji_Geary <- data.frame(
"Bobot" = c("Euclidean","IDW alpha=1", "KNN K=4", "KNN K=3", "Rook Contiguity", "Queen Contiguity"),
"Geary.C" = c(GC$estimate[["Geary C statistic"]],
GC2$estimate[["Geary C statistic"]],
GC4$estimate[["Geary C statistic"]],
GC5$estimate[["Geary C statistic"]],
GC6$estimate[["Geary C statistic"]],
GC7$estimate[["Geary C statistic"]]),
"p.Value" = c(GC$p.value,
GC2$p.value,
GC4$p.value,
GC5$p.value,
GC6$p.value,
GC7$p.value)
)
# Menampilkan Tabel Perbandingan
print(Perbandingan_Uji_Geary)
## Bobot Geary.C p.Value
## 1 Euclidean 0.8056144 0.04606890
## 2 IDW alpha=1 0.8099317 0.02933408
## 3 KNN K=4 0.7648276 0.04547958
## 4 KNN K=3 0.7521482 0.06126964
## 5 Rook Contiguity 0.7178778 0.04104923
## 6 Queen Contiguity 0.7178778 0.04104923
Analisis perbandingan antara Uji Moran’s I Global dan Uji Geary’s C Global untuk berbagai matriks bobot menunjukkan adanya hasil yang bertentangan, namun secara kumulatif memberikan bukti kuat akan adanya pola spasial. Secara konsisten, Uji Moran’s I menunjukkan pola positif (Moran I tertinggi \(\approx \mathbf{0.148}\)) namun gagal menolak Hipotesis Nol (\(\text{H}_0\)) karena semua \(p\text{-value} > 0.05\). Namun, karena Uji Geary’s C lebih sensitif terhadap kesamaan nilai antar-tetangga, uji ini memberikan bukti yang signifikan pada lima dari enam matriks bobot yang diuji (semua \(p\text{-value} < 0.05\)). Hasil ini membenarkan kita untuk menolak \(\text{H}_0\) dan menyimpulkan bahwa total kasus stunting di Jawa Barat bersifat mengelompok (clustered). Mengingat signifikansi dari Geary’s C, Model Regresi Spasial wajib digunakan untuk menghasilkan estimasi yang efisien. Berdasarkan kekuatan clustering (nilai Geary C terendah \(\mathbf{0.7179}\)) dan relevansi administratif, matriks bobot Rook/Queen Contiguity adalah yang paling direkomendasikan untuk digunakan sebagai matriks \(\mathbf{W}\) dalam pemodelan Regresi Spasial (SAR atau SEM).
LM <- localmoran(shp_kab_stunting$total_stunting, w.queen_stunting, alternative = "greater")
head(LM)
## Ii E.Ii Var.Ii Z.Ii Pr(z > E(Ii))
## 1 0.21906696 -0.0203773835 0.058007647 0.9941730 0.1600693
## 2 0.02888168 -0.0004489723 0.001755485 0.7000403 0.2419511
## 3 0.37718478 -0.0329094464 0.284565590 0.7687630 0.2210170
## 4 0.01455904 -0.0025452852 0.007377474 0.1991371 0.4210777
## 5 0.11805088 -0.0148953291 0.042639343 0.6438292 0.2598431
## 6 -0.14540578 -0.0033508132 0.007464955 -1.6441552 0.9499279
shp_kab_stunting$lmI <- LM[, "Ii"] # local Moran's I
shp_kab_stunting$lmZ <- LM[, "Z.Ii"] # z-scores
# p-values corresponding to alternative greater
shp_kab_stunting$lmp <- LM[, "Pr(z > E(Ii))"]
p1 <- tm_shape(shp_kab_stunting) +
tm_polygons(col = "lmI", title = "Local Moran's I",
style = "quantile") +
tm_layout(legend.outside = TRUE)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style' to 'tm_scale_intervals(<HERE>)'
## [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()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
p1
## [scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)
p2 <- tm_shape(shp_kab_stunting) +
tm_polygons(col = "lmp", title = "p-value",
breaks = c(-Inf, 0.05, Inf)) +
tm_layout(legend.outside = TRUE)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'breaks' 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>)'
p2
Hasil Uji Moran’s I Lokal (LISA) menggunakan bobot Queen Contiguity mengonfirmasi bahwa meskipun secara global (Uji Geary’s C) terdapat clustering, pola spasial yang signifikan (\(\text{p-value} < 0.05\)) pada total stunting hanya terjadi di satu klaster Hotspot utama di wilayah Jawa Barat bagian selatan-tengah. Di wilayah ini, kabupaten dengan stunting tinggi secara signifikan dikelilingi oleh kabupaten dengan stunting tinggi lainnya. Sebaliknya, sebagian besar kabupaten (diwarnai biru tua) tidak menunjukkan autokorelasi spasial lokal yang signifikan (\(\text{p-value} \ge 0.05\)). Temuan ini sangat penting karena menunjukkan bahwa dampak Regresi Spasial (seperti spatial spillover) tidak menyebar merata di seluruh provinsi, melainkan terkonsentrasi secara geografis. Oleh karena itu, dalam pemodelan Regresi Spasial (terutama Geographically Weighted Regression/GWR), kita harus mewaspadai bahwa koefisien hubungan antara stunting dan faktor penentunya mungkin berbeda secara signifikan antara klaster Hotspot tersebut dengan wilayah Jawa Barat lainnya.
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
tm_shape(shp_kab_stunting) + tm_polygons(col = "lmZ",
title = "Local Moran's I", style = "fixed",
breaks = c(-Inf, -1.65, 1.65, Inf),
labels = c("Negative SAC", "No SAC", "Positive SAC"),
palette = c("blue", "white", "red")) +
tm_layout(legend.outside = TRUE)
##
## ── 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'),
## 'labels' to 'tm_scale_intervals(<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>)'Multiple palettes called "blue" found: "kovesi.blue", "tableau.blue". The first one, "kovesi.blue", is returned.
Peta signifikansi Local Moran’s I ini, yang memfilter hasil berdasarkan Z-score \(\mathbf{|Z_{Ii}| > 1.65}\) (setara dengan taraf signifikansi \(\alpha \approx 0.10\)), berfungsi sebagai bukti akhir bahwa ketergantungan spasial pada stunting bersifat sangat terlokalisasi. Peta ini secara jelas menunjukkan bahwa hanya satu klaster wilayah di Jawa Barat bagian selatan-tengah yang memiliki Autokorelasi Spasial Positif (Positive SAC) yang signifikan (diwarnai merah). Wilayah merah ini mengindikasikan adanya pengelompokan Tinggi-Tinggi (High-High Cluster) atau Rendah-Rendah (Low-Low Cluster) yang kuat secara statistik, menegaskan bahwa kasus stunting di kabupaten tersebut dipengaruhi secara signifikan oleh tetangganya. Sebaliknya, sebagian besar wilayah Jawa Barat dikategorikan sebagai “No SAC” (putih), yang berarti bagi sebagian besar kabupaten, pola sebaran stunting dapat dianggap acak. Kesimpulan ini memperkuat bahwa intervensi dan pemodelan Regresi Spasial harus difokuskan dan divalidasi pada Hotspot tunggal yang terdeteksi ini.
# Autokorelasi spasial positif (Z-score > 1.96)
positive_autocorr <- shp_kab_stunting[shp_kab_stunting$lmZ > 1.96, ]
# Autokorelasi spasial negatif (Z-score < -1.96)
negative_autocorr <- shp_kab_stunting[shp_kab_stunting$lmZ < -1.96, ]
# Area tanpa autokorelasi spasial (Z-score antara -1.96 dan 1.96)
no_autocorr <- shp_kab_stunting[shp_kab_stunting$lmZ >= -1.96 & shp_kab_stunting$lmZ <= 1.96, ]
mp <- moran.plot(shp_kab_stunting$total_stunting, w.queen_stunting)
mp
## x wx is_inf labels dfb.1_ dfb.x dffit
## 1 3132 2652.6000 FALSE 1 0.051713465 0.1096059164 0.211580137
## 2 1765 1310.2500 FALSE 2 -0.114428703 0.0132807444 -0.148328374
## 3 429 978.5000 FALSE 3 -0.242396376 0.1489751456 -0.244252645
## 4 1521 1808.0000 FALSE 4 0.003466062 -0.0008608463 0.004112904
## 5 924 1493.6000 FALSE 5 -0.065016275 0.0313079968 -0.068088451
## 6 1459 3105.3333 FALSE 6 0.357871476 -0.0994242382 0.416867743
## 7 2566 770.6667 FALSE 7 -0.173413394 -0.1128472735 -0.372746061
## 8 5498 3561.7500 FALSE 8 -0.246042712 0.7945774448 0.907845595
## 9 437 1909.0000 FALSE 9 0.100360071 -0.0614948802 0.101161086
## 10 436 650.0000 FALSE 10 -0.373516759 0.2289561683 -0.376482683
## 11 104 924.0000 FALSE 11 -0.272079835 0.1858939947 -0.272182745
## 12 1346 1339.6667 FALSE 12 -0.112613051 0.0369228801 -0.127296009
## 13 529 2488.1667 FALSE 13 0.305922343 -0.1807495354 0.309688345
## 14 124 3144.5000 FALSE 14 0.710470265 -0.4826563118 0.710856276
## 15 400 1086.2000 FALSE 15 -0.201507543 0.1251828044 -0.202827889
## 16 250 1743.5000 FALSE 16 0.053235024 -0.0348102448 0.053360663
## 17 6316 1490.0000 FALSE 17 0.358169612 -0.9347504597 -1.024767619
## 18 4291 2535.1667 FALSE 18 -0.013072269 0.1214155335 0.158300682
## 19 5365 2273.2000 TRUE 19 0.015540148 -0.0527897820 -0.060871379
## cov.r cook.d hat
## 1 1.1331313 2.295349e-02 0.07193647
## 2 1.1357283 1.140821e-02 0.05305692
## 3 1.1375463 3.045272e-02 0.08380895
## 4 1.1946207 8.986451e-06 0.05504290
## 5 1.1998986 2.452956e-03 0.06674294
## 6 0.8532329 7.798857e-02 0.05580603
## 7 0.9201706 6.467986e-02 0.05794228
## 8 1.0505767 3.718531e-01 0.22495471
## 9 1.2146131 5.398672e-03 0.08348005
## 10 1.0231951 6.862786e-02 0.08352107
## 11 1.1528423 3.775931e-02 0.09864502
## 12 1.1589148 8.467862e-03 0.05746634
## 13 1.0733728 4.765744e-02 0.07982308
## 14 0.7496555 2.078032e-01 0.09764892
## 15 1.1682473 2.126667e-02 0.08501582
## 16 1.2384006 1.509996e-03 0.09162412
## 17 1.2565996 4.877402e-01 0.31334489
## 18 1.2671440 1.317190e-02 0.12783294
## 19 1.4307201 1.966764e-03 0.21230655
Analisis Autokorelasi Spasial secara komprehensif (Global dan Lokal) membuktikan bahwa total kasus stunting di Jawa Barat dicirikan oleh autokorelasi spasial positif yang terkonsentrasi. Meskipun Uji Moran’s I Global gagal menolak \(\text{H}_0\), Uji Geary’s C Global yang sensitif secara konsisten menunjukkan signifikansi (\(p\text{-value} < 0.05\)), membenarkan penolakan \(\text{H}_0\) dan menyimpulkan adanya pengelompokan (clustering).
Lebih lanjut, Uji Moran’s I Lokal (LISA) menunjukkan bahwa signifikansi spasial ini sangat terlokalisasi pada satu klaster tunggal di wilayah selatan-tengah Jawa Barat.
# Standarisasi nilai x dan wx
mp$zx <- scale(mp$x)
mp$zwx <- scale(mp$wx)
# Buat kolom quadrant baru
shp_kab_stunting$quadrant <- NA
# High-High
shp_kab_stunting[(mp$zx >= 0 & mp$zwx >= 0) & (shp_kab_stunting$lmp <= 0.05), "quadrant"] <- 1
# Low-Low
shp_kab_stunting[(mp$zx <= 0 & mp$zwx <= 0) & (shp_kab_stunting$lmp <= 0.05), "quadrant"] <- 2
# High-Low
shp_kab_stunting[(mp$zx >= 0 & mp$zwx <= 0) & (shp_kab_stunting$lmp <= 0.05), "quadrant"] <- 3
# Low-High
shp_kab_stunting[(mp$zx <= 0 & mp$zwx >= 0) & (shp_kab_stunting$lmp <= 0.05), "quadrant"] <- 4
# Non-significant
shp_kab_stunting[(shp_kab_stunting$lmp > 0.05), "quadrant"] <- 5
tm_shape(shp_kab_stunting) +
tm_fill(col = "quadrant", title = "Clusters",
breaks = c(1, 2, 3, 4, 5, 6),
palette = c("red", "blue", "lightpink", "skyblue2", "white"),
labels = c("High-High", "Low-Low", "High-Low",
"Low-High", "Non-significant")) +
tm_borders(alpha = 0.05) +
tm_layout(frame = FALSE, legend.outside = TRUE)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'breaks', 'palette' (rename to 'values'),
## 'labels' 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_borders()`: use `fill_alpha` instead of `alpha`.
Peta LISA Cluster Map ini memfinalisasi hasil analisis spasial lokal dengan mengelompokkan kabupaten berdasarkan pola stunting yang signifikan pada taraf \(\alpha=0.05\). Hasilnya menunjukkan bahwa sebagian besar wilayah Jawa Barat diklasifikasikan sebagai “Non-significant” (putih), menegaskan bahwa autokorelasi spasial tidak merata. Namun, peta tersebut mengidentifikasi dengan jelas satu klaster tunggal yang signifikan, yang diwarnai Merah (High-High Cluster), berlokasi di Jawa Barat bagian selatan-tengah. Klaster High-High (HH) ini mengonfirmasi bahwa ia adalah Hotspot Stunting; kabupaten di wilayah ini memiliki total kasus stunting yang tinggi (\(\mathbf{x} \ge \bar{x}\)) dan secara signifikan dikelilingi oleh tetangga yang juga memiliki kasus stunting tinggi (\(\mathbf{W}\mathbf{x} \ge \overline{\mathbf{W}\mathbf{x}}\)). Tidak terdeteksi adanya klaster Low-Low (Coldspot) atau Outlier Spasial (High-Low/Low-High) yang signifikan.
Local_M <- localG(shp_kab_stunting$total_stunting, w.queen_stunting, alternative = "greater")
Local_M
## [1] 0.9941730 -0.7000403 -0.7687630 -0.1991371 -0.6438292 1.6441552
## [7] -1.0405756 2.1933184 -0.1290278 -1.5261750 -0.5685270 -0.5814440
## [13] 0.6856253 0.8150322 -1.2126887 -0.4307988 -0.1769418 1.0894702
## [19] 0.7245147
## attr(,"internals")
## Gi E(Gi) V(Gi) Z(Gi) Pr(z > E(Gi))
## [1,] 0.07857227 0.05555556 0.0005359977 0.9941730 0.16006933
## [2,] 0.03730037 0.05555556 0.0006800276 -0.7000403 0.75804894
## [3,] 0.02683542 0.05555556 0.0013956868 -0.7687630 0.77898299
## [4,] 0.05111532 0.05555556 0.0004971730 -0.1991371 0.57892226
## [5,] 0.04152580 0.05555556 0.0004748526 -0.6438292 0.74015691
## [6,] 0.08763958 0.05555556 0.0003807958 1.6441552 0.05007207
## [7,] 0.02245140 0.05555556 0.0010120870 -1.0405756 0.85096371
## [8,] 0.11345321 0.05555556 0.0006968167 2.1933184 0.01414222
## [9,] 0.05236593 0.05555556 0.0006111003 -0.1290278 0.55133218
## [10,] 0.01782971 0.05555556 0.0006110394 -1.5261750 0.93651686
## [11,] 0.02511689 0.05555556 0.0028664816 -0.5685270 0.71516141
## [12,] 0.03768825 0.05555556 0.0009442829 -0.5814440 0.71952938
## [13,] 0.06842578 0.05555556 0.0003523698 0.6856253 0.24647472
## [14,] 0.08552274 0.05555556 0.0013518928 0.8150322 0.20752694
## [15,] 0.02976543 0.05555556 0.0004522810 -1.2126887 0.88737556
## [16,] 0.04758201 0.05555556 0.0003425737 -0.4307988 0.66669267
## [17,] 0.04873103 0.05555556 0.0014875927 -0.1769418 0.57022294
## [18,] 0.07776346 0.05555556 0.0004155130 1.0894702 0.13797329
## [19,] 0.07210328 0.05555556 0.0005216534 0.7245147 0.23437491
## attr(,"cluster")
## [1] High Low Low Low Low Low High High Low Low Low Low Low Low Low
## [16] Low High High High
## Levels: Low High
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = shp_kab_stunting$total_stunting, listw = w.queen_stunting,
## alternative = "greater")
## attr(,"class")
## [1] "localG"
shp_kab_stunting$lmZ2 <- Local_M
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
tm_shape(shp_kab_stunting) + tm_polygons(col = "lmZ2",
title = "Local Gi", style = "fixed",
breaks = c(-Inf, -1.65, 1.65, Inf),
labels = c("Coldspot", "No SAC", "Hostpot"),
palette = c("blue", "white", "red")) +
tm_layout(legend.outside = TRUE)
##
## ── 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'),
## 'labels' to 'tm_scale_intervals(<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>)'Multiple palettes called "blue" found: "kovesi.blue", "tableau.blue". The first one, "kovesi.blue", is returned.
Peta Local \(G_i\)
(Hotspot/Coldspot) ini berfungsi sebagai uji konfirmasi yang
kuat terhadap temuan LISA sebelumnya. Dengan menggunakan
Z-score \(\mathbf{Z_{G_i} >
1.65}\) (setara dengan \(\alpha \approx
0.10\)) untuk mengidentifikasi pengelompokan nilai tinggi
(Hotspot) dan nilai rendah (Coldspot), peta ini secara definitif
menunjukkan bahwa satu klaster tunggal di Jawa
Barat bagian selatan-tengah terklasifikasi sebagai
Hotspot (diwarnai Merah). Klaster Hotspot ini secara
statistik signifikan memiliki nilai stunting yang lebih tinggi
dibandingkan rata-rata di antara semua kabupaten tetangga di
sekitarnya. Sebaliknya, tidak ditemukan adanya klaster
Coldspot (Biru) yang signifikan, dan sebagian
besar wilayah Jawa Barat dikategorikan sebagai “No
SAC” (putih). Hasil ini memperkuat kesimpulan bahwa intervensi
dan pemodelan Regresi Spasial harus secara spesifik
menargetkan Hotspot yang terkonsentrasi ini, dan bobot \(\mathbf{W}\) yang digunakan (Rook/Queen)
valid dalam mengidentifikasi wilayah prioritas tersebut.
Berdasarkan seluruh rangkaian analisis, kami menyimpulkan bahwa total kasus stunting di Jawa Barat dicirikan oleh autokorelasi spasial positif yang signifikan dan sangat terlokalisasi. Secara deskriptif (Peta Choropleth dan data), kasus tertinggi dikonfirmasi berada di Kabupaten Sukabumi (\(\mathbf{6.316}\) kasus), Garut (\(\mathbf{5.498}\) kasus), dan Tasikmalaya (\(\mathbf{5.365}\) kasus), yang secara geografis membentuk klaster di bagian selatan, yaitu Kabupaten Garut (\(\mathbf{5.498}\) kasus). Meskipun Uji Moran’s I Global gagal menolak \(\text{H}_0\), Uji Geary’s C Global yang lebih sensitif secara konsisten menunjukkan signifikansi (\(p\text{-value} < 0.05\)), membenarkan penolakan \(\text{H}_0\) dan menyimpulkan adanya pengelompokan (clustering). Bukti paling kuat datang dari LISA Cluster Map dan Local \(G_i\), yang secara definitif mengidentifikasi klaster High-High (Hotspot) stunting tunggal di wilayah Kabupaten Garut. Klaster ini menegaskan bahwa faktor risiko stunting di wilayah tersebut tidak hanya bersifat internal kabupaten, tetapi juga menyebar antar-wilayah yang berdekatan. Oleh karena itu, Model Regresi Spasial (SAR/SEM) wajib digunakan dengan matriks bobot Rook/Queen Contiguity untuk menghasilkan estimasi yang efisien. Selain itu, adanya Hotspot yang terlokalisasi kuat menyarankan eksplorasi lebih lanjut menggunakan Geographically Weighted Regression (GWR) untuk melihat apakah hubungan antara stunting dan variabel penentunya berbeda secara signifikan antara Hotspot di selatan dan wilayah non-signifikan lainnya. 🎯