Bagian A. Konsep Teori

1. Jelaskan dengan kata-kata Anda apa yang dimaksud dengan autokorelasi spasial positif dan autokorelasi spasial negatif.

Autokorelasi spasial mengindikasikan bahwa nilai suatu variabel pada daerah tertentu terkait/dipengaruhi oleh nilai variabel tersebut pada daerah lain yang letaknya berdekatan (bertetangga). Autokorelasi spasial positif artinya wilayah yang berdekatan memiliki nilai data (misalnya jumlah kasus) yang serupa, tidak berbeda jauh/hampir sama. Sedangkan autokorelasi spasial negatif berarti lokasi yang berdekatan memiliki nilai data yang cenderung berbeda jauh atau berketerbalikan.

2. Berikan contoh sederhana dari fenomena yang bisa menunjukkan masingmasing pola tersebut.

Contoh autokorelasi spasial positif:

Contohnya adalah jumlah kasus COVID di kota yang padat penduduk. Kecamatan A memiliki jumlah kasus yang tinggi dan dikelilingi oleh kecamatan B, C, D dengan jumlah kasus tinggi pula. Hal ini mungkin terjadi karena kota tersebut memiliki tingkat kebersihan atau sanitasi yang mirip (biasanya rendah) sehingga kasus mudah menyebar luas. Kasus ini dapat menjadi salah satu contoh dari autokorelasi positif.
Contoh autokorelasi spasial negatif:

Misalnya, ada sebuah kawasan industri di Kalimantan dengan tingkat polusi tinggi namun industri ini dikelilingi oleh hutan-hutan lebat. Akibatnya, kualitas udara (misalnya kandungan CO2) yang tinggi dari industri langsung berbatasan dengan kualitas udara yang baik (CO2 rendah) dan membentuk pola yang berkebalikan. Kasus ini akan membentuk autokorelasi spasial negatif atau pola checkerboard.

3. Tuliskan rumus matematis dari:

Moran’s I

Rumus matematis untuk statistik Moran’s I adalah sebagai berikut:

\[ I = \frac{\frac{\sum_{i} \sum_{j} w_{ij} z_i \cdot z_j}{S_0}}{\frac{\sum_{i} z_i^2}{n}} \]

Keterangan:

\(I\) : Nilai Indeks Moran’s I.
\(n\) : Jumlah total observasi atau lokasi.
\(z_i = x_i - \bar{x}\) : Deviasi nilai observasi di lokasi \(i\) dari nilai rata-ratanya.
\(z_j = x_j - \bar{x}\) : Deviasi nilai observasi di lokasi \(j\) dari nilai rata-ratanya.
\(w_{ij}\) : Elemen dalam matriks pembobot spasial (spatial weights matrix).
\(S_0 = \sum_{i} \sum_{j} w_{ij}\) : Jumlah total dari semua bobot spasial.

Interpretasi:

\(I > 0\) : Autokorelasi spasial positif.
\(I < 0\) : Autokorelasi spasial negatif.
\(I \approx 0\) : Tidak ada pola (mendekati acak).

Geary’s C

Rumus matematis untuk statistik Geary’s C adalah sebagai berikut:

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

Keterangan:

\(C\) : Nilai Indeks Geary’s C.
\(N\) : Jumlah total unit spasial (lokasi).
\(x_i\) : Nilai variabel pada lokasi \(i\).
\(\bar{x}\) : Nilai rata-rata keseluruhan variabel \(x\).
\(w_{ij}\) : Elemen dalam matriks bobot spasial (spatial weights matrix), dengan \(w_{ii}=0\).
\(S_0 = \sum_{i} \sum_{j} w_{ij}\) : Jumlah total dari semua bobot spasial.

Interpretasi:

\(C \approx 1\) : Tidak ada autokorelasi spasial (nilai acak).
\(C < 1\) : Autokorelasi spasial positif.
\(C > 1\) : Autokorelasi spasial negatif.

Local Moran’s \(I_i\)

Rumus matematis untuk statistik Local Moran’s I adalah sebagai berikut:

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

dengan

\[ m_2 = \frac{1}{N}\sum_{k=1}^{N}(x_k - \bar{x})^2 \]

Keterangan:

\(I_i\) : Nilai Local Moran’s I pada lokasi \(i\).
\(N\) : Jumlah total observasi atau lokasi.
\(x_i\) : Nilai observasi di lokasi \(i\).
\(\bar{x}\) : Nilai rata-rata dari seluruh observasi.
\(w_{ij}\) : Elemen dalam matriks bobot spasial (bisa berupa queen, rook, atau k-nearest neighbors).
\(m_2 = \tfrac{1}{N}\sum_{k=1}^{N}(x_k - \bar{x})^2\) : Rata-rata kuadrat deviasi dari variabel \(x\).

Selain itu, Local Moran’s I memenuhi:

\[ I = \sum_{i=1}^{n} I_i \]

