Tugas Praktikum 9 : Spatial Autocorrelation
Rosita Ria Rusesta - M0501251016
2026-04-27
Pembelajaran Praktikum 9
library(sf)## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(readxl)
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(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, unionlibrary(ggplot2)
library(raster)## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectlibrary(sp)
library(spdep)
library(patchwork)##
## Attaching package: 'patchwork'## The following object is masked from 'package:raster':
##
## areaA. Joint Count Statistics
Joint count merupakan metode paling dasar dalam menentukan autokorelasi spasial antar area. Penjelasan lebih lengkap mengenai metode ini dapat dilihat pada Lizazaro (2016). Ilustrasi berikut ini diadaptasi dari Lizazaro (2016). Untuk melakukan analisis joint count, terlebih dulu kita load beberapa package pada R.
pri <- rep(1,12)
seg <- rep(0,4)
ter <- rep(1,2)
cua <- rep(0,4)
qui <- rep(1,2)
sex <- rep(0,12)
A <- matrix(c(pri, seg, ter, cua ,qui, sex), nrow=6, byrow=FALSE)
A## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 1 0 0 0 0
## [2,] 1 1 0 0 0 0
## [3,] 1 1 0 0 0 0
## [4,] 1 1 0 0 0 0
## [5,] 1 1 1 1 0 0
## [6,] 1 1 1 1 0 0Selanjutnya matriks A dikonversi menjadi raster
menggunakan fungsi raster().
rA <- raster(A)
rA## class : RasterLayer
## dimensions : 6, 6, 36 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : layer
## values : 0, 1 (min, max)Berikut ini adalah plot dari area A.
plot(rA)
text(coordinates(rA), labels=rA[ ], cex=1.5)
pA <- rasterToPolygons(rA, dissolve=FALSE)
pA## class : SpatialPolygonsDataFrame
## features : 36
## extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 1
## names : layer
## min values : 0
## max values : 1spA <- SpatialPolygons(pA@polygons)
nb1 <- poly2nb(spA, queen = T)
nb1## Neighbour list object:
## Number of regions: 36
## Number of nonzero links: 220
## Percentage nonzero weights: 16.97531
## Average number of links: 6.111111par(mai=c(0,0,0,0))
plot(spA, col='gray', border='blue')
xy <- coordinates(spA)
plot(nb1, xy, col='red', lwd=2, add=TRUE)
nb2 <- poly2nb(spA, queen = FALSE)
nb2## Neighbour list object:
## Number of regions: 36
## Number of nonzero links: 120
## Percentage nonzero weights: 9.259259
## Average number of links: 3.333333par(mai=c(0,0,0,0))
plot(spA, col='gray', border='blue')
xy <- coordinates(spA)
plot(nb2, xy, col='green', lwd=2, add=TRUE)
wl1 <- nb2listw(nb1, style='B')
wl2 <- nb2listw(nb2, style='B')
jc_test1 <- joincount.test(as.factor(pA$layer), wl1)
jc_test1##
## Join count test under nonfree sampling
##
## data: as.factor(pA$layer)
## weights: wl1
##
## Std. deviate for 0 = 5.1529, p-value = 1.282e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 53.00000 33.17460 14.80263
##
##
## Join count test under nonfree sampling
##
## data: as.factor(pA$layer)
## weights: wl1
##
## Std. deviate for 1 = 4.7634, p-value = 9.52e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 37.00000 20.95238 11.34999jc_test2 <- joincount.test(as.factor(pA$layer), wl2)
jc_test2##
## Join count test under nonfree sampling
##
## data: as.factor(pA$layer)
## weights: wl2
##
## Std. deviate for 0 = 5.4677, p-value = 2.28e-08
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 30.000000 18.095238 4.740611
##
##
## Join count test under nonfree sampling
##
## data: as.factor(pA$layer)
## weights: wl2
##
## Std. deviate for 1 = 5.1203, p-value = 1.525e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 22.000000 11.428571 4.262554Berdasarkan kedua output di atas, dengan p-value yang sangat kecil,
artinya kita dapat menolak hipotesis nol yang menyatakan bahwa tidak
terdapat autokorelasi. Sesuai dengan output di
atas, alternative hypothesis: greater, artinya kita dapat
menyimpulkan bahwa terdapat cukup bukti untuk menyatakan bahwa terdapat
autokorelasi pada taraf nyata 5%.
Pengujian hipotesis dapat pula dilakukan dengan melibatkan algoritma monte carlo seperti di bawah ini.
Global Autocorrelation
- Moran’s I
Sebagai ilustrasi, kita masih akan menggunakan data yang sama.
Seandainya kita ingin menguji autokorelasi menggunakan pendekatan indeks
moran, maka kita dapat menggunakan fungsi moran.test().
I1 <- moran.test(pA$layer,wl1)
I1##
## Moran I test under randomisation
##
## data: pA$layer
## weights: wl1
##
## Moran I statistic standard deviate = 7.4654, p-value = 4.151e-14
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.624090909 -0.028571429 0.007643049Sama halnya seperti yang telah diperlihatkan pada uji dengan metode joint count, uji moran juga dapat dilakukan dengan melibatkan simulasi monte carlo.
set.seed(123)
MC<- moran.mc(pA$layer, wl1, nsim=599)
# View results (including p-value)
MC##
## Monte-Carlo simulation of Moran I
##
## data: pA$layer
## weights: wl1
## number of simulations + 1: 600
##
## statistic = 0.62409, observed rank = 600, p-value = 0.001667
## alternative hypothesis: greater- Global Geary’s C
Geary’s C merupakan alternatif dari indeks Moran, yang memiliki nilai antara 0 s.d 2. Nilai 0 menunjukkan autokorelasi positif, 1 menunjukkan tidak ada autokorelasi, dan 2 menunjukkan autokorelasi negatif.
C1 <- geary.test(pA$layer,wl1)
C1##
## Geary C test under randomisation
##
## data: pA$layer
## weights: wl1
##
## Geary C statistic standard deviate = 7.4869, p-value = 3.526e-14
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.357954545 1.000000000 0.007354046GS1 <- geary.mc(pA$layer, wl1, nsim=599)
GS1##
## Monte-Carlo simulation of Geary C
##
## data: pA$layer
## weights: wl1
## number of simulations + 1: 600
##
## statistic = 0.35795, observed rank = 1, p-value = 0.001667
## alternative hypothesis: greaterLocal Autocorrelation
- Local Moran’s I
Pendekatan ini termasuk ke dalam Local Indicators for Spatial Association (LISA), yang mengindentifikasi autokorelasi pada tingkat lokal.
oid <- order(pA$layer)
resI <- localmoran(pA$layer, wl1)
head(resI)## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 3.75 -0.10714286 3.500600 2.0615528 0.03925033
## 2 1.75 -0.17857143 5.469688 0.8246211 0.40958672
## 3 0.40 -0.11428571 3.547275 0.2730593 0.78480760
## 4 4.00 -0.11428571 3.547275 2.1844747 0.02892738
## 5 4.00 -0.11428571 3.547275 2.1844747 0.02892738
## 6 2.40 -0.06857143 2.270256 1.6383560 0.10134744pA$z.li <- resI[,4]
pA$pvalue <- resI[,5]
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(pA, zcol="z.li", col.regions=lm.palette(20), main="Local Moran")
moran.plot(pA$layer,wl1)
- Getis-Ord Gi
Menurut Mendez (2020), pendekatan Getis-ord Gi dapat membantu mengidentifikasi pola penggerombolan berdasarkan ukuran autokorelasi pada level lokal.
local_g <- localG(pA$layer, wl1)
local_g## [1] 2.0615528 0.8246211 -0.2730593 -2.1844747 -2.1844747 -1.6383560
## [7] 2.7487371 1.2598378 -0.5233637 -2.9126330 -2.9126330 -2.1844747
## [13] 2.7487371 1.2598378 -0.5233637 -2.9126330 -2.9126330 -2.1844747
## [19] 2.7487371 2.0615528 1.0694824 -1.3197868 -2.1162099 -2.1844747
## [25] 2.7487371 2.8632678 2.0615528 -0.3435921 -1.3197868 -2.1844747
## [31] 2.0615528 2.7487371 2.7487371 0.8246211 -0.2730593 -1.6383560
## attr(,"internals")
## Gi E(Gi) V(Gi) Z(Gi) Pr(z != E(Gi))
## [1,] 0.2000000 0.08571429 0.003073229 2.0615528 0.039250330
## [2,] 0.2000000 0.14285714 0.004801921 0.8246211 0.409586724
## [3,] 0.1250000 0.14285714 0.004276711 -0.2730593 0.784807601
## [4,] 0.0000000 0.14285714 0.004276711 -2.1844747 0.028927383
## [5,] 0.0000000 0.14285714 0.004276711 -2.1844747 0.028927383
## [6,] 0.0000000 0.08571429 0.002737095 -1.6383560 0.101347443
## [7,] 0.3333333 0.14285714 0.004801921 2.7487371 0.005982535
## [8,] 0.3333333 0.22857143 0.006914766 1.2598378 0.207727870
## [9,] 0.1875000 0.22857143 0.006158463 -0.5233637 0.600721155
## [10,] 0.0000000 0.22857143 0.006158463 -2.9126330 0.003583956
## [11,] 0.0000000 0.22857143 0.006158463 -2.9126330 0.003583956
## [12,] 0.0000000 0.14285714 0.004276711 -2.1844747 0.028927383
## [13,] 0.3333333 0.14285714 0.004801921 2.7487371 0.005982535
## [14,] 0.3333333 0.22857143 0.006914766 1.2598378 0.207727870
## [15,] 0.1875000 0.22857143 0.006158463 -0.5233637 0.600721155
## [16,] 0.0000000 0.22857143 0.006158463 -2.9126330 0.003583956
## [17,] 0.0000000 0.22857143 0.006158463 -2.9126330 0.003583956
## [18,] 0.0000000 0.14285714 0.004276711 -2.1844747 0.028927383
## [19,] 0.3333333 0.14285714 0.004801921 2.7487371 0.005982535
## [20,] 0.4000000 0.22857143 0.006914766 2.0615528 0.039250330
## [21,] 0.3125000 0.22857143 0.006158463 1.0694824 0.284852347
## [22,] 0.1250000 0.22857143 0.006158463 -1.3197868 0.186906206
## [23,] 0.0625000 0.22857143 0.006158463 -2.1162099 0.034326960
## [24,] 0.0000000 0.14285714 0.004276711 -2.1844747 0.028927383
## [25,] 0.3333333 0.14285714 0.004801921 2.7487371 0.005982535
## [26,] 0.4666667 0.22857143 0.006914766 2.8632678 0.004192960
## [27,] 0.4000000 0.22857143 0.006914766 2.0615528 0.039250330
## [28,] 0.2000000 0.22857143 0.006914766 -0.3435921 0.731153040
## [29,] 0.1250000 0.22857143 0.006158463 -1.3197868 0.186906206
## [30,] 0.0000000 0.14285714 0.004276711 -2.1844747 0.028927383
## [31,] 0.2000000 0.08571429 0.003073229 2.0615528 0.039250330
## [32,] 0.3333333 0.14285714 0.004801921 2.7487371 0.005982535
## [33,] 0.3333333 0.14285714 0.004801921 2.7487371 0.005982535
## [34,] 0.2000000 0.14285714 0.004801921 0.8246211 0.409586724
## [35,] 0.1250000 0.14285714 0.004276711 -0.2730593 0.784807601
## [36,] 0.0000000 0.08571429 0.002737095 -1.6383560 0.101347443
## attr(,"cluster")
## [1] High High Low Low Low Low High High Low Low Low Low High High Low
## [16] Low Low Low High High Low Low Low Low High High High High Low Low
## [31] High High High High Low Low
## Levels: Low High
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = pA$layer, listw = wl1)
## attr(,"class")
## [1] "localG"pA$localg <- as.numeric(local_g)
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(pA, zcol="localg", col.regions=lm.palette(20), main="Local Gi")
Exercise
A. Persiapan Data
Data yang digunakan dalam praktik ini merupakan data sekunder berupa data kependudukan yang diperoleh dari Badan Pusat Statistik (BPS) dan peta shapefile berisi peta Provinsi Jawa Barat. Berdasarkan data tersebut, tahap awal yaitu proses menggabungkan data bentuk geometri wilayah (peta shp Jawa Barat) dengan data statistik atribut dalam format xlsx. Dataset mencakup informasi spasial dan tabular tingkat kabupaten/kota di Provinsi Jawa Barat untuk tahun 2015 dan 2016. Praktik ini melibatkan variabel utama, yaitu %Kemiskinan tahun 2015.
# load data
shp_jabar <- st_read("petaJabar2/Jabar2.shp")## Reading layer `Jabar2' from data source
## `/Users/rositariarusesta/Desktop/Tugas_spasial/petaJabar2/Jabar2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
## CRS: NAdata_jabar <- read_excel("data_jabar.xlsx")
data_jabar$KABKOTNO <- as.character(data_jabar$KABKOTNO)
# join
data_spasial <- shp_jabar %>%
left_join(data_jabar, by = "KABKOTNO") %>%
na.omit()Mengonstruksi matriks bobot spasial (W) yang mendefinisikan struktur interaksi antarwilayah berdasarkan kriteria jarak terdekat (k-Nearest Neighbors). Dengan menetapkan pusat gravitasi atau centroid pada setiap poligon dan mengidentifikasi k=4 tetangga terdekat, algoritma ini memastikan setiap unit observasi memiliki jumlah tetangga yang konsisten, sehingga menghindari masalah wilayah terisolasi yang sering terjadi pada data kepulauan atau wilayah dengan batas administratif yang tidak beraturan. Standardisasi baris (row-standardization) yang diterapkan pada tahap akhir berfungsi untuk menyeimbangkan pengaruh spasial antarwilayah, yang menjadi fondasi krusial bagi analisis autokorelasi seperti uji Moran’s I maupun pengembangan model regresi spasial guna memahami pola persebaran kemiskinan secara lebih akurat.
# spatial weight
coords <- st_coordinates(st_centroid(data_spasial))## Warning: st_centroid assumes attributes are constant over geometriesknn <- knearneigh(coords, k = 4)
nb_knn <- knn2nb(knn)
lw <- nb2listw(nb_knn, style = "W")Selain dalam jenis data asal berupa numerik, %kemiskinan dijadikan 2 kategori berdasarkan rata-rata. Bernilai 1 jika diatas rata-rata dan bernilai 0 untuk sebaliknya. Kemudian dilakukan visualisasi untuk mengeksplorasi sebaran spasial %kemiskinan.
threshold <- mean(data_spasial$p.miskin15)
data_spasial$miskin_bin <- ifelse(data_spasial$p.miskin15 > threshold, 1, 0)
x <- factor(data_spasial$miskin_bin)
# 1. Plot Data Asli (Spektrum Warna/Kontinu)
plot_asli <- ggplot(data = data_spasial) +
geom_sf(aes(fill = p.miskin15), color = "white", linewidth = 0.2) +
scale_fill_viridis_c(option = "magma", name = "Persentase (%)") + # Spektrum warna dari rendah ke tinggi
theme_minimal() +
labs(title = "Data Asli (Kontinu)",
subtitle = "Persentase Kemiskinan",
x = "Longitude", y = "Latitude") +
theme(legend.position = "bottom")
# 2. Plot Binary (Berdasarkan Threshold - Kode Anda sebelumnya)
plot_binary <- ggplot(data = data_spasial) +
geom_sf(aes(fill = factor(miskin_bin)), color = "white", linewidth = 0.2) +
scale_fill_manual(values = c("0" = "#E0E0E0", "1" = "#E41A1C"),
name = "Status",
labels = c("Aman", "Di Atas Rata-rata")) +
theme_minimal() +
labs(title = "Data Kategorik (Biner)",
subtitle = "Berdasarkan Threshold Rata-rata",
x = "Longitude", y = "Latitude") +
theme(legend.position = "bottom")
# 3. Menggabungkan kedua plot
(plot_asli | plot_binary) +
plot_annotation(title = "Visualisasi Sebaran %Kemiskinan Jawa Barat",
theme = theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5)))
B. Joint Count Statistics
Analisis Joint Count adalah metode statistik non-parametrik yang sangat spesifik digunakan untuk menguji autokorelasi spasial pada data biner (seperti 0 dan 1). Jika Moran’s I mengukur korelasi berdasarkan nilai kontinu, Joint Count mengukur korelasi berdasarkan kontiguitas (persentuhan) antar kategori. Sebelum dilakukan analisis, dapat dilihat visualisasi dari ketetanggaaan k-Nearest Neighbors (k=4) yang kemudian digunakan untuk matriks pembobot (W) pada Joint Count.
colors <- ifelse(data_spasial$miskin_bin == 0, "lightgrey", "red")
# 2. Plot peta dasar (poligon wilayah) dengan warna fill yang telah ditentukan
plot(st_geometry(data_spasial),
col = colors, # Mewarnai poligon sesuai vektor colors
border = "white", # Garis batas wilayah berwarna putih
lwd = 0.5, # Ketebalan garis batas wilayah
main = "Ketetanggaan KNN (k=4) Berdasarkan Status Kemiskinan")
# 3. Plot garis ketetanggaan (links) di atas peta dasar
# nb_knn: objek tetangga
# coords: koordinat centroid
plot(nb_knn, coords,
add = TRUE, # Menambahkan di atas plot yang sudah ada
col = "blue", # Warna garis tetangga
points = FALSE, # Tidak menampilkan titik centroid
lwd = 1.5) # Ketebalan garis tetangga
Dilakukan perhitungan berapa banyak sisi (edges) yang menghubungkan wilayah dengan kategori yang sama:
Same-color joins (1-1 dan 0-0): Jika jumlah hubungan 1-1 (wilayah miskin bertetangga dengan miskin) jauh lebih banyak daripada ekspektasi teoretis, maka terdapat klaster spasial positif.
Different-color joins (1-0): Jika hubungan ini yang dominan, maka terdapat pola negatif (sebaran catur/dispersi).
Dalam proses ini, akan didapatkan nilai Z-score dan P-value.
# joint count test
jc <- joincount.test(x, lw)
jc##
## Join count test under nonfree sampling
##
## data: x
## weights: lw
##
## Std. deviate for 0 = 1.3205, p-value = 0.09334
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 2.7500000 2.2500000 0.1433824
##
##
## Join count test under nonfree sampling
##
## data: x
## weights: lw
##
## Std. deviate for 1 = 2.4001, p-value = 0.008194
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 2.6250000 1.7500000 0.1329044Hasil uji Join Count menunjukkan adanya autokorelasi spasial positif yang signifikan pada kategori wilayah miskin (Kategori 1) dengan p-value=0.0082 < 0.05, yang berarti wilayah dengan tingkat kemiskinan di atas rata-rata cenderung mengelompok secara geografis membentuk klaster kemiskinan di Jawa Barat. Sebaliknya, pada kategori wilayah aman (Kategori 0), hasil uji tidak menunjukkan signifikansi statistik (p-value=0,093 > 0.05), mengindikasikan bahwa wilayah dengan kemiskinan rendah cenderung tersebar secara lebih acak tanpa pola pengelompokan yang kuat. Hal ini memberikan informasi bahwa fenomena kemiskinan di lokasi penelitian memiliki autokorelasi spasial.
# joint count dengan simulasi monte carlo
set.seed(123)
jc_mc <- joincount.mc(x, lw, nsim = 999)
jc_mc##
## Monte-Carlo simulation of join-count statistic
##
## data: x
## weights: lw
## number of simulations + 1: 1000
##
## Join-count statistic for 0 = 2.75, rank of observed statistic = 910.5,
## p-value = 0.0895
## alternative hypothesis: greater
## sample estimates:
## mean of simulation variance of simulation
## 2.2506256 0.1318101
##
##
## Monte-Carlo simulation of join-count statistic
##
## data: x
## weights: lw
## number of simulations + 1: 1000
##
## Join-count statistic for 1 = 2.625, rank of observed statistic = 990.5,
## p-value = 0.0095
## alternative hypothesis: greater
## sample estimates:
## mean of simulation variance of simulation
## 1.7374875 0.1350512Hasil simulasi Monte Carlo memperkuat uji sebelumnya dengan pendekatan yang lebih tangguh (robust), karena tidak bergantung pada asumsi distribusi normal melainkan pada 999 kali pengacakan data secara empiris.
C. Global Autocorrelation
1. Moran’s I
Setelah sebelumnya menggunakan Joint Count untuk data biner (0 dan 1), kini beralih ke Global Moran’s I, untuk menguji autokorelasi pada data kontinu (%kemiskinan tanpa dikategorikan)
Global Moran’s I mengukur hubungan antara nilai suatu variabel di sebuah lokasi dengan nilai variabel yang sama di lokasi tetangganya. Nilainya berkisar antara -1 hingga +1:
Nilai Positif (+): Menunjukkan adanya pengelompokan (Clustering). Wilayah dengan persentase kemiskinan tinggi dikelilingi wilayah yang juga tinggi.
Nilai Negatif (-): Menunjukkan pola menyebar (Dispersed). Wilayah dengan nilai tinggi dikelilingi oleh nilai rendah (seperti papan catur).
Nilai Nol (0): Menunjukkan pola acak (Random). Tidak ada keterkaitan spasial.
#========================================
#Global Moran I
#========================================
moran.test(data_spasial$p.miskin15, lw)##
## Moran I test under randomisation
##
## data: data_spasial$p.miskin15
## weights: lw
##
## Moran I statistic standard deviate = 2.8407, p-value = 0.002251
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.34343214 -0.06250000 0.02042008#uji moran dengan simulasi monte carlo
set.seed(123)
MC<- moran.mc(data_spasial$p.miskin15, lw, nsim=599)
MC##
## Monte-Carlo simulation of Moran I
##
## data: data_spasial$p.miskin15
## weights: lw
## number of simulations + 1: 600
##
## statistic = 0.34343, observed rank = 596, p-value = 0.006667
## alternative hypothesis: greaterHasil uji Global Moran’s I menunjukkan adanya autokorelasi spasial positif yang signifikan pada persentase kemiskinan dengan nilai indeks sebesar 0,3434 ( p-value =0,0022 < 0.05), yang secara empiris diperkuat oleh simulasi Monte Carlo dengan peringkat observasi 596 dari 600 percobaan (p-value=0,0067 < 0.05). Hal ini mengonfirmasi bahwa fenomena kemiskinan di Jawa Barat tidak tersebar secara acak, melainkan membentuk pola pengelompokan (clustering) yang kuat di mana wilayah dengan intensitas kemiskinan tinggi cenderung bertetangga dengan wilayah berkarakteristik serupa.
2. Global Geary’s C
Jika Moran’s I fokus pada penyimpangan nilai dari rata-rata (mirip korelasi), Geary’s C lebih fokus pada perbedaan kuadrat antar nilai tetangga (mirip prinsip varians atau semivariogram).
Nilai 0<C<1: Menunjukkan autokorelasi spasial positif (wilayah yang bertetangga memiliki nilai yang mirip/homogen).
Nilai C=1: Menunjukkan tidak ada autokorelasi (acak).
Nilai C>1: Menunjukkan autokorelasi spasial negatif (wilayah bertetangga memiliki nilai yang sangat berbeda/heterogen).
Menghitung statistik C untuk mengetahui tingkat kemiripan nilai
persentase kemiskinan antarwilayah yang bertetangga dan membandingkan
perbedaan nilai antar-poligon yang terhubung oleh matriks
bobot lw.
#========================================
#Global Geary’s C
#========================================
C1 <- geary.test(data_spasial$p.miskin15, lw)
C1##
## Geary C test under randomisation
##
## data: data_spasial$p.miskin15
## weights: lw
##
## Geary C statistic standard deviate = 2.9909, p-value = 0.001391
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.57006259 1.00000000 0.02066342#uji Geary’s C dengan simulasi monte carlo
set.seed(123)
GS1 <- geary.mc(data_spasial$p.miskin15, lw, nsim=599)
GS1##
## Monte-Carlo simulation of Geary C
##
## data: data_spasial$p.miskin15
## weights: lw
## number of simulations + 1: 600
##
## statistic = 0.57006, observed rank = 4, p-value = 0.006667
## alternative hypothesis: greaterHasil uji Global Geary’s C menghasilkan nilai statistik 0,5700 (p-value = 0,0013 < 0.05) , yang secara signifikan berada di bawah nilai harapan 1, sehingga mengonfirmasi adanya homogenitas atau kemiripan karakteristik kemiskinan yang kuat antarwilayah bertetangga. Validasi melalui simulasi Monte Carlo memperkuat temuan ini dengan peringkat observasi 4 dari 600 percobaan yang menghasilkan nilai statistik 0,57006 (p-value=0,0067 < 0.05), membuktikan bahwa pola klaster lokal ini bukan terjadi secara kebetulan. Secara statistik, hasil ini melengkapi Moran’s I dalam menegaskan bahwa fenomena kemiskinan di Jawa Barat memiliki autokorelasi spasial,
D. Local Autocorrelation
1. Local Moran’s I
Penerapan Local Moran’s I bertujuan untuk membedah
pola autokorelasi global ke dalam unit-unit observasi yang lebih kecil
guna mengidentifikasi heterogenitas spasial pada data kemiskinan di Jawa
Barat. Dengan melakukan kalkulasi
skor Z (z.li) dan P-value pada tiap wilayah, metode ini
mampu mendeteksi keberadaan indikator asosiasi spasial lokal
atau LISA, yang secara teknis memisahkan antara wilayah
yang menjadi bagian dari klaster signifikan dengan wilayah yang bersifat
acak.
Tahapan ekstraksi data ini sangat krusial dalam alur kerja sains data karena menjadi dasar dalam pengklasifikasian kuadran Hotspot dan Coldspot, yang memungkinkan analisis lebih mendalam mengenai faktor lokal apa yang menyebabkan sebuah wilayah memiliki tingkat kemiskinan yang secara signifikan lebih tinggi atau lebih rendah dibandingkan rata-rata wilayah di sekitarnya
#========================================
#local moran
#========================================
oid <- order(data_spasial$p.miskin15)
resI <- localmoran(data_spasial$p.miskin15, lw)
head(resI)## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 1.27585405 -0.0989976131 0.3032700914 2.4965545 0.01254064
## 2 0.29366453 -0.0170499022 0.0569812902 1.3016532 0.19303497
## 3 1.02723011 -0.1343387004 0.3953921675 1.8472723 0.06470769
## 4 0.41667348 -0.0359515480 0.1178407163 1.3185326 0.18732541
## 5 -0.01348771 -0.0001924233 0.0006541132 -0.5198420 0.60317369
## 6 -0.09087487 -0.0043770537 0.0148168433 -0.7106034 0.47733004data_spasial$z.li <- resI[,4]
data_spasial$pvalue <- resI[,5]Identifikasi Klaster Signifikan: Melalui p-value, dapat diidentifikasi unit wilayah yang memiliki ketergantungan spasial signifikan. Sebagai contoh, pada unit observasi pertama, diperoleh nilai Z-score sebesar 2,49 dengan p-value=0,0125 (<0,05). Hal ini menunjukkan bahwa wilayah tersebut merupakan pusat klaster yang signifikan secara statistik, di mana karakteristik tingkat kemiskinannya memiliki kemiripan yang kuat dengan wilayah-wilayah di sekitarnya.
Pola Pengelompokan (Indeks Ii): Nilai Ii yang positif pada beberapa baris awal (seperti pada observasi 1 hingga 4) mengindikasikan adanya asosiasi spasial positif, yang berarti wilayah-wilayah tersebut dikelilingi oleh tetangga yang memiliki nilai persentase kemiskinan yang mirip
Variansi dan Distribusi Lokal: Keberadaan wilayah dengan p-value > 0,05 menunjukkan bahwa tidak semua wilayah di Jawa Barat membentuk klaster. Wilayah-wilayah tersebut dikategorikan sebagai wilayah yang memiliki pola sebaran acak secara spasial, di mana tingkat kemiskinannya tidak dipengaruhi secara signifikan oleh kondisi wilayah tetangganya.
# Buat plot
ggplot(data = data_spasial) +
geom_sf(aes(fill = z.li), color = "black", linewidth = 0.3) +
scale_fill_stepsn(
colors = c("#ffffff", "#fef0d9", "#fdbb84", "#fc8d59", "#e34a33", "#b30000"),
n.breaks = 20,
name = NULL
) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
legend.position = "right",
legend.key.height = unit(1.5, "cm") # Memanjangkan legenda
) +
labs(title = "Local Moran")
moran.plot(data_spasial$p.miskin15,lw)
2. Getis-Ord Gi
localG(data_spasial$p.miskin15, lw):
Fungsi ini menghitung statistik Gi∗ lokal. Hasil yang dikeluarkan secara
otomatis dalam bentuk Z-score. Karena Z-score mengikuti
distribusi normal standar, kita bisa langsung menentukan
signifikansinya:
Z>1.96: Hotspot signifikan (Confidence Level 95%).
Z<−1.96: Coldspot signifikan (Confidence Level 95%).
−1.96<Z<1.96: Tidak signifikan secara spasial.
#========================================
#Getis-Ord Gi
#========================================
local_g <- localG(data_spasial$p.miskin15, lw)
local_g## [1] 2.49655451 1.30165318 1.84727228 1.31853263 0.51984202 -0.71060340
## [7] -1.62044928 -1.18559891 -2.26434951 -0.81197560 0.06131463 1.52570656
## [13] -1.86072899 -1.90537281 -0.07378118 1.20868799 2.03433502
## attr(,"internals")
## Gi E(Gi) V(Gi) Z(Gi) Pr(z != E(Gi))
## [1,] 0.09148252 0.0625 0.0001347691 2.49655451 0.01254064
## [2,] 0.07798074 0.0625 0.0001414467 1.30165318 0.19303497
## [3,] 0.08363659 0.0625 0.0001309206 1.84727228 0.06470769
## [4,] 0.07812600 0.0625 0.0001404475 1.31853263 0.18732541
## [5,] 0.06864298 0.0625 0.0001396417 0.51984202 0.60317369
## [6,] 0.05405104 0.0625 0.0001413682 -0.71060340 0.47733004
## [7,] 0.04450262 0.0625 0.0001233525 -1.62044928 0.10513580
## [8,] 0.04848862 0.0625 0.0001396646 -1.18559891 0.23578073
## [9,] 0.03618569 0.0625 0.0001350507 -2.26434951 0.02355263
## [10,] 0.05293281 0.0625 0.0001388298 -0.81197560 0.41680563
## [11,] 0.06316906 0.0625 0.0001190701 0.06131463 0.95110864
## [12,] 0.08064018 0.0625 0.0001413648 1.52570656 0.12708296
## [13,] 0.04173740 0.0625 0.0001245080 -1.86072899 0.06278246
## [14,] 0.04328474 0.0625 0.0001017027 -1.90537281 0.05673163
## [15,] 0.06166940 0.0625 0.0001267340 -0.07378118 0.94118450
## [16,] 0.07588858 0.0625 0.0001226990 1.20868799 0.22678274
## [17,] 0.08607714 0.0625 0.0001343189 2.03433502 0.04191783
## attr(,"cluster")
## [1] High High High High Low High Low High Low Low Low High Low Low Low
## [16] High Low
## Levels: Low High
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = data_spasial$p.miskin15, listw = lw)
## attr(,"class")
## [1] "localG"data_spasial$localg <- as.numeric(local_g)palet_krem_merah <- c("#ffffff", "#fef0d9", "#fdbb84", "#fc8d59", "#e34a33", "#b30000")
# 3. Membuat Plot
ggplot(data = data_spasial) +
geom_sf(aes(fill = localg), color = "black", linewidth = 0.3) +
# Menggunakan scale_fill_stepsn untuk membuat efek kotak-kotak pada legenda
scale_fill_stepsn(
colors = palet_krem_merah,
n.breaks = 18, # Menyesuaikan jumlah kotak warna pada legenda di gambar
name = NULL # Menghilangkan teks judul di atas legenda
) +
theme_void() + # Menghilangkan background abu-abu dan grid koordinat
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 20),
legend.position = "right",
legend.key.height = unit(1.2, "cm"), # Memanjangkan ukuran kotak legenda
legend.text = element_text(size = 10)
) +
labs(title = "Hot Spot Analysis (Gi*)")
Perbandingan antara peta Local Moran’s I dan Hot Spot Analysis(Gi∗) terlihat adanya konsistensi spasial yang kuat, di mana wilayah bagian timur dan tenggara teridentifikasi sebagai pusat klaster kemiskinan (hotspot) dengan nilai Z-score di atas 2,00. Sementara Local Moran’s I mampu menangkap heterogenitas lokal dan batas-batas klaster secara lebih tajam, analisis Getis-Ord Gi∗ memberikan gambaran yang lebih menyeluruh mengenai intensitas sebaran fenomena tersebut. Sinkronisasi kedua hasil ini membuktikan bahwa pengelompokan kemiskinan di wilayah tersebut bukan merupakan anomali acak, melainkan sebuah struktur spasial yang nyata, sehingga wilayah dengan gradasi warna merah pekat pada kedua peta tersebut harus menjadi prioritas utama dalam intervensi kebijakan berbasis kewilayahan.