Bagian A : Konsep Teori

1. Penjelasan singkat

Autokorelasi spasial positif berarti lokasi-lokasi yang berdekatan cenderung memiliki nilai variabel yang mirip (mis. daerah tetangga sama-sama tinggi atau sama-sama rendah). Dalam kata lain, ada klaster homogenitas ruang: nilai tinggi menempel pada nilai tinggi, nilai rendah menempel pada nilai rendah.

Autokorelasi spasial negatif berarti lokasi-lokasi yang berdekatan cenderung memiliki nilai variabel yang berbeda (mis. wilayah dengan nilai tinggi dikelilingi oleh wilayah bernilai rendah). Ini menunjukkan pola disparitas atau heterogenitas ruang — nilai tinggi berdekatan dengan nilai rendah.

2. Contoh sederhana fenomena

  • Autokorelasi spasial positif: harga tanah di sebuah kota — kawasan pusat bisnis dan kawasan elite biasanya berdekatan sehingga menimbulkan klaster harga tinggi; daerah permukiman padat miskin juga membentuk klaster harga rendah.
  • Autokorelasi spasial negatif: pola vegetasi akibat fragmentasi habitat — area sangat vegetasi tinggi diapit oleh lahan terdegradasi sehingga nilai tinggi dikelilingi oleh rendah; atau pola rumah mewah yang menyebar di antara kawasan sederhana (jika terjadi segregasi spasial khusus).

3. Rumus matematis

Sebelum rumus, definisikan beberapa notasi umum:

  • \(n\): jumlah lokasi (unit spasial).
  • \(y_i\): nilai variabel di lokasi $i$.
  • \(\bar{y} = \frac{1}{n}\sum_{i=1}^n y_i\) : (rata-rata sampel)
  • \(w_{ij}\): elemen matriks bobot spasial \(W\) (menunjukkan relasi antara lokasi \(i\) dan \(j\)). Biasanya \(w_{ii}=0\).
  • \(S_0 = \sum_{i=1}n\sum_{j=1}n w_{ij}\) : (jumlah semua bobot)

Moran’s I

Moran’s I adalah ukuran autokorelasi spasial global. Rumusnya:

\[ I = \frac{n}{S_0}\cdot \frac{\sum_{i=1}^n\sum_{j=1}^n w_{ij}(y_i-\bar{y})(y_j-\bar{y})}{\sum_{i=1}^n (y_i-\bar{y})^2}. \]

Interpretasi: nilai \(I\) positif menunjukkan autokorelasi positif; nilai negatif menunjukkan autokorelasi negatif. Nilai ekspektasi di bawah hipotesis nol (tidak ada autokorelasi) biasanya

\[E[I] = -\frac{1}{n-1}. \]

Geary’s C

Geary’s C lebih sensitif terhadap perbedaan lokal daripada Moran. Rumusnya:

\[ C = \frac{(n-1)}{2S_0} \cdot \frac{\sum_{i=1}^n\sum_{j=1}^n w_{ij}(y_i-y_j)^2}{\sum_{i=1}^n (y_i-\bar{y})^2}. \]

Interpretasi: untuk Geary’s C, nilai \(C<1\) menandakan autokorelasi positif (karena perbedaan antar tetangga kecil), \(C>1\) menandakan autokorelasi negatif; nilai \(C\approx 1\) berarti tidak ada autokorelasi.

Local Moran’s $I_i$ (Anselin)

Local Moran memberikan ukuran autokorelasi pada tiap lokasi \(i\) (analisis lokal). Sering dituliskan sebagai:

\[ I_i = (y_i-\bar{y}) \sum_{j=1}^n w_{ij}(y_j-\bar{y}) / m_2, \]

dengan

\[ m_2 = \frac{\sum_{k=1}^n (y_k-\bar{y})^2}{n}. \]

Interpretasi: nilai \(I_i\) tinggi positif menunjukkan lokasi \(i\) adalah bagian dari klaster nilai tinggi (high-high); nilai \(I_i\) tinggi negatif menunjukkan klaster low-low atau outlier high yang dikelilingi low (interpretasi tergantung tanda dan signifikansi).

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

Getis–Ord menguji klaster “hotspot” (nilai tinggi) dan “coldspot” (nilai rendah).

\(G_i\) (tanpa self)

Salah satu bentuk statistik lokal Getis-Ord (tanpa memasukkan self-weight) adalah:

\[ G_i = \frac{\sum_{j=1}^n w_{ij} y_j}{\sum_{j=1}^n y_j}. \]

\(G_i^*\) (dengan self)

Bentuk yang umum dipakai adalah \(G_i^*\) yang memasukkan elemen diagonal (self):

\[ G_i^* = \frac{\sum_{j=1}^n w_{ij} y_j}{\sum_{j=1}^n y_j} \quad\text{often standardized to a } z\text{-score for inference.} \]

