Autokorelasi Spasial Positif merupakan keadaan dimana unit-unit spasial cenderung memiliki atribut yang mirip satu dengan lainnya, yang berpotensi menghasilkan sebuah pola clustering antar data-data yang mirip, seperti geng ojek yang terbagi menjadi dua antara ojek pangkalan dengan ojek online/daring

Autokorelasi Spasial Positif merupakan keadaan dimana unit-unit spasial cenderung memiliki atribut yang berbeda satu dengan lainnya, yang berpotensi menghasilkan sebuah pola yang dispersing, mengakibatkan data saling mengelilingi dengan yang lain, seperti bidak-bidak catur, dimana buah catur seperti kuda dan peluncur (bishop) tidak saling berdekatan

Berikut adalah Rumus Matematis untuk Moran’s I, Geary’s C, Local Moran’s I, dan Getis-Ord Gi dan Gi*:

setwd("C:/Users/rasen/Downloads")
knitr::include_graphics("Spatial.jpg",error = FALSE)

Adapun beberapa perbedaan utama ukuran global dengan lokal pada autokorelasi spasial adalah sebagai berikut:

### Peta Kota Bandung
# Library
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(sp)
library(spData)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(spdep)
library(dplyr)
## 
## 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)

# Baca data GADM level 3 Indonesia (kecamatan)
Indo_Kec <- readRDS("gadm36_IDN_3_sp.rds")

# Ambil hanya Kota Bandung
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung", ]
Bandung$id <- seq_len(nrow(Bandung))

# Ubah ke format sf
Bandung_sf <- sf::st_as_sf(Bandung)

# ==== Opsi B: fallback (grid mock) bila Bandung_sf tidak tersedia ====
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
}

# Pastikan geometri valid
Bandung_sf <- sf::st_make_valid(Bandung_sf)

# Pastikan kolom id ada
Bandung_sf$id <- if (!"id" %in% names(Bandung_sf))
  seq_len(nrow(Bandung_sf)) else Bandung_sf$id

# ==== Plot peta sederhana ====
ggplot(Bandung_sf) +
  geom_sf(fill = "grey90", color = "white") +
  labs(title = "Kota Bandung - Unit (contoh)") +
  theme_minimal()

### Data yang Digunakan
# Memanggil Data 
data_csv <- read.csv("Luas Wilayah Kecamatan Kota Bandung.csv")

# Nama Kolom Plotting
data_csv <- data_csv %>%
  dplyr::rename(
    Luas_Km2 = "Luas.Wilayah..Km.persegi.",
    Kecamatan = "Kecamatan" # Kolom kunci
  )

# GABUNGKAN (JOIN) DATA ATRIBUT DENGAN DATA SPASIAL
Bandung_sf_LENGKAP <- Bandung_sf %>%
  dplyr::left_join(data_csv, by = c("NAME_3" = "Kecamatan")) 

# Periksa apakah proses join berhasil (opsional, tapi disarankan)
# print(head(Bandung_sf_LENGKAP))

