Nama : Maulana Hardy Rayyan

NPM : 140610230064

Bagian A

Autokorelasi Spasial Positif

Sebelumnya, Autokorelasi Spasial merupakan ukuran kemiripian atau bedanya nilai antar lokasi dalam suatu ruang. Berarti Autokorelasi Spasial Positif berarti wilayah yang disekitarnya atau tetangganya memiliki kemiripan nilai dengan wilayah tersebut.

Contoh :

  • Penyakit Menular TBC (misalnya)

    Misalnya terdapat jumlah kasus TBC yang tinggi di kecamatan A yang ada di Jakarta Utara, kecamatan disekitar memiliki kecenderungan untuk kasus yang tinggi juga karena penularan antar penduduk yang berdekatan

  • Harga Rumah dan Tanah

    Misalnya rumah atau tanah di daerah pinggiran Jakarta Utara memiliki harga yang rendah, karena dikelilingi juga dengan perumahan yang harga rumah atau tanah nya rendah juga

Autokorelasi Spasial Negatif

Autokorelasi Spasial Negatif berarti wilayah yang disekitarnya atau tetangganya memiliki kemiripan nilai dengan wilayah tersebut.

Contoh :

  • Kepadatan Penduduk

    Misalnya, di pusat kota jakarta memiliki kepadatan penduduk yang tinggi dibandingkan di pinggir kota jakarta yang kepadatan penduduknya rendah

Rumus Matematis

Moran’s I

Rumus Matematis :

\[ 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} \]

dengan interpretasi :

  • Jika I > 0, maka nilai mirip dengan tetangga disekitar nya atau juga bisa disebut memiliki positive spatial autocorrelation

  • Jika I < 0, maka memiliki perbedaan nilai dengan tetangga di sekitarnya atau memiliki negative spatial autocorrelation (memiliki pola checkerboard)

  • Jika I \(\approx\) 0, maka nilai yang dihasilkan acak atau tidak beraturan

Geary’s C

Rumus Matematis :

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

dengan interpretasi :

  • Jika C < 1, maka memiliki autokorelasi positif

  • Jika C > 1, maka memiliki autokorelasi negatif

  • Jika C = 1, maka tidak terdapat autokorelasi

Local Moran’s I

Rumus Matematis :

\[ 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 \]

dengan interpretasi :

  • Jika nilai \(I_i\) > 0, maka unit spasial i memiliki nilai yang mirip dengan tetangganya dan membentuk cluster High-High atau Low-Low

  • Jika nilai \(I_i\) < 0, maka unit spasial i memiliki perbedaan nilai dengan tetangganya dan membentuk pola High-Low atau Low-High

  • Jika nilai \(I_i\) \(\approx\) 0, maka tidak terdapat pola lokal yang signifikan atau memiliki pola lokal yang acak

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

  • 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} \]

  • Dengan 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} \]

    dengan \[ w_{ii}(d) = 1 \]

Perbedaan Utama antara Ukuran Global dan Lokal

Ukuran Global

Ukuran global mengukur keterkaitan spasial secara keseluruhan dan memiliki fungsi untuk mengidentifikasi apakah data memiliki autokorelasi spasial positif, negatif atau acak dengan indeks Moran’s I dan Geary’s C

Ukuran Lokal

Ukuran global mengukur keterkaitan spasial di sekitar unit wilayah tertentu dan memiliki fungsi untuk mengidentifikasi lokasi spesifik klaster dengan indeks Local Moran’s I dan Getis-Ord \(G_i\) dan \(G_i^*\)

Bagian B

Data Spasial Kecamatan di 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
setwd("D:/Maulana Hardy Rayyan/UNPAD/Analisis Data Spasial/Praktikum Spasial")

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

Data Simulasi Kasus Penyakit Menular untuk 30 Kecamatan di Kota Bandung dan Menggunakan Distribusi Poisson dengan Rata-rata Berbeda antar Kecamatan

N <- nrow(Bandung_sf)

