Konsep Teori

1. Penjelasan Autokorelasi Spasial Positif dan Negatif

Autokorelasi spasial positif terjadi ketika lokasi yang berdekatan memiliki nilai yang mirip atau cenderung serupa. Artinya, ada kecenderungan suatu variabel berkumpul dengan nilai yang sama (misalnya daerah dengan tingkat pendapatan tinggi dikelilingi oleh daerah dengan pendapatan tinggi juga).
Autokorelasi spasial negatif terjadi ketika lokasi yang berdekatan memiliki nilai yang berbeda secara signifikan. Artinya, ada kecenderungan suatu variabel tersebar dengan nilai berlawanan (misalnya daerah dengan angka kriminalitas tinggi berdekatan dengan daerah kriminalitas rendah).

2. Contoh Fenomena

Autokorelasi spasial positif: Penyebaran penyakit menular, misalnya kasus malaria sering terkonsentrasi di wilayah tertentu yang berdekatan.
Autokorelasi spasial negatif: Harga tanah di kota besar, di mana daerah pusat bisnis (harga tinggi) biasanya berdekatan dengan daerah pinggiran (harga rendah).

3. Rumus Matematis

Moran’s I
\[ I = \frac{n}{W} \cdot \frac{\sum_i \sum_j w_{ij} (x_i - \bar{x})(x_j - \bar{x})}{\sum_i (x_i - \bar{x})^2} \]
Mengukur autokorelasi spasial secara global dengan membandingkan kovariansi spasial terhadap varians total.
Geary’s C
\[ C = \frac{(n-1)}{2W} \cdot \frac{\sum_i \sum_j w_{ij} (x_i - x_j)^2}{\sum_i (x_i - \bar{x})^2} \]
Menekankan perbedaan lokal antara nilai tetangga; nilai kecil menunjukkan kesamaan tinggi, nilai besar menunjukkan perbedaan besar.
Local Moran’s \(I_i\)
\[ I_i = \frac{(x_i - \bar{x})}{m_2} \sum_j w_{ij}(x_j - \bar{x}) \]
Digunakan untuk mengidentifikasi klaster atau outlier pada lokasi tertentu dengan melihat hubungan nilai sebuah titik dengan tetangganya.

dengan
\[ m_2 = \frac{\sum_i (x_i - \bar{x})^2}{n} \]

Getis–Ord \(G_i\) dan \(G_i^*\)
\[ G_i = \frac{\sum_j w_{ij} x_j}{\sum_j x_j}, \quad G_i^* = \frac{\sum_j w_{ij} x_j}{\sum_j x_j} \]
Digunakan untuk mengukur konsentrasi nilai tinggi (hot spot) atau rendah (cold spot) di sekitar lokasi \(i\).

Perbedaan utama adalah pada definisi bobot \(w_{ij}\), apakah titik \(i\) sendiri ikut dihitung atau tidak.

4. Perbedaan Ukuran Global dan Lokal

Ukuran global (Moran’s I dan Geary’s C) memberikan gambaran umum pola spasial untuk seluruh wilayah studi. Cocok untuk mengetahui apakah secara keseluruhan terdapat autokorelasi spasial.
Ukuran lokal (Local Moran’s I dan Getis–Ord Gi) digunakan untuk mendeteksi adanya klaster atau outlier lokal. Dengan ukuran lokal, kita bisa mengetahui lokasi spesifik yang memiliki pengelompokan nilai tinggi atau rendah.

Analisis Data (Simulasi)

1. Data Spasial Kota Bandung

library(sf)

## Warning: package 'sf' was built under R version 4.4.3

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

library(sp)

## Warning: package 'sp' was built under R version 4.4.3

library(spdep)

## Warning: package 'spdep' was built under R version 4.4.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.4.3

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

## Warning: package 'dplyr' was built under R version 4.4.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

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.3

library(mapview)

## Warning: package 'mapview' was built under R version 4.4.3

library(spData)

# Arahkan ke folder kerja
setwd("C:/Kuliah/Semester 5/Spasial")

# Baca data RDS level kecamatan Indonesia
Indo_Kec <- readRDS("gadm36_IDN_3_sp.rds")