yang menunjukkan bahwa jumlah dari seluruh Local Moran’s I sama dengan Global Moran’s I.

Interpretasi:

\(I_i\) digunakan untuk mengukur tingkat autokorelasi spasial lokal di sekitar lokasi \(i\).
Nilai \(I_i\) yang signifikan positif menunjukkan adanya autokorelasi spasial positif.
Nilai \(I_i\) yang signifikan negatif menunjukkan adanya autokorelasi spasial negatif.

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

Definisi asli dari statistik lokal Getis–Ord tidak terstandar adalah rasio massa tetangga dengan massa total (Ord & Getis, 1995). Untuk suatu band jarak \(d\) (atau bobot ketetanggaan biner) dengan \(w_{ij}(d) \in \{0,1\}\):

Tanpa memasukkan lokasi \(i\) (local \(G_i\)):

\[ G_i(d) = \frac{\sum_{j \neq i} w_{ij}(d)\, x_j}{\sum_{j \neq i} x_j} \]

Memasukkan lokasi \(i\) (local \(G_i^*\)):

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

Keterangan:

\(G_i, G_i^*\) : Nilai indeks lokal Getis–Ord untuk lokasi \(i\).
\(x_j\) : Nilai variabel yang diamati pada lokasi \(j\).
\(w_{ij}(d)\) : Elemen dalam matriks bobot spasial (spatial weights matrix) yang menunjukkan kedekatan lokasi \(i\) dengan \(j\).
\(d\) : Jarak atau band jarak yang digunakan untuk menentukan tetangga.

Interpretasi:

Nilai \(G_i\) atau \(G_i^*\) yang lebih besar dari rata-rata menunjukkan adanya hot spot (high-high cluster).
Nilai yang lebih kecil dari rata-rata menunjukkan adanya cold spot (low-low cluster).
Signifikansi statistik biasanya diuji menggunakan uji analitik atau permutation test.

4. Jelaskan perbedaan utama antara ukuran global dan lokal.

Perbedaan antara ukuran global dan lokal disajikan dalam tabel berikut:

Aspek	Ukuran Global	Ukuran Lokal
Cakupan Analisis	Menggambarkan pola keseluruhan wilayah	Mengidentifikasi pola pada setiap lokasi
Hasil	Satu nilai ringkas (mis. Moran’s I, Geary’s C)	Nilai berbeda untuk tiap unit wilayah (mis. Local Moran’s I_i, Getis–Ord G_i^*)
Pertanyaan	“Apakah ada autokorelasi spasial secara umum?”	“Di mana lokasi klaster, outlier, hot spot, atau cold spot?”
Sensitivitas	Kurang sensitif terhadap variasi lokal	Lebih peka menangkap heterogenitas spasial
Interpretasi	Menyatakan ada/tidak dan arah autokorelasi (positif/negatif/acak)	Menunjukkan pola detail: High–High, Low–Low, High–Low, Low–High, hot spot, cold spot
Signifikansi	Satu uji statistik untuk seluruh wilayah	Uji statistik dilakukan per lokasi (perlu hati-hati dengan multiple testing)

Bagian B. Analisis Data (Simulasi)

1. Gunakan data spasial kecamatan di Kota Bandung (atau gunakan grid simulasi jika data asli tidak tersedia).

Data yang digunakan dalam analisis ini adalah data simulasi Jumlah Balita Pneumonia.

# Persiapan library
library(sf)
library(dplyr)
library(tmap)
library(sp)
library(spdep)
library(ggplot2)
library(mapview)
library(spData)
setwd("D:/OneDrive - Universitas Padjadjaran/SMT 5/SPASIAL/Data Spatial")

# Peta administratif Indonesia level kecamatan (gadm36 level 3)
Indo_Kec <- readRDS('gadm36_IDN_3_sp.rds')
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung",]
Bandung$id <- c(1:30)
BDG_sf <- st_as_sf(Bandung)
plot(Bandung, main = "Batasan Kecamatan di Kota Bandung")

2. Buat data simulasi kasus penyakit menular (misalnya diare per 10.000 penduduk) untuk 30 kecamatan di Kota Bandung. Gunakan distribusi Poisson dengan rata-rata berbeda antar kecamatan.

Berikut pembuatan data simulasi jumlah kasus balita pneumonia yang diatur agar ada autokorelasi signifikan untuk tujuan pembelajaran.