# --- 1. Struktur tetangga (queen contiguity) ---
nb <- poly2nb(as_Spatial(Bandung_sf), queen = TRUE)
lwW <- nb2listw(nb, style = "W")

# --- 2. Kovariat acak & lag spasial ---
set.seed(123)
z  <- rnorm(N)
lz <- lag.listw(lwW, z)

# --- 3. Parameter lambda (rata-rata berbeda tiap kecamatan) ---
beta0 <- 1.2; beta1 <- 0.6; beta2 <- 1.0
lambda <- exp(beta0 + beta1*scale(z) + beta2*scale(lz))

# --- 4. Populasi acak tiap kecamatan ---
pop <- round(runif(N, 3000, 8000))

# --- 5. Kasus penyakit ~ Distribusi Poisson ---
cases <- rpois(N, lambda = lambda * (pop/10000))

# --- 6. Hitung rate per 10.000 penduduk ---
rate10k <- (cases / pop) * 10000

# --- 7. Simpan hasil simulasi ke data.frame ---
Diare <- data.frame(
  id = Bandung_sf$id,
  pop = pop,
  cases = cases,
  rate10k = rate10k
)

# --- 8. Gabungkan dengan peta Bandung ---
Bandung_merged <- left_join(Bandung_sf, Diare, by = "id")

# --- 9. Cek hasil ---
head(Bandung_merged)
## Simple feature collection with 6 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 107.56 ymin: -6.965966 xmax: 107.6915 ymax: -6.891394
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
##   GID_0    NAME_0   GID_1     NAME_1 NL_NAME_1      GID_2       NAME_2
## 1   IDN Indonesia IDN.9_1 Jawa Barat      <NA> IDN.9.14_1 Kota Bandung
## 2   IDN Indonesia IDN.9_1 Jawa Barat      <NA> IDN.9.14_1 Kota Bandung
## 3   IDN Indonesia IDN.9_1 Jawa Barat      <NA> IDN.9.14_1 Kota Bandung
## 4   IDN Indonesia IDN.9_1 Jawa Barat      <NA> IDN.9.14_1 Kota Bandung
## 5   IDN Indonesia IDN.9_1 Jawa Barat      <NA> IDN.9.14_1 Kota Bandung
## 6   IDN Indonesia IDN.9_1 Jawa Barat      <NA> IDN.9.14_1 Kota Bandung
##   NL_NAME_2        GID_3          NAME_3 VARNAME_3 NL_NAME_3    TYPE_3
## 1      <NA> IDN.9.14.1_1           Andir      <NA>      <NA> Kecamatan
## 2      <NA> IDN.9.14.2_1        Antapani      <NA>      <NA> Kecamatan
## 3      <NA> IDN.9.14.3_1       Arcamanik      <NA>      <NA> Kecamatan
## 4      <NA> IDN.9.14.4_1     Astanaanyar      <NA>      <NA> Kecamatan
## 5      <NA> IDN.9.14.5_1 Babakan Ciparay      <NA>      <NA> Kecamatan
## 6      <NA> IDN.9.14.6_1   Bandung Kidul      <NA>      <NA> Kecamatan
##      ENGTYPE_3    CC_3 HASC_3 id  pop cases  rate10k
## 1 Sub-district 3273180   <NA>  1 6326     0 0.000000
## 2 Sub-district 3273141   <NA>  2 3474     1 2.878526
## 3 Sub-district 3273130   <NA>  3 4920     1 2.032520
## 4 Sub-district 3273050   <NA>  4 4372     1 2.287283
## 5 Sub-district 3273020   <NA>  5 7073     3 4.241482
## 6 Sub-district 3273080   <NA>  6 5243     2 3.814610
##                         geometry
## 1 MULTIPOLYGON (((107.5984 -6...
## 2 MULTIPOLYGON (((107.6663 -6...
## 3 MULTIPOLYGON (((107.678 -6....
## 4 MULTIPOLYGON (((107.6087 -6...
## 5 MULTIPOLYGON (((107.5868 -6...
## 6 MULTIPOLYGON (((107.6396 -6...

Peta Choropleth