Namun dalam praktik statistik spasial, \(G_i^*\) sering dinormalisasi menjadi statistik \(z\) dengan bentuk umum:

\[ z(G_i^*) = \frac{\sum_j w_{ij} y_j - \mu_i}{\sigma_i}, \]

dengan \(\mu_i\) dan \(\sigma_i\) adalah ekspektasi dan standar deviasi dari jumlah tertimbang di bawah hipotesis nol (bervariasi menurut definisi \(W\) dan asumsi).

Catatan: ada beberapa variasi notasi untuk \(G_i\) dan \(G_i^*\) dalam literatur. Kunci penggunaannya: Getis–Ord menyorot klaster nilai tinggi (hotspot) dan nilai rendah (coldspot), sedangkan Local Moran menilai keterkaitan nilai suatu lokasi dengan tetangganya secara kuadrat (berkaitan tanda).

4. Perbedaan utama antara ukuran global dan lokal

  • Skala pengukuran: ukuran global (mis. Moran’s I, Geary’s C) merangkum seluruh pola autokorelasi dalam satu nilai ringkasan untuk seluruh kawasan studi. Ukuran lokal (mis. Local Moran \(I_i\), Getis–Ord \(G_i^*\)) menghitung statistik untuk tiap lokasi sehingga bisa mendeteksi di mana klaster atau outlier berada.

  • Deteksi pola lokal: global dapat menunjukkan ada autokorelasi secara umum, tetapi tidak menginformasikan lokasi klaster atau outlier. Lokal mengidentifikasi titik atau wilayah spesifik (hotspot, coldspot, high–low outliers, dll.).

  • Interpretasi: global bersifat agregat; nilai positif/negatif memberi gambaran keseluruhan. Lokal memerlukan koreksi multipel dan uji signifikansi spasial (mis. permutasi) karena banyak uji simultan.

  • Sensitivitas: beberapa statistik berbeda sensitivitasnya terhadap pola lokal vs global (Geary lebih sensitif terhadap perbedaan lokal; Moran menangkap struktur global lebih baik).

Bagian B : Analis Data (Simulasi)

Simulasi Data TBC Kota Bandung

Import Library dan Tentukan set.seed

library(sf)
## 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
library(sp)
## Warning: package 'sp' was built under R version 4.3.3
library(spdep)
## 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')`
library(dplyr)
## 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
library(ggplot2)
library(mapview)
library(spData)
set.seed(57)

Membuat Peta Spasial Kota Bandung

setwd("C:/Users/Johanes Credo/Downloads/Spatial")
Indo_Kec <- readRDS("gadm36_IDN_3_sp.rds")
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung", ]
Bandung$id <- seq_len(nrow(Bandung))
Bandung_sf <- sf::st_as_sf(Bandung)

if(!exists("Bandung_sf")) {
  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
}

Bandung_sf <- sf::st_make_valid(Bandung_sf)
Bandung_sf$id <- if(!"id" %in% names(Bandung_sf)) seq_len(nrow(Bandung_sf)) else Bandung_sf$id
ggplot(Bandung_sf) + geom_sf(fill="grey90", color="white") +
labs(title="Kota Bandung") + theme_minimal()

### Objek neighbor, bobot spasial, dan simulasi data

# Membuat objek neighbor
nb <- spdep::poly2nb(sf::as_Spatial(Bandung_sf), queen = TRUE)
# Membuat bobot spasial
lwW <- spdep::nb2listw(nb, style="W") # row-standardized
lwB <- spdep::nb2listw(nb, style="B") # binary (untuk raw G)
# Simulasi data
N <- nrow(Bandung_sf)
z <- rnorm(N)
lz <- spdep::lag.listw(lwW, z)

Membuat simulasi data dan visualisasi peta spasial

# Parameter model simulasi
beta0 <- 1.3
beta1 <- 0.5
beta2 <- 1.1

# Variabel acak untuk heterogenitas spasial
N <- nrow(Bandung_sf)
z  <- rnorm(N)
lz <- rnorm(N)

# Rata-rata (lambda) untuk kasus TBC
lambda <- exp(beta0 + beta1*scale(z) + beta2*scale(lz))

# Populasi simulasi (acak antara 3.000–8.000 per kecamatan)
pop <- round(runif(N, 3000, 8000))

# Kasus TBC per kecamatan (Poisson)
cases <- rpois(N, lambda = lambda * (pop/10000))

# Buat tabel data TBC
TBC <- dplyr::tibble(
  id     = Bandung_sf$id,
  pop    = pop,
  cases  = cases,
  rate10k = (cases/pop) * 10000
)

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

# Plot peta hasil simulasi
ggplot(Bandung_merged) +
  geom_sf(aes(fill = rate10k), color = "white", size = 0.2) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Rate TBC (per 10.000) – Simulasi",
       fill = "Rate") +
  theme_minimal()