set.seed(77)
BDG_sf <- st_make_valid(BDG_sf)
BDG_sp <- as(BDG_sf, "Spatial")
nb_queen <- spdep::poly2nb(BDG_sp, queen = TRUE)
lw_queen <- spdep::nb2listw(nb_queen, style = "W", zero.policy = TRUE)
cent <- st_centroid(BDG_sf)
xy   <- st_coordinates(cent)
x    <- as.numeric(scale(xy[,1]))
y    <- as.numeric(scale(xy[,2]))
POP_U5 <- round(runif(nrow(BDG_sf), min = 2000, max = 8000))
z0 <- rnorm(nrow(BDG_sf), 0, 1)
z1 <- lag.listw(lw_queen, z0, zero.policy = TRUE)
z2 <- lag.listw(lw_queen, z1, zero.policy = TRUE)
z3 <- lag.listw(lw_queen, z2, zero.policy = TRUE)
grad   <- 0.8*x - 0.6*y
smooth <- 0.7*z1 + 0.6*z2 + 0.5*z3
s_base <- as.numeric(scale(grad + smooth))
make_cases <- function(alpha = 42, delta_LH = 2.4, delta_HL = 1.6, n_LH = 3, n_HL = 2) {
  rate_base <- 150
  rate_raw  <- rate_base + alpha * s_base
  lag_s  <- lag.listw(lw_queen, s_base, zero.policy = TRUE)
  idx_LH <- order(lag_s, decreasing = TRUE)[1:n_LH]                
  rate_raw[idx_LH] <- rate_raw[idx_LH] - delta_LH * sd(rate_raw)
  idx_HL <- setdiff(order(lag_s, decreasing = FALSE)[1:n_HL], idx_LH)
  rate_raw[idx_HL] <- rate_raw[idx_HL] + delta_HL * sd(rate_raw)
  RATE_U5_10K <- pmax(30, rate_raw + rnorm(length(rate_raw), 0, 4))
  lambda      <- RATE_U5_10K/10000 * POP_U5
  CASES_U5    <- rpois(length(lambda), lambda)
  list(RATE = RATE_U5_10K, CASES = CASES_U5, idx_LH = idx_LH, idx_HL = idx_HL)
}
classify_lisa <- function(y) {
  z     <- as.numeric(scale(y))
  lag_z <- lag.listw(lw_queen, z, zero.policy = TRUE)
  lm_df <- as.data.frame(localmoran(z, lw_queen, zero.policy = TRUE))
  if ("Pr(z > 0)" %in% names(lm_df)) {
    p_right <- lm_df[,"Pr(z > 0)"]; p_left <- 1 - p_right
    p_two   <- 2*pmin(p_right, p_left)
  } else if ("Pr(z <= 0)" %in% names(lm_df)) {
    p_left  <- lm_df[,"Pr(z <= 0)"]; p_right <- 1 - p_left
    p_two   <- 2*pmin(p_right, p_left)
  } else if ("Pr(z != E(Ii))" %in% names(lm_df)) {
    p_two   <- lm_df[,"Pr(z != E(Ii))"]
  } else {
    zval <- lm_df[,"Z.Ii"]; p_two <- 2*pnorm(-abs(zval))
  }
  quad <- ifelse(z >= 0 & lag_z >= 0, "High-High",
                 ifelse(z < 0 & lag_z < 0, "Low-Low",
                        ifelse(z >= 0 & lag_z < 0, "High-Low", "Low-High")))
  cluster <- ifelse(p_two <= 0.05, quad, "Not Significant")
  list(cat = cluster, p = p_two, Ii = lm_df$Ii, z = z, lag_z = lag_z)
}
best <- NULL
for (a in c(42, 48, 54, 62)) {
  sim <- make_cases(alpha = a, delta_LH = 2.4, delta_HL = 1.6, n_LH = 3, n_HL = 2)
  cls <- classify_lisa(sim$CASES)
  tab <- table(factor(cls$cat, levels=c("High-High","Low-Low","High-Low","Low-High","Not Significant")))
  n_sig <- sum(tab[c("High-High","Low-Low","High-Low","Low-High")], na.rm=TRUE)
  n_LH  <- if ("Low-High" %in% names(tab)) as.integer(tab["Low-High"]) else 0
  best <- list(alpha=a, sim=sim, cls=cls, tab=tab, n_sig=n_sig, n_LH=n_LH)
  if (n_sig >= 12 && n_LH >= 2) break
}
BDG_sf <- BDG_sf |>
  dplyr::mutate(
    POP_U5      = POP_U5,
    RATE_U5_10K = best$sim$RATE,
    CASES_U5    = best$sim$CASES)
z_cases <- as.numeric(scale(BDG_sf$CASES_U5))
lag_z   <- lag.listw(lw_queen, z_cases, zero.policy = TRUE)
lisa_df <- as.data.frame(localmoran(z_cases, lw_queen, zero.policy = TRUE))
if ("Pr(z > 0)" %in% names(lisa_df)) {
  p_right <- lisa_df[,"Pr(z > 0)"]; p_left <- 1 - p_right
  p_two   <- 2*pmin(p_right, p_left)
} else if ("Pr(z <= 0)" %in% names(lisa_df)) {
  p_left  <- lisa_df[,"Pr(z <= 0)"]; p_right <- 1 - p_left
  p_two   <- 2*pmin(p_right, p_left)
} else if ("Pr(z != E(Ii))" %in% names(lisa_df)) {
  p_two   <- lisa_df[,"Pr(z != E(Ii))"]
} else {
  zval <- lisa_df[,"Z.Ii"]; p_two <- 2*pnorm(-abs(zval))
}
quad_raw <- ifelse(z_cases >= 0 & lag_z >= 0, "High-High",
                   ifelse(z_cases < 0 & lag_z < 0, "Low-Low",
                          ifelse(z_cases >= 0 & lag_z < 0, "High-Low", "Low-High")))
LISA_cat <- ifelse(p_two <= 0.05, quad_raw, "Not Significant")
BDG_sf <- BDG_sf |>
  mutate(
    z_cases  = z_cases,
    lag_z    = lag_z,
    Ii       = lisa_df$Ii,
    p_lisa   = p_two,
    LISA_cat = factor(LISA_cat,
                      levels = c("High-High","Low-Low","High-Low","Low-High","Not Significant")))
head(BDG_sf %>%
       st_drop_geometry() %>%
       dplyr::select(NAME_3, POP_U5, CASES_U5, RATE_U5_10K))

##               NAME_3 POP_U5 CASES_U5 RATE_U5_10K
## 6199           Andir   3748       20    93.36067
## 6210        Antapani   6305      118   181.16525
## 6221       Arcamanik   7174      161   200.69077
## 6223     Astanaanyar   7700      108   125.94976
## 6224 Babakan Ciparay   6433       89   145.39484
## 6225   Bandung Kidul   4745       76   173.67572

3. Buat peta choropleth dari data simulasi tersebut.

Pemetaan Kasus (Continuous Mapping)

#Pemetaan kontinu
ggplot(BDG_sf) +
  geom_sf(aes(fill = RATE_U5_10K), color = "white", size = 0.2) +
  scale_fill_gradient(
    low  = "#f5f5dc", 
    high = "#08306b") +
  labs(
    title   = "Rate Balita Pneumonia per 10.000 Balita",
    fill    = "Kasus / 10.000",
    caption = "Catatan: Data simulasi berdistribusi Poisson dengan pola spasial halus + hotspot") +
  theme_minimal(base_size = 12)

Pemetaan Diskrit (Discrete Mapping)

# Interval kategori
breaks <- c(-Inf, 100, 150, 200, Inf)
labels <- c("Very Low", "Low", "High", "Very High")
# Buat variabel diskrit
BDG_sf$Pneumonia_Discrete <- cut(
  BDG_sf$RATE_U5_10K,
  breaks = breaks,
  labels = labels,
  right = TRUE)
# Pemetaan diskrit
ggplot() +
  geom_sf(data = BDG_sf, aes(fill = Pneumonia_Discrete), color = "white") +
  theme_bw() +
  scale_fill_manual(
    values = c(
      "Very Low"  = "#f5f5dc",
      "Low"       = "#d0dbe8",
      "High"      = "#5b9bd5",
      "Very High" = "#08306b" )) +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    legend.position = "right") +
  labs(
    title = "Kategori Rate Pneumonia Balita per 10.000",
    fill  = "Kategori")

4. Apa pola spasial yang terlihat secara visual?

Berdasarkan peta choropleth yang dihasilkan, dapat diketahui beberapa hal diantaranya:

Bagian tenggara terindikasi memiliki tingkat jumlah kasus balita pneumonia yang tinggi (Very High) dibandingkan wilayah lainnya. Hal ini mungkin mengindikasikan daerah hot-spot.
Daerah barat dan barat laut menunjukkan jumlah balita pneumonia yang sedikit (Low dan Very Low).
Meskipun hampir dikelilingi wilayah dengan kategori tinggi, daerah timur terindikasi memiliki jumlah balita pneumonia yang rendah (Very Low) sehingga mungkin menunjukkan autokorelasi negatif.
Distribusi spasial menunjukkan adanya pola pengelompokan (clustering), bukan sebaran acak.
Warna peta memperlihatkan adanya gradien spasial, di mana kasus meningkat dari barat laut menuju tenggara.

Bagian C. Pengukuran Autokorelasi

1. Hitung Moran’s I untuk data kasus balita pneumonia di Kota Bandung.

Berapa nilai Moran’s I?

# Queen contiguity
y <- BDG_sf$RATE_U5_10K
moran_queen <- moran.test(y, lw_queen)
moran_queen

## 
##  Moran I test under randomisation
## 
## data:  y  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 4.0131, p-value = 2.997e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.42306745       -0.03448276        0.01299941

# Rook contiguity
W <- poly2nb(BDG_sf, row.names = BDG_sf$NAME_3, queen = FALSE)
WL <- nb2listw(W)
moran_rook <- moran.test(y, WL)
moran_rook