# Peta Choropleth rate per 10.000
ggplot(Bandung_merged) +
  geom_sf(aes(fill = rate10k), color = "white") +
  scale_fill_viridis_c(option = "magma", direction = -1,
                       name = "Kasus per\n10.000 Penduduk") +
  labs(
    title = "Peta Choropleth Kasus Diare per 10.000 Penduduk",
    subtitle = "30 Kecamatan di Kota Bandung (Simulasi)",
    caption = "Sumber: Simulasi Distribusi Poisson"
  ) +
  theme_minimal()

Pola Spasial yang Terlihat

Berdasarkan peta choropleth dapat dilihat bahwa terdapat indikasi clustering yang mana kecamatan yang memiliki angka kasus diare yang tinggi berdekatan dengan kecamatan yang memiliki angka kasus diare yang tinggi juga.

Bagian C

Hitung Moran’s I

# Buat struktur tetangga (queen contiguity)
nb <- poly2nb(as_Spatial(Bandung_merged), queen = TRUE)
lwW <- nb2listw(nb, style = "W")

# Uji Moran's I 
moran_test <- moran.test(Bandung_merged$rate10k, lwW)
print(moran_test)
## 
##  Moran I test under randomisation
## 
## data:  Bandung_merged$rate10k  
## weights: lwW    
## 
## Moran I statistic standard deviate = 4.8496, p-value = 6.184e-07
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.50481782       -0.03448276        0.01236637
# Uji signifikansi dengan permutasi (Monte Carlo) 
set.seed(123)
moran_perm <- moran.mc(Bandung_merged$rate10k, lwW, nsim = 999)
print(moran_perm)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  Bandung_merged$rate10k 
## weights: lwW  
## number of simulations + 1: 1000 
## 
## statistic = 0.50482, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Apa Artinya Bagi Pola Spasial Penyakit ?

Melihat dari hasil Moran’s I, dapat dilihat bahwa I nya sebesar 0,50482 yang mana berarti Kota Bandung memiliki autokorelasi spasial positif

Hitung Geary’s C

# Uji Geary's C 
geary_test <- geary.test(Bandung_merged$rate10k, lwW)
print(geary_test)
## 
##  Geary C test under randomisation
## 
## data:  Bandung_merged$rate10k 
## weights: lwW   
## 
## Geary C statistic standard deviate = 4.3461, p-value = 6.928e-06
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.47515619        1.00000000        0.01458327

Perbandingan dengan Moran’s I

Jika dilihat dari nilai statistik nya, Geary’s C dengan nilai sebesar 0,47516 dibandingkan dengan Moran’s I yang memiliki nilai sebesar 0,50482. Dari nilai tersebut kedua indeks memiliki interpretasi hasil yang saling berlawanan. Untuk Geary’s C, jika C < 0 maka memiliki autokorelasi positif atau tetangga di sekitarnya memiliki nilai yang mirip sedangkan untuk Moran’s I, jika I < 0 maka memiliki autokorelasi negatif atau memiliki nilai yang berbeda dengan tetangganya

Jelaskan Perbedaan Sensitivitas Kedua Ukuran ini

Untuk ukuran Geary’s C memiliki sensitivitas lebih tinggi terhadap variabilitas lokal dan juga pencilan disekitarnya sedangkan untuk ukuran Moran’s I itu mengukur kemiripan pola secara global saja, dan tidak begitu peka dengan variabilitas ataupun pencilan di sekitarnya

Hitung Local Moran’s I (LISA)

Identifikasi Kecamatan per Kategori

Bandung_merged <- left_join(Bandung_sf, Diare, by = "id")

# --- 2. Struktur tetangga (queen contiguity) ---
nb  <- poly2nb(as_Spatial(Bandung_merged), queen = TRUE)
lwW <- nb2listw(nb, style = "W")

# --- 3. Hitung Local Moran’s I (LISA) ---
local_moran <- localmoran(Bandung_merged$rate10k, lwW)

