library(sf)
library(sp)
library(spdep)
library(dplyr)
library(ggplot2)
library(viridis)
library(stringr)

A - Teori dan Rumus

1. Autokorelasi spasial positif vs negatif

  • Autokorelasi spasial positif: kondisi ketika nilai suatu variabel pada satu lokasi cenderung mirip dengan nilai pada lokasi-lokasi di sekitarnya. Misalnya, daerah dengan pendapatan tinggi biasanya berdekatan dengan daerah lain berpendapatan tinggi.
  • Autokorelasi spasial negatif: kondisi ketika nilai suatu variabel pada suatu lokasi berlawanan dengan nilai tetangganya. Contohnya, area dengan kepadatan tinggi dikelilingi area dengan kepadatan rendah, membentuk pola papan catur.

2. Contoh fenomena nyata

  • Positif: Tingkat kejadian demam berdarah di suatu kecamatan cenderung tinggi jika kecamatan sekitarnya juga tinggi.
  • Negatif: Harga tanah di batas kota-desa; kota mahal, desa sekitarnya rendah.

3. Rumus matematis

Moran’s I (global)

\[ I = \frac{N}{S_0} \cdot \frac{\sum_{i=1}^N \sum_{j=1}^N w_{ij}(x_i - \bar{x})(x_j - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2}, \quad S_0 = \sum_i \sum_j w_{ij} \]

Geary’s C (global)

\[ C = \frac{(N-1)}{2S_0} \cdot \frac{\sum_{i=1}^N \sum_{j=1}^N w_{ij}(x_i - x_j)^2}{\sum_{i=1}^N (x_i - \bar{x})^2} \]

Local Moran’s I (LISA)

\[ I_i = \frac{(x_i - \bar{x})}{m_2} \sum_{j=1}^N w_{ij}(x_j - \bar{x}), \quad m_2 = \frac{1}{N}\sum_{k=1}^N (x_k - \bar{x})^2 \]

Getis–Ord \(G_i\) dan \(G_i^*\)

\[ G_i(d) = \frac{\sum_{j \neq i} w_{ij}(d) x_j}{\sum_{j \neq i} x_j}, \qquad G_i^*(d) = \frac{\sum_{j=1}^N w_{ij}(d) x_j}{\sum_{j=1}^N x_j}, \quad w_{ii}(d) = 1 \]

4. Perbedaan ukuran global vs lokal

  • Ukuran global (Moran’s I, Geary’s C): ringkasan untuk keseluruhan wilayah.
  • Ukuran lokal (Local Moran’s I, Getis Ord): fokus pada lokasi tertentu.

B - Simulasi data

1. Memuat Data Spasial dan Persiapan

Indo_Kec <- readRDS("gadm36_IDN_3_sp.rds")
Bandung_sf <- st_as_sf(Indo_Kec) %>%
  filter(NAME_2 == "Kota Bandung") %>%
  mutate(NAME_3_clean = str_trim(NAME_3) %>% str_to_lower())

2. Simulasi Data Kasus DBD dengan Hotspot

set.seed(0078)
n <- nrow(Bandung_sf)
kec_list <- Bandung_sf$NAME_3
pop <- sample(20000:90000, size = n, replace = TRUE)
centroids <- st_centroid(Bandung_sf)
coords <- st_coordinates(centroids)

center1 <- coords[round(n * 0.25), ] 
center2 <- coords[round(n * 0.7), ] 

# jarak ke pusat hotspot
dist1 <- sqrt((coords[, 1] - center1[1])^2 + (coords[, 2] - center1[2])^2)
dist2 <- sqrt((coords[, 1] - center2[1])^2 + (coords[, 2] - center2[2])^2)

# baseline & hotspot
lambda_base <- 10 
hot1 <- 30
hot2 <- 25

lambda_per10k <- lambda_base +
  hot1 * exp(-(dist1^2) / (2 * 0.04^2)) +
  hot2 * exp(-(dist2^2) / (2 * 0.06^2))

expected_cases <- lambda_per10k * (pop / 10000)
cases <- rpois(n = n, lambda = expected_cases)

DBD_dummy <- data.frame(
  Kecamatan = kec_list,
  Kecamatan_clean = str_trim(kec_list) %>% str_to_lower(),
  Population = pop,
  lambda_per10k = lambda_per10k,
  expected_cases = expected_cases,
  DBD_cases = cases,
  rate_per10k = cases / pop * 10000,
  stringsAsFactors = FALSE
)

3. Gabungkan Data dan Buat Peta Choropleth

Bandung_joined <- left_join(Bandung_sf, DBD_dummy, by = c("NAME_3_clean" = "Kecamatan_clean"))

p_counts <- ggplot(Bandung_joined) +
  geom_sf(aes(fill = DBD_cases), color = "grey30") +
  scale_fill_viridis_c(option = "plasma", name = "Kasus DBD", na.value = "grey90") +
  ggtitle("Choropleth: Jumlah Kasus DBD - Kota Bandung (Simulasi)") +
  theme_minimal()

p_rate <- ggplot(Bandung_joined) +
  geom_sf(aes(fill = rate_per10k), color = "grey30") +
  scale_fill_viridis_c(option = "plasma", name = "Rate per 10k", na.value = "grey90") +
  ggtitle("Choropleth: Rate Kasus DBD per 10.000 - Kota Bandung (Simulasi)") +
  theme_minimal()

print(p_counts)

print(p_rate)

Bandung_joined %>%
  st_drop_geometry() %>%
  select(NAME_3, Population, DBD_cases, rate_per10k) %>%
  arrange(desc(DBD_cases)) %>%
  head(10)
##              NAME_3 Population DBD_cases rate_per10k
## 1     Bandung Wetan      84190       537    63.78430
## 2       Batununggal      82383       452    54.86569
## 3     Sumur Bandung      62117       402    64.71658
## 4      Kiaracondong      70479       352    49.94395
## 5             Andir      74733       329    44.02339
## 6  Cibeunying Kidul      56187       311    55.35088
## 7         Rancasari      77076       308    39.96056
## 8       Mandalajati      71855       307    42.72493
## 9          Buahbatu      70505       298    42.26651
## 10 Cibeunying Kaler      46997       276    58.72715

4. Interpretasi

Berdasarkan visualisasi pada kedua peta, pola spasial sebaran kasus DBD di Kota Bandung secara jelas menunjukkan adanya pengelompokan (klaster). Terlihat sebuah konsentrasi kasus yang tinggi atau hotspot di wilayah tengah hingga timur kota. Pola ini sangat konsisten, baik pada peta yang menampilkan jumlah kasus absolut maupun pada peta tingkat risiko (SIR), yang menandakan bahwa area tersebut memang memiliki risiko penyebaran yang lebih tinggi dan bukan hanya karena populasi yang padat. Sebaliknya, wilayah di bagian pinggir kota seperti area barat, utara, dan selatan secara umum menunjukkan jumlah kasus dan tingkat risiko yang lebih rendah. Pola yang tidak acak ini mengindikasikan adanya faktor-faktor spasial tertentu yang berkontribusi terhadap tingginya kasus di area hotspot tersebut.

C - Pengukuran Autokorelasi

1. Moran’s I

Bandung_sp <- as(Bandung_joined, "Spatial")
nb <- poly2nb(Bandung_sp, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

mi <- moran.test(Bandung_joined$rate_per10k, lw, randomisation = TRUE)
mi
## 
##  Moran I test under randomisation
## 
## data:  Bandung_joined$rate_per10k  
## weights: lw    
## 
## Moran I statistic standard deviate = 5.6635, p-value = 7.417e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.61681121       -0.03448276        0.01322483

2. Geary’s C

gc <- geary.test(Bandung_joined$rate_per10k, lw, randomisation = TRUE)
gc
## 
##  Geary C test under randomisation
## 
## data:  Bandung_joined$rate_per10k 
## weights: lw   
## 
## Geary C statistic standard deviate = 5.4992, p-value = 1.907e-08
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.35759621        1.00000000        0.01364629

3. Local Moran’s I (LISA)

lisa <- localmoran(scale(Bandung_joined$rate_per10k)[,1], lw, zero.policy = TRUE)
Bandung_joined$LISA_I <- lisa[,1]
Bandung_joined$LISA_p <- lisa[,5]

z <- scale(Bandung_joined$rate_per10k)[,1]
lagz <- lag.listw(lw, z)
Bandung_joined$quad <- case_when(
  z >= 0 & lagz >= 0 & Bandung_joined$LISA_p <= 0.05 ~ "High-High",
  z < 0 & lagz < 0 & Bandung_joined$LISA_p <= 0.05 ~ "Low-Low",
  z >= 0 & lagz < 0 & Bandung_joined$LISA_p <= 0.05 ~ "High-Low",
  z < 0 & lagz >= 0 & Bandung_joined$LISA_p <= 0.05 ~ "Low-High",
  TRUE ~ "Not sig"
)

ggplot(Bandung_joined) +
  geom_sf(aes(fill = quad), color = "grey30") +
  scale_fill_manual(values = c("High-High"="#d73027","Low-Low"="#4575b4",
                               "High-Low"="#fdae61","Low-High"="#74add1",
                               "Not sig"="grey80")) +
  ggtitle("Peta Local Moran's I (LISA) - Kota Bandung") +
  theme_minimal()

4. Getis–Ord Gi*

Gi <- localG(Bandung_joined$rate_per10k, listw = lw)
Bandung_joined$GiZ <- as.numeric(Gi)

Bandung_joined$Gi_cat <- case_when(
  Bandung_joined$GiZ >= 1.96 ~ "Hotspot",
  Bandung_joined$GiZ <= -1.96 ~ "Coldspot",
  TRUE ~ "Not sig"
)

ggplot(Bandung_joined) +
  geom_sf(aes(fill = Gi_cat), color = "grey30") +
  scale_fill_manual(values = c("Hotspot"="#b2182b",
                               "Coldspot"="#2166ac",
                               "Not sig"="grey80")) +
  ggtitle("Peta Getis–Ord Gi* - Hotspot & Coldspot") +
  theme_minimal()

D - Diskusi Kritis

1. Pemanfaatan Hasil Analisis Autokorelasi Spasial

Hasil analisis autokorelasi spasial dapat memberikan informasi yang sangat berguna bagi Dinas Kesehatan Kota Bandung dalam menyusun strategi pencegahan dan intervensi penyakit menular seperti DBD. Dengan mengetahui lokasi hotspot, dinas kesehatan bisa lebih tepat sasaran dalam melakukan intervensi berbasis wilayah, misalnya penyemprotan fogging, edukasi masyarakat, perbaikan sanitasi, serta peningkatan layanan kesehatan di daerah rawan. Sementara itu, identifikasi coldspot membantu memastikan bahwa wilayah dengan risiko rendah tetap dipertahankan melalui monitoring rutin. Secara umum, analisis spasial membantu mengalokasikan sumber daya yang terbatas secara lebih efektif dan efisien, dengan prioritas pada kecamatan yang terbukti memiliki konsentrasi kasus tinggi.

2. Keterbatasan Analisis Autokorelasi Spasial

  • MAUP (Modifiable Areal Unit Problem): hasil analisis bisa berubah jika unit wilayah analisis diubah (misalnya menggunakan kelurahan alih-alih kecamatan). Dengan kata lain, pola yang terlihat bisa dipengaruhi oleh bagaimana batas wilayah didefinisikan.
  • Ukuran bobot spasial: pilihan metode bobot spasial (rook, queen, atau k-nearest neighbors) dapat menghasilkan hasil yang berbeda. Queen biasanya lebih inklusif karena mempertimbangkan semua sisi dan sudut, sementara rook hanya mempertimbangkan sisi. Pemilihan bobot spasial yang tidak tepat bisa mengubah interpretasi.
  • Masalah multiple testing pada analisis lokal: pada Local Moran’s I atau Getis–Ord Gi*, setiap wilayah diuji secara terpisah. Hal ini meningkatkan risiko kesalahan false positive, yaitu beberapa wilayah terlihat signifikan padahal sebenarnya tidak. Oleh karena itu, hasil lokal sebaiknya diinterpretasikan dengan hati-hati dan bila memungkinkan disertai dengan koreksi untuk multiple testing.