# Filter hanya Kota Bandung
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung", ]
Bandung$id <- seq_len(nrow(Bandung))
Bandung_sf <- sf::st_as_sf(Bandung)

# Jika data asli tidak ada, siapkan grid simulasi (30 poligon)
if(!exists("Bandung_sf") || nrow(Bandung_sf) == 0){
  bb <- sf::st_as_sfc(sf::st_bbox(c(xmin=107.55, ymin=-6.98,
                                    xmax=107.72, ymax=-6.85), crs=4326))
  grid <- sf::st_make_grid(bb, n=c(6,5)) |> sf::st_as_sf() |>
          dplyr::mutate(id = dplyr::row_number())
  Bandung_sf <- grid
}

2. Simulasi Kasus Penyakit Menular

# Bangun matriks tetangga (queen)
nb <- spdep::poly2nb(sf::as_Spatial(Bandung_sf), queen = TRUE)

# Listw untuk analitik global/lokal
lwW <- spdep::nb2listw(nb, style="W") # row-standardized weights
lwB <- spdep::nb2listw(nb, style="B") # binary weights (untuk raw G)

# ==== Simulasi variabel kovariat spasial ====
set.seed(123)   # FIX agar hasil tidak berubah
N <- nrow(Bandung_sf)
z <- rnorm(N)
lz <- spdep::lag.listw(lwW, z)

# Model rata-rata kasus (log-linear)
beta0 <- 1.2; beta1 <- 0.6; beta2 <- 1.0
lambda <- exp(beta0 + beta1*scale(z) + beta2*scale(lz))

# Simulasi populasi dan kasus
pop <- round(runif(N, 3000, 8000))       # penduduk tiap unit
cases <- rpois(N, lambda = lambda * (pop/10000))  # kasus Poisson
Diare <- dplyr::tibble(id = Bandung_sf$id, pop = pop, cases = cases,
                       rate10k = (cases/pop)*10000)

# Gabungkan ke shapefile
Bandung_merged <- dplyr::left_join(Bandung_sf, Diare, by="id")

3. Peta Choropleth

ggplot(Bandung_merged) +
  geom_sf(aes(fill = rate10k), color="white", size=0.2) +
  scale_fill_viridis_c(option="magma") +
  labs(title="Rate Diare (per 10.000 penduduk) — Simulasi",
       fill="Rate") +
  theme_minimal()

4. Pola Spasial yang Terlihat Secara Visual

Berdasarkan peta choropleth di atas, terlihat adanya distribusi spasial yang tidak merata pada rate diare per 10.000 penduduk. Kecamatan di bagian utara cenderung membentuk klaster dengan rate tinggi (hotspot), yang ditunjukkan dengan warna kuning hingga merah. Sebaliknya, sebagian kecamatan di bagian selatan dan timur tampak memiliki rate rendah (coldspot), dengan warna ungu hingga hitam. Sementara itu, kecamatan lain berada pada level menengah dan tersebar di sekitar wilayah tengah. Pola ini menunjukkan adanya indikasi autokorelasi spasial positif, yaitu nilai tinggi cenderung berdekatan dengan nilai tinggi, dan nilai rendah cenderung berdekatan dengan nilai rendah.

Pengukuran Autokorelasi

Moran’s I

moran_res <- spdep::moran.test(Bandung_merged$rate10k, lwW,
randomisation = TRUE, alternative = "two.sided")
moran_res

## 
##  Moran I test under randomisation
## 
## data:  Bandung_merged$rate10k  
## weights: lwW    
## 
## Moran I statistic standard deviate = 4.8496, p-value = 1.237e-06
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.50481782       -0.03448276        0.01236637

Nilai Moran’s I:
Hasil uji memberikan nilai Moran’s I = 0.5048, lebih besar dari ekspektasi acaknya (–0.0345).

Signifikansi Statistik (uji randomisasi):
Nilai p-value = 1.237 × 10⁻⁶, jauh lebih kecil daripada 0.05. Hal ini menunjukkan pola spasial yang terdeteksi sangat signifikan secara statistik.