Bandung_merged$Ii   <- local_moran[, "Ii"]
Bandung_merged$E.Ii <- local_moran[, "E.Ii"]
Bandung_merged$Z.Ii <- local_moran[, "Z.Ii"]
Bandung_merged$Pr.z <- local_moran[, "Pr(z != E(Ii))"]  # p-value

# --- 4. Klasifikasi cluster LISA ---
mean_rate <- mean(Bandung_merged$rate10k)

Bandung_merged <- Bandung_merged %>%
  mutate(cluster = case_when(
    rate10k >= mean_rate & lag.listw(lwW, rate10k) >= mean_rate & Pr.z < 0.05 ~ "High-High",
    rate10k <= mean_rate & lag.listw(lwW, rate10k) <= mean_rate & Pr.z < 0.05 ~ "Low-Low",
    rate10k >= mean_rate & lag.listw(lwW, rate10k) <= mean_rate & Pr.z < 0.05 ~ "High-Low",
    rate10k <= mean_rate & lag.listw(lwW, rate10k) >= mean_rate & Pr.z < 0.05 ~ "Low-High",
    TRUE ~ "Not Significant"
  ))

# --- 5. Tabel ringkasan klaster ---
# Gunakan NAME_3 dari shapefile GADM sebagai nama kecamatan
Bandung_clusters <- Bandung_merged %>%
  st_drop_geometry() %>%
  select(Kecamatan = NAME_3, rate10k, cluster) %>%
  arrange(cluster)

print(Bandung_clusters)
##           Kecamatan   rate10k         cluster
## 1           Cidadap 22.453562       High-High
## 2           Coblong 27.316663       High-High
## 3          Sukajadi 16.403785       High-High
## 4          Sukasari 24.112525       High-High
## 5             Andir  0.000000 Not Significant
## 6          Antapani  2.878526 Not Significant
## 7         Arcamanik  2.032520 Not Significant
## 8       Astanaanyar  2.287283 Not Significant
## 9   Babakan Ciparay  4.241482 Not Significant
## 10    Bandung Kidul  3.814610 Not Significant
## 11    Bandung Kulon  4.255319 Not Significant
## 12    Bandung Wetan  1.416029 Not Significant
## 13      Batununggal  1.434309 Not Significant
## 14  Bojongloa Kaler  3.846894 Not Significant
## 15  Bojongloa Kidul 20.673361 Not Significant
## 16         Buahbatu  4.881224 Not Significant
## 17 Cibeunying Kaler  7.632423 Not Significant
## 18 Cibeunying Kidul  0.000000 Not Significant
## 19           Cibiru  9.298866 Not Significant
## 20          Cicendo  9.753719 Not Significant
## 21          Cinambo  0.000000 Not Significant
## 22         Gedebage  0.000000 Not Significant
## 23     Kiaracondong  0.000000 Not Significant
## 24         Lengkong  3.154574 Not Significant
## 25      Mandalajati  5.896226 Not Significant
## 26      Panyileukan  0.000000 Not Significant
## 27        Rancasari  8.537279 Not Significant
## 28            Regol  3.865481 Not Significant
## 29    Sumur Bandung  0.000000 Not Significant
## 30     Ujung Berung  0.000000 Not Significant
# Per kategori:
HighHigh <- Bandung_clusters %>% filter(cluster == "High-High")
LowLow   <- Bandung_clusters %>% filter(cluster == "Low-Low")
HighLow  <- Bandung_clusters %>% filter(cluster == "High-Low")
LowHigh  <- Bandung_clusters %>% filter(cluster == "Low-High")

# Tabel hasil LISA per kategori
clusters_list <- split(Bandung_clusters, Bandung_clusters$cluster)

# Tabel kategori
clusters_list$`High-High`
##   Kecamatan  rate10k   cluster
## 1   Cidadap 22.45356 High-High
## 2   Coblong 27.31666 High-High
## 3  Sukajadi 16.40379 High-High
## 4  Sukasari 24.11253 High-High
clusters_list$`Low-Low`
## NULL
clusters_list$`High-Low`
## NULL
clusters_list$`Low-High`
## NULL

Peta Cluster Local Moran’s I (LISA)