ggplot(Bandung_sf_LENGKAP) +
  # aes(fill = Luas_Km2) mengisi warna berdasarkan kolom Luas_Km2
  geom_sf(aes(fill = Luas_Km2), color = "black", linewidth = 0.2) +
  
  # Skala warna (contoh menggunakan palet viridis)
  scale_fill_viridis_c(
    option = "plasma", # Pilih palet warna yang menarik
    name = expression("Luas Wilayah (km"^2*")") # Label legenda dengan notasi LateX
  ) +
  
  # Tambahkan Label
  labs(
    title = "Peta Luas Wilayah Kecamatan di Kota Bandung",
    caption = "Sumber Data: BPS Kota Bandung"
  ) +
  
  # Tema visual
  theme_minimal() +
  theme(
    axis.text = element_blank(), # Hilangkan koordinat sumbu
    axis.ticks = element_blank(),
    panel.grid.major = element_line(colour = "transparent"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

### Uji Moran's I

nb <- spdep::poly2nb(sf::as_Spatial(Bandung_sf), queen = TRUE)
lwW <- spdep::nb2listw(nb, style="W") # row-standardized
lwB <- spdep::nb2listw(nb, style="B") # binary (untuk raw G)

moran_res <- spdep::moran.test(data_csv$Luas_Km2, lwW, randomisation = TRUE,
                               alternative = "two.sided");moran_res
## 
##  Moran I test under randomisation
## 
## data:  data_csv$Luas_Km2  
## weights: lwW    
## 
## Moran I statistic standard deviate = -1.1532, p-value = 0.2488
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.16624656       -0.03448276        0.01305452

Hasil uji Moran’s I tersebut menunjukkan nilai statistik sebesar −0.16624656 dengan p-value 0.2488. Nilai negatif ini secara deskriptif mengarah pada adanya sedikit kecenderungan autokorelasi spasial negatif (dispersi), yang berarti unit-unit wilayah dengan luas besar cenderung berdekatan dengan unit-unit wilayah dengan luas kecil, dan sebaliknya. Namun, kecenderungan ini tidak terbukti signifikan secara statistik (dengan taraf signifikansi sebesar 0.05) karena p-value hasil uji lebih besar dari nilai taraf signifikansinya.

Oleh karena itu, dapat disimpulkan bahwa tidak ada bukti signifikan adanya ketergantungan spasial (spatial autocorrelation) pada luas wilayah kecamatan di Bandung berdasarkan matriks bobot yang digunakan. Distribusi luas wilayah unit-unit spasial tersebut dianggap acak (random). Artinya, dengan mengetahui seberapa besar luas suatu wilayah (dalam km persegi) di suatu kecamatan tidak memberikan informasi yang berarti untuk memprediksi apakah wilayah tetangganya juga akan berukuran besar atau kecil.

### Uji Geary's C

geary_res <- spdep::geary.test(data_csv$Luas_Km2, lwW, randomisation = TRUE,
                               alternative = "two.sided");geary_res
## 
##  Geary C test under randomisation
## 
## data:  data_csv$Luas_Km2 
## weights: lwW   
## 
## Geary C statistic standard deviate = -1.1357, p-value = 0.2561
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        1.13357119        1.00000000        0.01383218

Hasil Uji Geary’s C yang teramati adalah 1.133 dengan nilai p-value 0.2561. Nilai C lebih besar dari ekspektasi (1.0) menunjukkan adanya kecenderungan ke arah autokorelasi spasial negatif (dispersi), di mana wilayah yang luasnya sangat berbeda cenderung berdekatan. Namun, karena p-value ini jauh lebih tinggi dari 0.05, hasil ini tidak signifikan secara statistik. Artinya, pola dispersi yang teramati pada luas wilayah kecamatan di bandung adalah pola yang acak.

### Uji Local Moran's I (LISA)

# 1. Standardisasi Variabel
# PERBAIKAN: Menggunakan objek spasial gabungan (Bandung_data)
x <- scale(data_csv$Luas_Km2)[,1]

# 2. Hitung Nilai Lag
lagx <- spdep::lag.listw(lwW, x)

# 3. Hitung Local Moran's I (LISA)
lisa <- spdep::localmoran(x, lwW, 
                         alternative = "two.sided", # Kesalahan ketik ('alternavite') diperbaiki di sini
                         zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
alpha <- 0.05

# 4. Tentukan Kuadran (High-High, Low-Low, Outliers)
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)"
)

# 5. Gabungkan Hasil LISA ke Objek Spasial dan Tentukan Signifikansi
# PERBAIKAN: Menggabungkan hasil LISA ke objek SPASIAL (Bandung_data), bukan data_csv
Bandung_LISA <- dplyr::bind_cols(Bandung_sf, lisa_df) |>
  dplyr::mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))

# 6. Plot Peta LISA
ggplot(Bandung_LISA) +
  # Menggunakan linewidth untuk konsistensi
  geom_sf(aes(fill = quad), color="white", linewidth=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): Luas Wilayah", fill="Kategori") +
  theme_minimal()


Dari peta choropleth hasil uji Local Moran’s I (LISA) ini, terlihat bahwa sebagian besar wilayah kecamatan didominasi oleh kategori “Not significant” (Berwarna abu-abu), sehingga dapat disimpulkan bahwa luas suatu kecamatan tidak memiliki korelasi spasial yang signifikan dengan luas wilayah tetangganya, serta tidak terbukti adanya pengklasteran di mayoritas wilayah

Namun, berbeda dengan Kecamatan Babakan Ciparay dan Kecamatan Kiaracondong (warna biru) yang masuk dalam kategori Low-High. Ini berarti Kecamatan Babakan Ciparay dan Kecamatan Kiaracondong memiliki luas wilayah yang relatif rendah (kecil), tetapi mereka secara langsung dikelilingi oleh kecamatan-kecamatan tetangga yang memiliki luas wilayah yang jauh lebih tinggi (besar). Keberadaan outlier Low-High ini sangat penting bagi dinas setempat karena mengidentifikasi titik fokus heterogenitas spasial yang ekstrem.

x_raw <- data_csv$Luas_Km2
sum_x <- sum(x_raw)

Wb <- spdep::listw2mat(lwB)