Interpretasi:
Nilai Moran’s I yang positif dan signifikan menandakan adanya autokorelasi spasial positif. Artinya, kecamatan dengan angka kejadian diare yang tinggi cenderung berdekatan dengan kecamatan lain yang juga memiliki angka tinggi, begitu pula kecamatan dengan angka rendah cenderung berkelompok dengan kecamatan lain yang rendah. Dengan kata lain, terdapat pola spasial berupa cluster hotspot (tinggi–tinggi) dan cluster coldspot (rendah–rendah) di wilayah studi.

Geary’s C

geary_res <- spdep::geary.test(Bandung_merged$rate10k, lwW,
randomisation = TRUE, alternative = "two.sided")
geary_res

## 
##  Geary C test under randomisation
## 
## data:  Bandung_merged$rate10k 
## weights: lwW   
## 
## Geary C statistic standard deviate = 4.3461, p-value = 1.386e-05
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.47515619        1.00000000        0.01458327

Perbandingan dengan Moran’s I:
Hasil uji memberikan nilai Geary’s C = 0.475 dengan p-value = 1.39 × 10⁻⁵. Nilai ini jauh lebih kecil daripada ekspektasi acaknya (1.0), sehingga menunjukkan adanya autokorelasi spasial positif. Temuan ini konsisten dengan hasil Moran’s I (0.505, signifikan), sehingga kedua ukuran sama-sama mendeteksi adanya clustering kasus diare di Kota Bandung.

Perbedaan sensitivitas Moran’s I vs Geary’s C:
- Moran’s I lebih menekankan pada hubungan global antara nilai suatu lokasi dengan rata-rata tetangganya. Ia peka terhadap pola umum di seluruh wilayah.
- Geary’s C lebih menekankan pada perbedaan lokal antar pasangan wilayah yang bersebelahan, sehingga lebih sensitif terhadap ketidakseragaman pada skala lokal. Dengan kata lain, jika terdapat perbedaan tajam antar kecamatan yang bertetangga, Geary’s C lebih cepat mendeteksinya dibanding Moran’s I.

Local Moran’s I (LISA)

x <- scale(Bandung_merged$rate10k)[,1]
lagx <- spdep::lag.listw(lwW, x)

# Hitung Local Moran's I
lisa <- spdep::localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")

alpha <- 0.05
quad <- dplyr::case_when(
  x >= 0 & lagx >= 0 ~ "High-High",
  x < 0 & lagx < 0 ~ "Low-Low",
  x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
  x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
)

# Gabungkan hasil dengan data spasial
Bandung_LISA <- dplyr::bind_cols(Bandung_merged, lisa_df) |>
  dplyr::mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))

# Ambil tabel tanpa geometry, hanya kecamatan yang signifikan
hasil_LISA_signif <- Bandung_LISA |>
  dplyr::select(NAME_3, quad, Ii, Pi.two.sided) |>
  sf::st_set_geometry(NULL) |>
  dplyr::filter(Pi.two.sided <= 0.05)

hasil_LISA_signif

Berdasarkan hasil perhitungan Local Moran’s I pada tingkat signifikansi 5%:

High-High (Hotspot): Teridentifikasi pada Kecamatan Cidadap, Coblong, Sukajadi, dan Sukasari. Artinya, keempat kecamatan ini memiliki angka kasus diare tinggi dan dikelilingi oleh tetangga dengan angka tinggi pula. Hal ini menunjukkan adanya klaster konsentrasi kasus tinggi di wilayah barat–utara Kota Bandung.
Low-Low (Coldspot): Tidak ada kecamatan yang masuk kategori ini pada taraf signifikansi 5%.
Low-High (Outlier): Tidak ada kecamatan dengan nilai rendah yang dikelilingi tetangga bernilai tinggi yang signifikan.
High-Low (Outlier): Tidak ditemukan kecamatan dengan angka tinggi yang justru dikelilingi tetangga rendah.

ggplot(Bandung_LISA) +
geom_sf(aes(fill = quad), color="white", size=0.2) +
scale_fill_manual(values=c(
"High-High"="#d73027","Low-Low"="#4575b4",
"High-Low (Outlier)"="#fdae61","Low-High (Outlier)"="#74add1","Not significant"="grey85"
)) +
labs(title="Local Moran's I (LISA)", fill="Kategori") +
theme_minimal()

Pola spasial signifikan:
- Terlihat adanya klaster High-High (hotspot) di wilayah utara Kota Bandung (Cidadap, Coblong, Sukajadi, Sukasari). Kecamatan-kecamatan ini menunjukkan angka kasus diare tinggi yang berdekatan secara spasial, sehingga berpotensi menjadi pusat konsentrasi penularan penyakit.
Wilayah lain:
- Sebagian besar kecamatan lainnya berwarna abu-abu (Not significant), artinya tidak menunjukkan pola autokorelasi lokal yang signifikan. Dengan kata lain, tidak ditemukan Low-Low (coldspot) maupun outlier Low-High atau High-Low pada tingkat signifikansi 5%.

Getis-Ord \(G_i^{\ast}\)

x_raw <- Bandung_merged$rate10k
sum_x <- sum(x_raw)

# Matriks bobot biner (tidak distandarisasi)
Wb <- spdep::listw2mat(lwB) # diag = 0

# G_i (tanpa i)
num_G <- as.numeric(Wb %*% x_raw) 
den_G <- (sum_x - x_raw) 
G_raw <- num_G / den_G

# G_i^* (dengan i)
Wb_star <- Wb; diag(Wb_star) <- 1
num_Gs <- as.numeric(Wb_star %*% x_raw) 
den_Gs <- sum_x 
G_star_raw <- num_Gs / den_Gs

# Z-skor (spdep::localG)
Gz <- spdep::localG(x_raw, listw = lwW) 

Bandung_G <- dplyr::mutate(Bandung_merged,
  G_raw = G_raw,
  G_star_raw = G_star_raw,
  z_Gistar = as.numeric(Gz),
  hotcold = dplyr::case_when(
    z_Gistar >= 1.96 ~ "Hot spot (p<0.05)",
    z_Gistar <= -1.96 ~ "Cold spot (p<0.05)",
    TRUE ~ "Not significant"
  )
)

# Ringkasan hasil
summary(dplyr::select(Bandung_G, G_raw, G_star_raw, z_Gistar))

##      G_raw           G_star_raw         z_Gistar                geometry 
##  Min.   :0.00000   Min.   :0.01069   Min.   :-1.7571   MULTIPOLYGON :30  
##  1st Qu.:0.05115   1st Qu.:0.07107   1st Qu.:-1.1364   epsg:NA      : 0  
##  Median :0.09574   Median :0.11959   Median :-0.9334   +proj=long...: 0  
##  Mean   :0.13832   Mean   :0.16419   Mean   :-0.2353                     
##  3rd Qu.:0.16788   3rd Qu.:0.18923   3rd Qu.: 0.1827                     
##  Max.   :0.48942   Max.   :0.53346   Max.   : 4.1541

hasil_Gi_star <- Bandung_G |>
  dplyr::select(NAME_3, z_Gistar, hotcold) |>
  sf::st_set_geometry(NULL) |>
  dplyr::filter(hotcold != "Not significant") |>
  dplyr::arrange(desc(z_Gistar))

hasil_Gi_star

Hasil Identifikasi Hot Spot dan Cold Spot (Gi*, p < 0.05)

Hot spot (z > 1.96):
Terdeteksi pada empat kecamatan, yaitu:
- Cidadap (z = 4.15)
- Sukajadi (z = 3.47)
- Sukasari (z = 2.76)
- Coblong (z = 2.09)
Keempat kecamatan ini membentuk klaster konsentrasi kasus diare tinggi di wilayah utara–barat Kota Bandung, yang dapat menjadi fokus utama intervensi kesehatan.
Cold spot (z < -1.96):
Tidak ditemukan kecamatan dengan z-score signifikan negatif pada taraf 5%.
Tidak signifikan:
Sebagian besar kecamatan lainnya tidak menunjukkan pola hot spot maupun cold spot secara statistik.

ggplot(Bandung_G) +
  geom_sf(aes(fill = G_star_raw), color="white", size=0.2) +
  scale_fill_viridis_c() +
  labs(title="Raw Getis–Ord Gi* (proporsi massa tetangga)", 
       fill="G*_raw") +
  theme_minimal()

Berdasarkan peta Raw Getis–Ord Gi* yang ditampilkan:

Nilai Gi* yang tinggi (ditunjukkan dengan warna kuning–hijau) terlihat terkonsentrasi di wilayah utara Kota Bandung. Hal ini mengindikasikan adanya potensi hotspot, yakni kecamatan dengan tingkat kasus diare relatif tinggi yang berdekatan satu sama lain.

Nilai Gi* yang rendah (ditunjukkan dengan warna biru tua–ungu) banyak muncul di wilayah selatan, yang mencerminkan kecenderungan ke arah coldspot atau daerah dengan tingkat kasus rendah.
Pola spasial ini menunjukkan adanya heterogenitas antar kecamatan, di mana bagian utara menjadi pusat intensitas kasus, sedangkan bagian selatan relatif lebih rendah.

ggplot(Bandung_G) +
  geom_sf(aes(fill = hotcold), color="white", size=0.2) +
  scale_fill_manual(values=c("Hot spot (p<0.05)"="#b2182b",
                             "Cold spot (p<0.05)"="#2166ac",
                             "Not significant"="grey85")) +
  labs(title="Getis–Ord Gi* — Hot/Cold Spots (z-skor)", fill=NULL) +
  theme_minimal()

Berdasarkan peta Getis–Ord Gi* (z-skor) yang ditampilkan:

Hot spot signifikan (p < 0.05): Terlihat jelas di wilayah utara Kota Bandung, meliputi kecamatan seperti Cidadap, Sukajadi, Sukasari, dan Coblong. Kawasan ini menunjukkan konsentrasi kasus diare tinggi yang signifikan secara statistik.
Tidak signifikan (abu-abu): Sebagian besar kecamatan lainnya tidak masuk kategori hot spot maupun cold spot, sehingga tidak teridentifikasi sebagai klaster signifikan.
Cold spot signifikan: Tidak ditemukan pada peta ini, berbeda dengan hasil LISA yang sempat mengindikasikan adanya pola Low-Low di bagian tertentu.

Diskusi Kritis

Hasil analisis autokorelasi spasial (misalnya Getis–Ord Gi*) dapat membantu Dinas Kesehatan dalam:

- Identifikasi klaster signifikan: Hot spot menunjukkan kecamatan dengan tingkat kasus penyakit menular yang tinggi dan signifikan secara statistik, sedangkan cold spot menunjukkan wilayah dengan tingkat kasus rendah.

- Prioritas intervensi: Sumber daya kesehatan (misalnya tenaga medis, vaksinasi, sosialisasi) dapat difokuskan terlebih dahulu ke wilayah hot spot di utara Bandung.

- Perencanaan pencegahan: Wilayah cold spot perlu tetap dijaga agar tidak terjadi peningkatan kasus, misalnya dengan monitoring berkala dan pencegahan berbasis komunitas.

- Evidensi berbasis data spasial: Hasil peta memberikan dasar yang lebih kuat dalam pengambilan keputusan berbasis lokasi (place-based policy).

Meskipun bermanfaat, analisis autokorelasi spasial memiliki sejumlah keterbatasan, di antaranya:

MAUP (Modifiable Areal Unit Problem)
Hasil analisis sangat dipengaruhi oleh skala dan cara pengelompokan wilayah administrasi. Misalnya, jika batas kecamatan berubah, hasil hot/cold spot bisa berbeda.
Ukuran bobot spasial
Pilihan metode tetangga (rook, queen, atau k-nearest neighbors) dapat menghasilkan hasil yang berbeda. Pemilihan skema bobot perlu disesuaikan dengan konteks epidemiologi dan pola penyebaran penyakit.
Masalah multiple testing
Analisis lokal seperti LISA atau Gi* melakukan uji statistik pada banyak lokasi sekaligus. Hal ini meningkatkan risiko false positive (wilayah dianggap signifikan padahal tidak). Koreksi seperti FDR (False Discovery Rate) kadang diperlukan.

Kesimpulan

Analisis autokorelasi spasial penting untuk memetakan distribusi penyakit menular di Kota Bandung. Namun, hasilnya perlu ditafsirkan dengan hati-hati karena ada keterbatasan metodologis seperti MAUP, pemilihan bobot spasial, dan isu multiple testing.

Tugas 1 - Spatial Statistics

Ahmad Firmansyah