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.
Sebelum rumus, definisikan beberapa notasi umum:
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 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 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 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).
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.
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)
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)
# 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.
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_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)
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.
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.
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.
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.