## 
##  Moran I test under randomisation
## 
## data:  y  
## weights: WL    
## 
## Moran I statistic standard deviate = 3.8071, p-value = 7.031e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.41266955       -0.03448276        0.01379514

Hasil uji menggunakan bobot matriks Queen menghasilkan nilai Moran’s I sebesar 0.423 sedangkan jika menggunakan bobot matriks Rook menghasilkan nilai Moran’s I sebesar 0.413. Perbedaan ini disebabkan karena Rook contiguity menganggap dua wilayah bertetangga jika berbagi sisi saja sedangkan Queen contiguity menganggap dua wilayah bertetangga jika berbagi sisi atau sudut.

Nilai Moran’s I sebesar 0.423 (Queen) yang lebih besar dari 1 menunjukkan bahwa ada autokorelasi spasial positif dan cenderung tidak acak (ada pola clustering).

Untuk kedepannya akan digunakan Queen contiguity karena lebih dapat menangkap interaksi antarwilayah.

Apakah signifikan secara statistik (uji permutasi)?

Hipotesis

H₀: Tidak terdapat autokorelasi spasial pada data pneumonia balita (pola sebaran acak).

H₁: Terdapat autokorelasi spasial positif pada data pneumonia balita (pola clustering).

Taraf Signifikansi

α = 0.05

Statistik Uji

# Uji permutasi (Monte Carlo)
moran_perm <- moran.mc(y, lw_queen, nsim = 999)
moran_perm

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  y 
## weights: lw_queen  
## number of simulations + 1: 1000 
## 
## statistic = 0.42307, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Kriteria Uji

Jika p-value ≤ α, maka H₀ ditolak.

Keputusan

Nilai p-value (0.001) < α (0.05) maka H₀ ditolak.

Kesimpulan

Berdasarkan taraf signifikansi 5%, terdapat autokorelasi spasial positif yang signifikan pada data jumlah balita pneumonia di Kota Bandung, sehingga menunjukkan adanya pola pengelompokan (clustering).

Apa artinya bagi pola spasial penyakit?

Hasil uji autokorelasi spasial menunjukkan bahwa pola sebaran kasus pneumonia balita di Kota Bandung tidak bersifat acak, melainkan membentuk pola pengelompokan (clustering). Artinya, wilayah dengan angka kejadian tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki angka tinggi, demikian pula wilayah dengan angka rendah cenderung berdekatan dengan wilayah berangka rendah. Kondisi ini mengindikasikan adanya faktor lingkungan, sosial, atau demografis yang serupa antarwilayah bertetangga yang memengaruhi risiko penyakit, sehingga penyakit tidak tersebar merata tetapi terkonsentrasi pada area tertentu.

2. Hitung Geary’s C.

Nilai Geary’s C dihitung dengan cara berikut.

y <- BDG_sf$RATE_U5_10K
# Uji Asimtotik
geary_asym <- geary.test(y, lw_queen, randomisation = TRUE)
geary_asym

## 
##  Geary C test under randomisation
## 
## data:  y 
## weights: lw_queen   
## 
## Geary C statistic standard deviate = 3.4952, p-value = 0.0002369
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.58803993        1.00000000        0.01389233

# Uji Permutasi (Monte Carlo)
set.seed(77)
geary_perm <- geary.mc(y, lw_queen, nsim = 999)
geary_perm

## 
##  Monte-Carlo simulation of Geary C
## 
## data:  y 
## weights: lw_queen  
## number of simulations + 1: 1000 
## 
## statistic = 0.58804, observed rank = 2, p-value = 0.002
## alternative hypothesis: greater

Berdasarkan uji diperoleh nilai Geary C sebesar 0.588 dengan p-value sebesar 0.002 yang kurang dari α (0.05) sehingga diperoleh keputusan tolak H₀. Artinya terdapat autokorelasi spasial positif yang signifikan.

Bagaimana perbandingannya dengan Moran’s I?

Nilai Moran’s I dan Geary C menunjukkan hasil yang konsisten, yaitu adanya autokorelasi spasial positif yang signifikan yang artinya terdapat pola clustering pada data jumlah balita pneumonia. Bedanya, Moran’s I berfokus pada pola secara global sedangkan Geary C lebih berfokus pada pola antar wilayah yang bertetangga namun keduanya masih sama-sama merupakan ukuran global.

Jelaskan perbedaan sensitivitas kedua ukuran ini.

Keduanya sama sama merupakan ukuran global namun berbeda di tingkat sensitivitasnya.

Moran’s I dihitung dengan membandingkan deviasi nilai suatu wilayah terhadap rata-rata keseluruhan, sehingga secara statistik mirip dengan koefisien korelasi Pearson. Rumusnya berbasis kovarians antara suatu lokasi dengan tetangganya terhadap rata-rata global. Oleh karena itu, Moran’s I lebih menekankan pada pola spasial secara menyeluruh, sehingga sensitif untuk mengukur autokorelasi global tetapi kurang sensitif terhadap variasi kecil atau heterogenitas lokal.
Geary’s C menghitung selisih kuadrat antar nilai wilayah yang bertetangga langsung. Rumus ini membuat Geary’s C lebih peka terhadap perbedaan lokal karena menilai kesamaan atau ketidaksamaan nilai antar pasangan tetangga. Secara statistik menyerupai ukuran mean squared difference, sehingga lebih efektif mendeteksi heterogenitas lokal atau variasi kecil antar wilayah dibanding Moran’s I.

3. Hitung Local Moran’s I (LISA).

Nilai Local Moran’s (LISA) dihitung dengan cara berikut.

y      <- BDG_sf$RATE_U5_10K
z      <- as.numeric(scale(y))
lag_z  <- lag.listw(lw_queen, z, zero.policy = TRUE)
# Local Moran's I
lisa_df <- as.data.frame(localmoran(z, lw_queen, zero.policy = TRUE))
if ("Pr(z > 0)" %in% names(lisa_df)) {
  p_right <- lisa_df[,"Pr(z > 0)"]
  p_left  <- 1 - p_right
  p_two   <- 2*pmin(p_right, p_left)} else if ("Pr(z <= 0)" %in% names(lisa_df)) {
  p_left  <- lisa_df[,"Pr(z <= 0)"]
  p_right <- 1 - p_left
  p_two   <- 2*pmin(p_right, p_left)} else if ("Pr(z != E(Ii))" %in% names(lisa_df)) {
  p_two   <- lisa_df[,"Pr(z != E(Ii))"]} else {
  zval <- lisa_df[,"Z.Ii"]
  p_two <- 2*pnorm(-abs(zval))}
# Klasifikasi kuadran LISA
quad_raw <- ifelse(z >= 0 & lag_z >= 0, "High-High",
                   ifelse(z < 0 & lag_z < 0, "Low-Low",
                          ifelse(z >= 0 & lag_z < 0, "High-Low", "Low-High")))
sig_thr  <- 0.05
cluster  <- ifelse(p_two <= sig_thr, quad_raw, "Not Significant")
BDG_sf <- BDG_sf |>
  dplyr::mutate(
    z_cases   = z,
    lag_z     = lag_z,
    Ii        = lisa_df$Ii,
    Z_Ii      = lisa_df$Z.Ii,
    p_lisa    = p_two,
    LISA_cat  = factor(cluster,
                       levels = c("High-High","Low-Low","High-Low","Low-High","Not Significant")))
lisa_tabel <- BDG_sf %>%
  st_drop_geometry() %>%
  group_by(LISA_cat) %>%
  summarise(
    Kecamatan = paste(NAME_3, collapse = ","),
    Jumlah    = n())

Identifikasi kecamatan yang masuk kategori High-High, Low-Low, High-Low, dan Low-High.

Kecamatan berdasarkan kategorinya dapat dilihat melalui tabel berikut.

lisa_tabel

## # A tibble: 4 × 3
##   LISA_cat        Kecamatan                                               Jumlah
##   <fct>           <chr>                                                    <int>
## 1 High-High       Antapani,Arcamanik,Buahbatu,Mandalajati,Rancasari            5
## 2 Low-Low         Bandung Wetan,Cicendo,Cidadap,Coblong,Sukajadi               5
## 3 Low-High        Cinambo                                                      1
## 4 Not Significant Andir,Astanaanyar,Babakan Ciparay,Bandung Kidul,Bandun…     19

Interpretasi

High-High (HH) adalah wilayah dengan nilai tinggi yang dikelilingi oleh wilayah bernilai tinggi.
Low-Low (LL) adalah wilayah dengan nilai rendah yang dikelilingi oleh wilayah bernilai rendah.
High-Low (HL) adalah wilayah dengan nilai tinggi namun dikelilingi oleh wilayah bernilai rendah, disebut outlier positif.
Low-High (LH) adalah wilayah dengan nilai rendah tetapi dikelilingi oleh wilayah bernilai tinggi, disebut outlier negatif.

Buat peta cluster LISA.

peta_lisa <- ggplot(BDG_sf) +
  geom_sf(aes(fill = LISA_cat), color = "white", size = 0.2) +
  scale_fill_manual(
    values = c(
      "High-High"      = "#D7191C",
      "Low-Low"        = "#1B9E77",
      "Low-High"       = "#2C7FB8",
      "Not significant"= "grey99")) +
  labs(title   = "LISA Plot – Rate Balita Pneumonia per 10.000",
    fill    = "Cluster") + theme_minimal(base_size = 12); peta_lisa

Apa interpretasi hasil ini untuk kasus penyakit menular?

Wilayah kategori High-High yaitu Antapani, Arcamanik, Buahbatu, Mandalajati, dan Rancasari perlu menjadi prioritas dalam penanganan dan upaya penanggulangan penyebaran penyakit pneumonia pada balita karena jumlah kasusnya lebih marak dan berpusat di area ini. Pemerintah disarankan untuk lebih fokus dalam mengimplementasikan upaya pencegahan dan penularan seperti peningkatan pelayanan kesehatan dan edukasi masyarakat di kelima kecamatan ini.

Wilayah kategori Low-Low yaitu Bandung Wetan, Cicendo, Cidadap, Coblong, dan Sukajadi memiliki jumlah kasus balita pneumonia yang sedikit sehingga relatif aman namun tetap harus dijaga baik dari akses layanan kesehatannya maupun perilaku hidup bersih di masyarakat.

Wilayah kategori Low-High yaitu Cinambo perlu diteliti lebih lanjut karena memiliki tingkat kasus yang rendah padahal dikelilingi oleh wilayah dengan kasus tinggi. Hal ini mungkin dikarenakan adanya perbedaan yang cukup besar dari segi fasilitas kesehatan maupun perilaku masyarakatnya dan tetap harus dijaga agar tidak tertular dari wilayah tetangganya.

Wilayah yang termasuk kategori Not Significant artinya tidak memiliki autokorelasi spasial yang signifikan atau dengan kata lain wilayah tersebut tidak dapat membentuk klaster yang jelas.

4. Hitung Getis–Ord 𝐺∗𝑖

Nilai Getis–Ord 𝐺∗𝑖 dihitung dengan cara berikut ini.

# Identifikasi Variabel
x_raw <- BDG_sf$RATE_U5_10K
sum_x <- sum(x_raw)
Wb <- spdep::listw2mat(lw_queen)
# 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) — set w_ii = 1
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
Gz <- spdep::localG(x_raw, listw = lw_queen, zero.policy = TRUE)
BDG_sf <- dplyr::mutate(BDG_sf,
  G_raw      = G_raw,
  G_star_raw = G_star_raw,
  z_Gistar   = as.numeric(Gz),
  GSTAR_cat  = dplyr::case_when(
    z_Gistar >=  1.96 ~ "Hot spot",
    z_Gistar <= -1.96 ~ "Cold spot",
    TRUE              ~ "Not Significant"))
# Ringkasan Nilai
summary(dplyr::select(BDG_sf, G_raw, G_star_raw, z_Gistar))

##      G_raw           G_star_raw         z_Gistar                 geometry 
##  Min.   :0.02013   Min.   :0.03645   Min.   :-2.75934   MULTIPOLYGON :30  
##  1st Qu.:0.02890   1st Qu.:0.05551   1st Qu.:-1.29368   epsg:NA      : 0  
##  Median :0.03395   Median :0.06344   Median :-0.12187   +proj=long...: 0  
##  Mean   :0.03443   Mean   :0.06655   Mean   : 0.01837                     
##  3rd Qu.:0.03879   3rd Qu.:0.08013   3rd Qu.: 0.78893                     
##  Max.   :0.05549   Max.   :0.10780   Max.   : 3.49744

Interpretasi

Hasil ringkasan menunjukkan bahwa nilai Getis–Ord 𝐺∗𝑖 bervariasi antar kecamatan, dengan Z-score berkisar dari -2.76 hingga 3.50. Adanya nilai Z positif tinggi menandakan hotspot signifikan (wilayah dengan kasus pneumonia tinggi berdekatan dengan wilayah tinggi), sedangkan nilai Z negatif rendah menandakan coldspot signifikan (wilayah kasus rendah berdekatan dengan wilayah rendah). Secara umum, distribusi kasus tidak acak, melainkan membentuk pola kluster panas dan dingin.

Tentukan kecamatan yang termasuk hot spot dan cold spot.

Kecamatan yang termasuk hot spot dan cold spot beserta nilai Getis–Ord 𝐺∗𝑖-nya dapat dilihat melalui tabel berikut.

# Tabel Identifikasi Kecamatan
BDG_sf %>% st_drop_geometry() %>% filter(GSTAR_cat %in% c("Hot spot", "Cold spot")) %>%
  transmute(
    Kecamatan = NAME_3,
    Kategori  = GSTAR_cat,
    Gi_star   = round(z_Gistar, 3)) %>%
  arrange(Kategori, dplyr::desc(if_else(Kategori == "Hot spot", Gi_star, -Gi_star)))

##        Kecamatan  Kategori Gi_star
## 1        Cicendo Cold spot  -2.759
## 2       Sukajadi Cold spot  -2.262
## 3        Coblong Cold spot  -2.243
## 4  Bandung Wetan Cold spot  -2.160
## 5        Cidadap Cold spot  -2.154
## 6      Arcamanik  Hot spot   3.497
## 7      Rancasari  Hot spot   3.274
## 8       Buahbatu  Hot spot   2.630
## 9       Antapani  Hot spot   2.313
## 10       Cinambo  Hot spot   2.135
## 11   Mandalajati  Hot spot   2.029