Interpretasi: Dapat dilihat bahwa sebagian besar daerah terindikasi rendah dengan ditandai warna ungu tua. Untuk daerah dengan kasus simulasi tinggi berada di utara sebelah kiri dan daerah dengan kasus sedang beberapa tersebar di tengah peta. Dengan begitu cukup terlihat pola dimana bagian bawah didominasi oleh daerah kasus rendah, bagian tengah tersebar daerah dengan kasus sedang, dan di atas terdapat bagian peta yang terindikasi kasus tinggi. Hal ini mengindikasikan adanya pola autokorelasi spasial.

Bagian C : Perhitungan Autokorelasi

Global 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 = 1.7084, p-value = 0.08756
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.115255472      -0.034482759       0.007682036

Interpretasi: Distribusi rate TBC per 10.000 pada simulasi ini menunjukkan Moran’s I positif kecil (0.115) dengan p-value 0.0876, sehingga tidak signifikan pada taraf 5%, tetapi ada indikasi lemah adanya pola spasial (wilayah dengan angka mirip cenderung berdekatan).

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 = 1.9933, p-value = 0.04623
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.72026131        1.00000000        0.01969604

Interpretasi: Distribusi rate TBC per 10.000 pada simulasi ini menunjukkan lebih kecil dari 1 (0.7203) dengan p-value 0.04623, sehingga signifikan pada taraf 5%, sehingga ada bukti bahwa nilai rate TBC di wilayah yang bertetangga memang cenderung mirip.

Perbandingan Global Moran’s I dan Geary’s C Hasil Moran’s I menunjukkan bahwa terdapat indikasi lemah adanya hubungan autokorelasi spasial, tapi Geary’s C menunjukkan ada bukti bahwa rate TBC di wilayah yang bertetangga cenderung mirip. Hal ini menunjukkan bahwa Global Moran’s I lebih sensitif pada pola global, sedangkan Geary’s C lebih sensitif pada pola tetangga (perbedaan tetangga dekat)

Local Moran’s I

x <- scale(Bandung_merged$rate10k)[,1]
lagx <- spdep::lag.listw(lwW, x)
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)"
)
Bandung_LISA <- dplyr::bind_cols(Bandung_merged, lisa_df) |>
dplyr::mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))

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

Interpretasi: Hasil Local Moran’s I (LISA) menununjukkan hanya ada 1 daerah yang signifikan ditandai dengan warna biru muda (Low-High) dan untuk sisanya berwarna yang artinya tidak signifikan. Hal ini menunjukkan indikasi hubungan tetangga yang rendah. Untuk daera signifikan (biru muda), terindikasi Low-High artinya daerah dengan kasus rendah dikelilingi daerah kasus tinggi.

Getis-Ord

x_raw <- Bandung_merged$rate10k
sum_x <- sum(x_raw)
# Matriks bobot biner (tidak distandarisasi)
Wb <- spdep::listw2mat(lwB) # diag=0 secara default
# G_i (tanpa i)
num_G <- as.numeric(Wb %*% x_raw) # Σ_{ji} w_ij x_j
den_G <- (sum_x - x_raw) # Σ_{ji} x_j
G_raw <- num_G / den_G
# G_i^* (dengan i) — set w_ii = 1
Wb_star <- Wb; diag(Wb_star) <- 1
num_Gs <- as.numeric(Wb_star %*% x_raw) # Σ_j w_ij^* x_j
den_Gs <- sum_x # Σ_j x_j
G_star_raw <- num_Gs / den_Gs
# Z-skor (yang dihitung spdep::localG)
Gz <- spdep::localG(x_raw, listw = lwW) # z(G_i^*) dengan listw yang diberikan
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 (p0.05)",
z_Gistar <= -1.96 ~ "Cold spot (p0.05)",
TRUE ~ "Not significant"
))
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.03323   Min.   :-1.30861   MULTIPOLYGON :30  
##  1st Qu.:0.05254   1st Qu.:0.06900   1st Qu.:-0.87301   epsg:NA      : 0  
##  Median :0.10099   Median :0.11082   Median :-0.38842   +proj=long...: 0  
##  Mean   :0.13963   Mean   :0.16725   Mean   :-0.16243                     
##  3rd Qu.:0.21477   3rd Qu.:0.25612   3rd Qu.:-0.01006                     
##  Max.   :0.41915   Max.   :0.45633   Max.   : 3.36635

Interpretasi: - G_raw (Gertis-Ord G), artinya mengukur intensitas clustering nilai tinggi/rendah di sekitar setiap wilayah.

  • G_star_raw, artinya G yang mempertimbangkan bobot spasial termasuk dirinya sendiri.

  • z_Gistar, digunakan untuk uji signifikansi dengan z > 1.96 terindikasi hotspot, z < -1.96 terindikasi coldspot, dan -1.96 < z < 1.96 tidak signifikan.