num_G <- as.numeric(Wb %*% x_raw)
den_G <- (sum_x - x_raw)
G_raw <- num_G/den_G
Wb_star <- Wb; diag(Wb_star) <- 1
num_Gs <- as.numeric(Wb_star %*% x_raw) 
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) 
Bandung_G <- dplyr::mutate(Bandung_sf,
                           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.05792   Min.   :0.07883   Min.   :-1.505604   MULTIPOLYGON :30  
##  1st Qu.:0.11657   1st Qu.:0.14584   1st Qu.:-0.547311   epsg:NA      : 0  
##  Median :0.16397   Median :0.19328   Median : 0.002264   +proj=long...: 0  
##  Mean   :0.16430   Mean   :0.19211   Mean   : 0.133310                     
##  3rd Qu.:0.20991   3rd Qu.:0.22892   3rd Qu.: 0.944439                     
##  Max.   :0.27930   Max.   :0.30540   Max.   : 2.257632
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()


Dari Hasil peta Choropleth Uji Getis-Ord G-star, terlihat adanya variasi spasial dari warna-warna kecamatan dalam distribusi nilai G-star_raw. Semakin rendah nilai G-star_rownya, semakin ia jauh dari wilayah lain yang setipe dengan wilayar tersebut, dan berlaku juga sebaliknya.

Salah satu hal yang paling terlihat adalah kecamatan yang diwarnai dengan warna biru tua, seperti Kecamatan Bandung Kulon, Sukasari, dan Cibiru yang dimana 3 kecamatan tersebut tidak dikelilingi atau tidak dekat dengan kecamatan yang memiliki luas mirip dengan kecamatan tersebut. Di lain sisi, kecamatan Bandung Wetan terindikasi dekat dengan wilayah-wilayah lain yang memiliki luas kecamatan yang mirip dengan kecamatan lainnya, dilihat dari peta choropleth yang berwarna kuning

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


Dari Hasil Uji Getis-Ord Gi-star pada peta Choropleth menunjukkan bahwa hampir semua wilayah tidak menunjukkan adanya pola pengelompokan kuat antar wilayah yang seragam, kecuali pada Kecamatan Babakan Ciparay dan Kecamatan Kiaracondong yang terindikasi dikelilingi oleh wilayah dengan nilai variabel tinggi.

Hasil peta ini terbilang sangat mirip dengan hasil tes Local Moran’s I (LISA). Perbedaannya bahwa uji LISA mendeteksi bahwa 2 kecamatan tersebut merupakan kecamatan yang memiliki nilai variabel rendah, namun dikelilingi oleh kecamatan dengan nilai variabel yang tinggi, sedangkan uji ini mendeteksi 2 kecamatan tersebut memiliki nilai variabel tinggi yang dikelilingi oleh kecamatan yang cenderung luas.

Terlepas dari perbedaan tersebut, persamaan dari 2 uji tersebut adalah kedua kecamatan tersebut dikelilingi oleh kecamatan-kecamatan dengan luas kecamatan yang cenderung luas.

# Rook: buat ulang neighbor dengan queen = FALSE
nb_rook <- spdep::poly2nb(sf::as_Spatial(Bandung_sf), queen = FALSE)
lw_rook <- spdep::nb2listw(nb_rook, style="W")
# kNN (k = 4)
coords <- sf::st_coordinates(sf::st_centroid(Bandung_sf))
## Warning: st_centroid assumes attributes are constant over geometries
nb_knn <- spdep::knn2nb(spdep::knearneigh(coords, k = 4))
lw_knn <- spdep::nb2listw(nb_knn, style="W")

c(mean_z_queen = mean(spdep::localG(data_csv$Luas_Km2, lwW)),
  mean_z_rook = mean(spdep::localG(data_csv$Luas_Km2, lw_rook)),
  mean_z_knn4 = mean(spdep::localG(data_csv$Luas_Km2, lw_knn))
)
## mean_z_queen  mean_z_rook  mean_z_knn4 
##   0.13331031   0.09777828   0.25563958

Hasil analisis autokorelasi spasial luas wilayah ini dapat memberikan pandangan baru untuk pemerintah setempat, baik jika hasilnya signifikan ataupun acak. Jika analisis menunjukkan adanya klasterisasi (Hot Spot ataupun cold spot) luas wilayah, pemerintah dapat menggunakan informasi ini untuk alokasi sumber daya dan perencanaan infrastruktur yang lebih terfokus dan bisa menjadi acuan untuk perbandingan kualitas infrastruktur yang ideal satu kecamatan dengan kecamatan lainnya

Selain itu, apabila analisis menunjukkan bahwa luas wilayah bersifat acak (tidak signifikan - dilihat dari uji Moran’s I dan Geary’s C), dapat dikatakan bahwa luas wilayah kecamatan tidak ditentukan berdasarkan kecamatan-kecamatan yang ada di sekitarnya. Pembagian kecamatan tersebut mungkin berpatok kepada faktor-faktor lain, seperti ekonomi, infrastruktur, kepadatan penduduk, dll.

Meskipun begitu, analisis autokorelasi spasial tetap memiliki keterbatasan dari hasil ujinya, seperti pada aspek: