Bagian A. Konsep Teori
1. Pengertian autokorelasi spasial positif & negatif
Autokorelasi spasial positif berarti nilai variabel yang mirip (misal: tinggi) cenderung berdekatan secara geografis. Dalam kata lain, lokasi dengan nilai tinggi dikelilingi oleh lokasi lain yang juga bernilai tinggi, dan lokasi dengan nilai rendah cenderung dikelilingi lokasi rendah. Autokorelasi spasial negatif berarti pola “checkerboard”: nilai tinggi cenderung dikelilingi oleh nilai rendah dan sebaliknya. Jadi tetangga sering berbeda nilai.
Contoh sederhana:
Positif: wilayah A = 30 kasus, tetangganya B = 28 dan C = 32.
Negatif: wilayah A = 30 kasus, tetangganya B = 5 dan C = 6 (kontras tinggi-rendah).
2. Contoh fenomena
Autokorelasi positif: Penyebaran penyakit menular di kota; area padat penduduk terinfeksi lebih banyak di sekitar area sejenis, maka terbentuk hotspot. Misalnya saat pandemi covid-19 di Indonesia yang terjadi sekitar tahun 2019-2022
Autokorelasi negatif: Pola penggunaan lahan yang bergantian (checkerboard), atau lokasi industri yang sengaja memisah kualitas lingkungan rendah di area industri bersebelahan dengan permukiman bernilai lingkungan tinggi.
3. Rumus matematis dan penjelasan variabel
3.1 Moran’s I (Global)
Rumus (ditampilkan):
\[ I = \frac{N}{S_0} \cdot \frac{\sum_i \sum_j w_{ij}(x_i - \bar{x})(x_j - \bar{x})}{\sum_i (x_i - \bar{x})^2} \]
dengan:
\(N\): jumlah unit/area (mis. jumlah kecamatan).
\(x_i\): nilai variabel di unit \(i\) (mis. jumlah kasus).
\(x_bar\): rata-rata \(x\) di semua unit.
\(w_ij\): elemen matriks bobot spasial W yang menyatakan hubungan i–j (mis. 1 jika bertetangga, 0 jika tidak, atau bobot jarak).
\(S_0 = \sum_i \sum_j w_ij\): jumlah total semua bobot. Interpretasi nilai Moran’s I:
\(I > 0\): autokorelasi positif (clustering: high–high atau low–low).
\(I < 0\): autokorelasi negatif (checkerboard).
\(I \approx 0\): pola acak. Inferensi: gunakan uji permutasi atau z-score/p-value untuk menentukan signifikansi.
3.2 Geary’s C (Global)
Rumus: \[ C = \frac{(N-1)\sum_i \sum_j w_{ij}(x_i - x_j)^2}{2S_0 \sum_i (x_i - \bar{x})^2} \]
Keterangan simbol:
\(N\): jumlah unit/area (mis. jumlah kecamatan).
\(x_i\): nilai variabel di unit \(i\) (mis. jumlah kasus).
\(x_bar\): rata-rata \(x\) di semua unit.
\(w_ij\): elemen matriks bobot spasial W yang menyatakan hubungan i–j (mis. 1 jika bertetangga, 0 jika tidak, atau bobot jarak).
\(S_0 = \sum_i \sum_j w_ij\): jumlah total semua bobot. Interpretasi Geary’s C:
\(C < 1\): autokorelasi positif (tetangga mirip).
\(C > 1\): autokorelasi negatif (tetangga berbeda).
\(C = 1\): tidak ada autokorelasi (acak). Geary’s C lebih sensitif terhadap perbedaan lokal (contrast antar tetangga) dibanding Moran.
3.3 Local Moran’s I (LISA)
Rumus untuk unit i:
\[ I_i = \frac{(x_i - \bar{x})}{m_2} \sum_j w_{ij}(x_j - \bar{x}), \quad m_2 = \frac{1}{N}\sum_k (x_k - \bar{x})^2 \] dengan
\(m_2 = \frac{1}{N} * \sum_k (x_k - x_{bar})^2\).
Keterangan:
\(I_i\) mengukur kontribusi lokal unit \(i\) terhadap autokorelasi global.
Jika \(I_i > 0\) dan \(x_i\) tinggi serta tetangga tinggi, High–High.
Jika \(I_i > 0\) dan \(x_i\) rendah serta tetangga rendah, Low–Low.
Jika \(I_i < 0\), outlier (High–Low atau Low–High).
3.4 Getis–Ord \(G_i\) dan \(G_i^*\)
Rumus raw (proporsi massa tetangga):
Tanpa memasukkan lokasi i sendiri (local G_i):
\[ G_i = \frac{\sum_{j=1}^{n} w_{ij}(d) x_j}{\sum_{j=1}^{n} x_j} \]
Memasukkan lokasi i (\(G_i^*\)):
\[ G_i^* = \frac{\sum_{j=1}^{n} w_{ij}(d) x_j}{\sum_{j=1}^{n} x_j}, \quad \text{dengan } w_{ii} \neq 0 \]
Keterangan:
\(w_{ij}(d)\): bobot biner pada band jarak d (1 jika jarak <= d atau ditetapkan sebagai tetangga).
Interpretasi: nilai tinggi \(G_i^*\), hot spot; nilai rendah, cold spot.
Untuk inferensi biasanya distandarisasi menjadi z-score, z(\(G_i^*\)).
4. Perbedaan utama antara ukuran global dan lokal
Ukuran Global
Ukuran global dalam autokorelasi spasial, seperti Moran’s I dan Geary’s C, memberikan satu nilai ringkasan yang mewakili kecenderungan spasial pada keseluruhan area. Indeks ini menjawab pertanyaan umum apakah data cenderung terdistribusi secara clustered, acak, atau berpola menyebar (dispersed). Misalnya, nilai Moran’s I yang positif signifikan menunjukkan bahwa unit dengan nilai serupa (tinggi dengan tinggi, rendah dengan rendah) cenderung berdekatan, sedangkan nilai negatif mengindikasikan pola papan catur (Cliff & Ord, 1981). Kelebihan ukuran global adalah kesederhanaannya dalam mendeteksi tren umum, namun kelemahannya adalah ketidakmampuannya untuk mengidentifikasi lokasi spesifik di mana klaster atau outlier terjadi. Dengan kata lain, ukuran global hanya menangkap pola rata-rata, sehingga bisa menutupi heterogenitas antarwilayah (Anselin, 1995).Ukuran Lokal
Sebaliknya, ukuran lokal seperti Local Moran’s I (LISA) dan Getis–Ord \(G_i^*\) menghitung autokorelasi pada tingkat unit individu, sehingga dapat mengungkapkan lokasi spesifik di mana terdapat klaster (High–High, Low–Low) atau outlier (High–Low, Low–High). Analisis ini memungkinkan identifikasi pola spasial heterogen yang tidak terlihat dengan ukuran global, misalnya mendeteksi hotspot atau coldspot penyakit dalam suatu kota. Kelebihan pendekatan lokal adalah kemampuannya memberikan detail spasial, yang sangat penting untuk kebijakan berbasis wilayah. Sebagai contoh, Raza et al. (2020) menemukan bahwa meskipun prevalensi diare di Mozambik secara global terklaster, hanya provinsi tertentu seperti Tete dan Zambezia yang teridentifikasi sebagai klaster High–High melalui analisis lokal. Hal ini menunjukkan bahwa ukuran lokal melengkapi ukuran global dengan menyoroti di mana pola signifikan terjadi di ruang.
B. Analisis Data
## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE## Warning: package 'spdep' was built under R version 4.3.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.3.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Warning: package 'ggplot2' was built under R version 4.3.3## Warning: package 'viridis' was built under R version 4.3.3## Loading required package: viridisLite## Warning: package 'viridisLite' was built under R version 4.3.31. Data
Sumber data DBD per kecamatan : KASUS_DBD_DI_KECAMATAN_KOTA_BANDUNG
## Warning: package 'readxl' was built under R version 4.3.3## # A tibble: 6 Ă— 2
## Kecamatan JumlahKasus
## <chr> <dbl>
## 1 Andir 137
## 2 Antapani 259
## 3 Arcamanik 314
## 4 Astana Anyar 157
## 5 Babakan Ciparay 237
## 6 Bandung Kidul 1562. Data Spasial
Sumber data geografis : SHP_KECAMATAN_KOTA_BANDUNG
# Path shapefile
shp_path <- "C:/Users/ACER/Downloads/JAWABARAT_ADM_KEC_2021/JAWABARAT_ADM_KEC.shp"
# Baca shapefile
jabar_kec <- st_read(shp_path)## Reading layer `JAWABARAT_ADM_KEC' from data source
## `C:\Users\ACER\Downloads\JAWABARAT_ADM_KEC_2021\JAWABARAT_ADM_KEC.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 627 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: 106.3703 ymin: -7.82099 xmax: 108.8471 ymax: -5.806538
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84## [1] "Cianjur" "Kota Bandung" "Indramayu" "Majalengka"
## [5] "Bandung" "Cirebon" "Bogor" "Purwakarta"
## [9] "Bekasi" "Garut" "Kota Banjar" "Ciamis"
## [13] "Sukabumi" "Kota Bekasi" "Tasikmalaya" "Karawang"
## [17] "Kota Sukabumi" "Bandung Barat" "Kota Depok" "Subang"
## [21] "Kota Bogor" "Sumedang" "Kota Tasikmalaya" "Kuningan"
## [25] "Pangandaran" "Kota Cimahi" "Kota Cirebon"# Filter hanya data Kota Bandung
bandung_kec <- jabar_kec %>%
mutate(WADMKK = trimws(WADMKK)) %>%
filter(WADMKK == "Kota Bandung")
# Join shapefile dengan data kasus (kecocokan nama : WADMKC vs Kecamatan)
library(stringr)## Warning: package 'stringr' was built under R version 4.3.3# Samakan kapitalisasi + perbaiki penamaan
bandung_data <- bandung_kec %>%
mutate(
WADMKC = str_to_title(WADMKC),
WADMKC = ifelse(WADMKC == "astanaanyar", "Astana Anyar", WADMKC),
WADMKC = ifelse(WADMKC == "Ujungberung", "Ujung Berung", WADMKC),
WADMKC = ifelse(WADMKC == "Ujung Berung", "Ujung Berung", WADMKC) # pastikan spasi benar
) %>%
left_join(
dbd_data %>%
mutate(Kecamatan = str_to_title(Kecamatan)),
by = c("WADMKC" = "Kecamatan")
)
head(bandung_data)## Simple feature collection with 6 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 107.5597 ymin: -6.965704 xmax: 107.6913 ymax: -6.890994
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
## OBJECTID WADMKC WADMKK WADMPR SHAPE_Leng SHAPE_Area
## 1 2 Andir Kota Bandung Jawa Barat 0.12725056 0.0003467384
## 2 4 Antapani Kota Bandung Jawa Barat 0.09587712 0.0003452283
## 3 6 Arcamanik Kota Bandung Jawa Barat 0.12073520 0.0006001662
## 4 10 Astana Anyar Kota Bandung Jawa Barat 0.10506945 0.0002160417
## 5 13 Babakan Ciparay Kota Bandung Jawa Barat 0.12740210 0.0005753570
## 6 20 Bandung Kidul Kota Bandung Jawa Barat 0.12181502 0.0004082145
## JumlahKasus geometry
## 1 137 MULTIPOLYGON Z (((107.5638 ...
## 2 259 MULTIPOLYGON Z (((107.6624 ...
## 3 314 MULTIPOLYGON Z (((107.6695 ...
## 4 157 MULTIPOLYGON Z (((107.5939 ...
## 5 237 MULTIPOLYGON Z (((107.5864 ...
## 6 156 MULTIPOLYGON Z (((107.6398 ...## [1] 0## [1] "Andir" "Antapani" "Arcamanik" "Astana Anyar"
## [5] "Babakan Ciparay" "Bandung Kidul" "Bandung Kulon" "Bandung Wetan"
## [9] "Batununggal" "Bojongloa Kaler" "Bojongloa Kidul" "Buahbatu"
## [13] "Cibeunying Kaler" "Cibeunying Kidul" "Cibiru" "Cicendo"
## [17] "Cidadap" "Cinambo" "Coblong" "Gedebage"
## [21] "Kiaracondong" "Lengkong" "Mandalajati" "Panyileukan"
## [25] "Rancasari" "Regol" "Sukajadi" "Sukasari"
## [29] "Sumur Bandung" "Ujung Berung"3. Peta Chloropleth
ggplot(bandung_data) +
geom_sf(aes(fill = JumlahKasus), color = "white") +
scale_fill_viridis_c(option = "plasma") +
theme_minimal() +
labs(title = "Peta Kasus DBD Kota Bandung (2024)", fill = "Jumlah Kasus")
4. Pola Spasial yang terlihat secara Visual
plot(bandung_data["JumlahKasus"],
main = "Peta Kasus DBD Kota Bandung (2024)",
pal = viridis::viridis)
## â„ą tmap modes "plot" - "view"
## â„ą toggle with `tmap::ttm()`tm_shape(bandung_data) +
tm_polygons(
col = "JumlahKasus",
fill.scale = tm_scale(values = "Reds"),
border.col = "gray40"
) +
tm_title("Peta Kasus DBD di Kecamatan di Kota Bandung (2024)")##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [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').Bagian C. Pengukuran Autokorelasi
library(spdep)
# Buat neighbors dan spatial weights
# nb: daftar tetangga, lw: bobot spasial
nb <- poly2nb(bandung_data, queen = TRUE)## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate1. Moran’s I
moran_res <- moran.test(bandung_data$JumlahKasus, lw, randomisation = TRUE, zero.policy = TRUE)
moran_res##
## Moran I test under randomisation
##
## data: bandung_data$JumlahKasus
## weights: lw
##
## Moran I statistic standard deviate = 2.8627, p-value = 0.0021
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.29059926 -0.03448276 0.012895262. Geary’s C
geary_res <- geary.test(bandung_data$JumlahKasus, lw, randomisation = TRUE, zero.policy = TRUE)
geary_res##
## Geary C test under randomisation
##
## data: bandung_data$JumlahKasus
## weights: lw
##
## Geary C statistic standard deviate = 2.3506, p-value = 0.00937
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.72180841 1.00000000 0.014006013. Local Moran’s I (LISA)
local_res <- localmoran(bandung_data$JumlahKasus, lw, zero.policy = TRUE)
# Simpan hasil
bandung_data$lisa_I <- local_res[,1]
bandung_data$lisa_p <- local_res[,5]
# Hitung rata-rata kasus
mean_rate <- mean(bandung_data$JumlahKasus)
# Kategorisasi cluster
bandung_data$quad_sig <- NA
for (i in 1:nrow(bandung_data)) {
if (bandung_data$lisa_p[i] <= 0.05) {
if (bandung_data$JumlahKasus[i] > mean_rate & lag.listw(lw, bandung_data$JumlahKasus, zero.policy=TRUE)[i] > mean_rate)
bandung_data$quad_sig[i] <- "High-High"
else if (bandung_data$JumlahKasus[i] < mean_rate & lag.listw(lw, bandung_data$JumlahKasus, zero.policy=TRUE)[i] < mean_rate)
bandung_data$quad_sig[i] <- "Low-Low"
else if (bandung_data$JumlahKasus[i] > mean_rate & lag.listw(lw, bandung_data$JumlahKasus, zero.policy=TRUE)[i] < mean_rate)
bandung_data$quad_sig[i] <- "High-Low"
else if (bandung_data$JumlahKasus[i] < mean_rate & lag.listw(lw, bandung_data$JumlahKasus, zero.policy=TRUE)[i] > mean_rate)
bandung_data$quad_sig[i] <- "Low-High"
} else {
bandung_data$quad_sig[i] <- "Not Sig"
}
}
# Visualisasi peta LISA
ggplot(bandung_data) +
geom_sf(aes(fill = quad_sig), color="white", size=0.2) +
scale_fill_manual(values=c("High-High"="red","Low-Low"="blue",
"High-Low"="orange","Low-High"="lightblue",
"Not Sig"="grey80")) +
labs(title="Peta Cluster LISA (Local Moran’s I)") +
theme_minimal()
# Buat vektor nama untuk tiap kategori (jika kolom quad_sig ada)
hh <- bandung_data$WADMKC[bandung_data$quad_sig == "High-High"]
ll <- bandung_data$WADMKC[bandung_data$quad_sig == "Low-Low"]
hl <- bandung_data$WADMKC[bandung_data$quad_sig == "High-Low"]
lh <- bandung_data$WADMKC[bandung_data$quad_sig == "Low-High"]
# Tampilkan ringkasan jumlah tiap kategori dan nama kecamatan yang termasuk
cat("Ringkasan jumlah kategori LISA:\n")## Ringkasan jumlah kategori LISA:##
## High-High Low-Low Not Sig
## 3 3 24##
## Nama kecamatan per kategori (jika ada):##
## High-High (Hotspot):## [1] "Antapani" "Buahbatu" "Rancasari"##
## Low-Low (Coldspot):## [1] "Astana Anyar" "Regol" "Sumur Bandung"##
## High-Low (Outlier - High surrounded by Low):## Tidak ada kecamatan High-Low yang signifikan.##
## Low-High (Outlier - Low surrounded by High):## Tidak ada kecamatan Low-High yang signifikan.4. Getis–Ord \(G_i^*\)
gstar <- localG(bandung_data$JumlahKasus, lw, zero.policy = TRUE)
bandung_data$zGstar <- as.numeric(gstar)
# Visualisasi peta Hotspot
ggplot(bandung_data) +
geom_sf(aes(fill = zGstar), color="white", size=0.2) +
scale_fill_gradient2(low="blue", mid="white", high="red", midpoint=0,
name="Gi* Z-score") +
labs(title="Peta Hotspot Getis–Ord Gi*") +
theme_minimal()