Dari hasil di atas, nilai z berkisar di rentang -1.3 sampai -0.01 yang menunjukkan sebagian daerah tidak signifikan atau tidak terindikasi hotspot/coldspot. Namun, terdapat 1 nilai z = 3.36635 yang artinya terdapat 1 daerah yang terindikasi sebagai hotspot.

Visualisasi Getis-Ord

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

Intrpretasi: Gertis Ord G* digunakan untuk mengukur wilayah terindikasi hotspot atau coldspot dengan indikasi 0 - 0.1 (ungu tua) artinya tidak ada indikasi, 0.2 - 0.3 (hijau) artinya kontribusi sedang, dan 0.4 - 0.45 (kuning) artinya kontribusi tinggi atau ada indikasi hotspot atau coldspot.

Dari hasil di atas, didapatkan bahwa daerah bawah dan kanan didominasi oleh daerah dengan kontribusi rendah, daerah tengah didominasi oleh daerah dengan kontribusi sedang, dan daerah atas terdapat daerah dengan kontribusi tinggi. Hal ini menunjukkan bahwa daerah yang berpotensi untuk terindikasi sebagai hotspot/coldspot adalah bagian atas.

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

Interpretasi: Dari hasil di atas, terlihat terdapat satu daerah yang terindikasi blok merah (hotspot), dan sisanya menunjukkan warna abu artinya tidak signifikan (tidak ada indikasi). Daerah hotspot memiliki arti bahwa daerah tersebut dikelilingi oleh daerah yang memiliki nilai rate TBC yang tinggi. Hasil hotspot ini tidak serta-merta menyatakan bahwa daerah tersebut memiliki nilai rate TBC yang tinggi. Namun, daerah tersebut dikelilingi oleh daerah sekitarnya yang memiliki kasus tinggi sehingga ikut masuk.

Hal ini sedikit berbeda dengan hasil LISA yang menunjukkan daerah yang sama sebagai Low-High (daerah dengan kasus TBC rendah dikelilingi daerah kasus TBC tinggi). Perbedaan ini terjadi dikarenakan perbedaan fokus kedua metode. LISA lebih fokus dan sensitif pada outlier, sedangkan Gertis-Ord G* lebih fokus terhadap pendeteksian hotspot.

Bagian D : Diskusi Kritis

Bagaimana hasil analisis autokorelasi spasial bisa membantu dinas kesehatan dalam menyusun strategi pencegahan dan intervensi penyakit menular di Kota Bandung?

Hasil analisis autokorelasi spasial bisa membantu dinas kesehatan dalam menyusun strategi pencegahan dan intervensi penyakit menular di Kota Bandung. Dengan analisis autokorelasi spasial, dinas kesehatan dapat melihat wilayah mana yang merupakan wilayah kasus rendah/tinggi. Dan juga dapat dilihat lebih mengenai pola peta spasial apakah berpola saling mempengaruhi (autokorelasi) atau tidak, serta dapat mendeteksi wilayah outlier (High-High, Low-Low, High-Low/Low-High) dan juga hotspot atau coldspot.

Dengan dapat mendeteksi hal-hal tersebut, dinas kesehatan dapat melakukan pencegahan dan intervensi lebih mudah. Pencegahan dan intervensi dapat lebih difokuskan pada wilayah yang terdeteksi sebagai outlier maupun hotspot/coldspot. Sehingga, tidak perlu mengeluarkan banyak waktu, biaya dan tenaga untuk melakukan pencegahan dan intervensi di semua wilayah. Namun, tentu tetap perlu dilakukan pengawasan pada wilayah lainnya agar tidak menjadi outlier maupun hotspot/coldspot baru.

Keterbatasan dari analisis autokorelasi spasial

  • MUAP (Modifiable Areal Unit Problem) Hasil analisis dapat berubah apabila batas wilayah terjadi perubahan. Hal ini artinya analisis lebih dipengaruhi oleh unit analisis dibandingkan dengan fenomena asli yang terjadi.

  • Ukuran bobot spasial (Spatial Weights Matrix) Hasil sensitif pada asumsi pendekatan yang dipilih. Pada kasus ini digunakan pendekatan queen neighbours, tapi hasil analisis dapat berubah apabila asusmsi pendekatan diubah, misal rook atau k-nearest.

  • Masalah Multiple Testing pada analisis lokal Pada perhitungan autokorelasi LISA, dilakukan cukup banyak uji statistik di tiap wilayahnya. Hal ini cukup besar memungkinkan terjadinya false positive atau kesalahan pendeteksian hootspot/otlier. Dengan begitu disarakan untuk melakukan koreksi misal degan Bonferroni atau FDR.