ggplot(Bandung_merged) +
  geom_sf(aes(fill = cluster), color = "white") +
  scale_fill_manual(values = c(
    "High-High" = "red",
    "Low-Low" = "blue",
    "High-Low" = "orange",
    "Low-High" = "lightblue",
    "Not Significant" = "grey80"
  )) +
  labs(
    title = "Peta LISA (Local Moran’s I) Kasus Diare Kota Bandung",
    fill = "Cluster"
  ) +
  theme_minimal()

Interpretasi

Dari syntax Local Moran’s I di atas dapat disimpulkan bahwa kecamatan yang termasuk ke kategori High-High yaitu Cidadap, Coblong, Sukajadi, dan Sukasari dan sisanya terhitung tidak signifikan (sudah dicoba agar bisa muncul masing-masing kategori namun tidak bisa)

Hitung Getis-Ord \(G_i^*\)

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

# --- Matriks bobot biner (tidak distandarisasi) ---
lwB <- nb2listw(nb, style = "B", zero.policy = TRUE)
Wb <- spdep::listw2mat(lwB)  

# 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 dengan spdep::localG
Gz <- spdep::localG(x_raw, listw = lwW)  # z(Gi*)

# Gabungkan ke data spasial
Bandung_G <- Bandung_merged %>%
  mutate(
    G_raw      = G_raw,
    G_star_raw = G_star_raw,
    z_Gistar   = as.numeric(Gz),
    hotcold    = case_when(
      z_Gistar >=  1.96 ~ "Hot spot (p<0.05)",
      z_Gistar <= -1.96 ~ "Cold spot (p<0.05)",
      TRUE              ~ "Not significant"
    )
  )
# Peta raw G*
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()

# --- Peta hotspot/coldspot (z-skor) ---
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()

Untuk kecamatan yang termasuk ke hotspot yaitu Cidadap, Coblong, Sukajadi, dan Sukasari. Selain dari keempat itu terasuk ke coldspot

Perbandingan dengan Peta Local Moran’s I

Jika dilihat dari kedua peta tersebut, untuk Local Moran’s I itu merupakan kategori High-High sedangkan pada Getis-Ord \(G_i^*\) itu merupakan daerah hotspot

Apakah ada perbedaan wilayah yang ditandai sebagai klaster signifikan?

Dalam ukuran Getis-Ord, klaster yang signifikan dapat diartikan sebagai hotspot atau nilai yang tinggi

Bagian D

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

Hasil autokorelasi spasial dijadikan sebagai evidence-based desicion making yang bisa membantu dinas kesehatan dengan berbagai ukuran yang digunakan seperti Moran’s I dan Geary’s C yang dapat digunakan untuk melihat pola cluster penyakit yang sedang diteliti dan juga dinas kesehatan dapat mengidentifikasi lokasi yang perlu ditangani terlebih dahulu dengan melihat Hotspot \(G_i^*\)

Sebutkan keterbatasan dari analisis autokorelasi spasial, misalnya terkait dengan :

  • MAUP (Modifiable Areal Unit Problem)

    MAUP dapat menjadi masalah pada hasil analisis spasial yang bergantung pada skala tertentu. Sebagai contoh, hasil autokorelasi spasial bisa memberikan hasil yang berbeda jika membandingkan hasil yang menggunakan data level kecamatan dan data level kelurahan

  • Ukuran Bobot Spasial (rook, queen, KNN)

    Pemilihan bobot spasial ini akan memengaruhi analisis spasial yang dihasilkan. Dalam praktek nya, Local Moran’s I dan Getis-Ord akan menghasilkan hasil yang berbeda bergantung pada pemilihan bobot spasial nya

  • Masalah multiple testing pada analisis lokal

    Pada analsis lokal seperti Local Moran’s I dan Getis-Ord pastinya akan dilakukan banyak uji signifikansi secara simultan untuk tiap wilayah yang diuji. Hal ini bisa menimbulkan masalah seperti false positive sehingga perlu kehati hatian dalam menggunakan analisis lokal agar mengurangi resiko terjadinya false positive yang dihasilkan