Dimana apabila suatu wilayah memiliki nilai Z-Score > 1.96 maka akan diidentifikasi sebagai hot spot sedangkan apabila nilai Z-Score < -1.96 maka diidentifikasi sebagai cold spot. Kemudian untuk memperjelas informasi, berikut choropleth untuk wilayah hot spot maupun cold spot.

cols_gstar <- c("Hot spot" = "#D7191C", "Cold spot" = "#1B9E77", "Not significant" = "grey99")
ggplot(BDG_sf) +
  geom_sf(aes(fill = GSTAR_cat), color = "white", size = 0.2) +
  scale_fill_manual(
    values = cols_gstar,
    breaks = c("Hot spot", "Cold spot", "Not significant"), drop   = FALSE) +
  labs(
    title   = "Getis–Ord G* Plot – Rate Balita Pneumonia per 10.000",
    fill    = "Kategori") +  theme_minimal(base_size = 12)

Bandingkan hasilnya dengan peta LISA.

Lihat kembali peta LISA yang sebelumnya dihasilkan berikut ini.

peta_lisa

Terlihat ada satu wilayah yang teridentifikasi berbeda, untuk lebih jelasnya hasil Getis–Ord 𝐺∗𝑖 dengan LISA disajikan pada tabel berikut.

table(GSTAR = BDG_sf$GSTAR_cat, LISA = BDG_sf$LISA_cat)

##                  LISA
## GSTAR             High-High Low-Low High-Low Low-High Not Significant
##   Cold spot               0       5        0        0               0
##   Hot spot                5       0        0        1               0
##   Not Significant         0       0        0        0              19

Interpretasi

Tabel diatas menunjukkan bahwa uji Getis–Ord 𝐺∗𝑖 dan LISA menghasilkan pengelompokkan wilayah yang konsisten dimana sebanyak tepat 5 wilayah yang dikategorikan High-High pada LISA juga dikategorikan sebagai Hot spot pada Getis–Ord 𝐺∗𝑖.

Begitupun 5 wilayah Low-Low pada LISA dikategorikan sebagai Cold spot pada Getis–Ord 𝐺∗𝑖.

Namun terdapat satu wilayah yang dikategorikan sebagai Low-High pada LISA justru dikategorikan sebagai Hot spot pada Getis–Ord 𝐺∗𝑖. Hal ini disebabkan oleh perbedaan sensitivitas pengelompokkan oleh LISA serta Getis–Ord 𝐺∗𝑖 dimana LISA cenderung mengidentifikasi pola lokal antar tetangga sedangkan Getis–Ord 𝐺∗𝑖 berfokus pada besar nilai tinggi/rendah di sekitar unit observasi.

Adapun kedua metode dapat dengan tepat mengidentifikasi 19 wilayah yang tidak memiliki pola spasial dan tidak signifikan secara statistik.

Apakah ada perbedaan wilayah yang ditandai sebagai klaster signifikan?

Kecamatan yang ditandai secara berbeda di kedua metode dapat diketahui dengan cara berikut.

BDG_sf %>% st_drop_geometry() %>% filter(LISA_cat == "Low-High", GSTAR_cat == "Hot spot") %>% select(Kecamatan = NAME_3, LISA_cat, GSTAR_cat)

##   Kecamatan LISA_cat GSTAR_cat
## 1   Cinambo Low-High  Hot spot

Interpretasi

Ya. Kecamatan Cinambo oleh LISA dikategorikan sebagai Low-High, sedangkan oleh Getis–Ord G* justru teridentifikasi sebagai Hot spot. Hal ini menunjukkan adanya perbedaan sensitivitas antara kedua metode dalam mendeteksi pola klaster signifikan.

Tugas 1 Autokorelasi Spasial

Humaira Zeanova

2025-09-26

Bagian A. Konsep Teori

Moran’s I

Geary’s C

Local Moran’s \(I_i\)

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

Bagian B. Analisis Data (Simulasi)

Pemetaan Kasus (Continuous Mapping)

Pemetaan Diskrit (Discrete Mapping)

Bagian C. Pengukuran Autokorelasi

Berapa nilai Moran’s I?

Apakah signifikan secara statistik (uji permutasi)?

Apa artinya bagi pola spasial penyakit?

Bagaimana perbandingannya dengan Moran’s I?

Jelaskan perbedaan sensitivitas kedua ukuran ini.

Identifikasi kecamatan yang masuk kategori High-High, Low-Low, High-Low, dan Low-High.

Buat peta cluster LISA.

Apa interpretasi hasil ini untuk kasus penyakit menular?

Tentukan kecamatan yang termasuk hot spot dan cold spot.

Bandingkan hasilnya dengan peta LISA.

Apakah ada perbedaan wilayah yang ditandai sebagai klaster signifikan?

Bagian D. Diskusi Kritis

R Pubs