Autokorelasi Spasial, merupakan indikator yang melihat hubungan antara sebuah wilayah dengan wilayah lain yang berdekatan atau bertetanggaan. Autokorelasi spasial positif menunjukkan hubungan antara wilayah lain yang memiliki hubungan positif, atau sebuah hubungan yang memiliki karakteristik yang sama atau berdekatan. Walaupun autokorelasi spasial negatif menunjukkan sebaliknya, dimana sebuah wilayah dengan wilayah berdekatan memiliki karakteristik yang berbeda atau tidak berdekatan.
\[ I = \frac{n}{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij}} \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} \]
\[ C = \frac{(n - 1) \sum_{i} \sum_{j} w_{ij} (x_i - x_j)^2} {2n S^2 \sum_{i} \sum_{j} w_{ij}} \]
dengan varians sampel: \[ S^2 = \frac{1}{n} \sum_{i} (x_i - \bar{x})^2 \] - Local Moran’s I
\[ I_i = \frac{(x_i - \bar{x}) \sum_{j=1}^{n} w_{ij} (x_j - \bar{x})} {\dfrac{\sum_{i=1}^{n} (x_i - \bar{x})^2}{n}} \]
Tanpa memasukkan lokasi i \[ G_i(d) = \frac{\sum_{j \ne i} w_{ij}(d)\, x_j} {\sum_{j \ne i} x_j} \]
Dengan memasukkan lokasi i \[ G_i^*(d) = \frac{\sum_{j=1}^{N} w_{ij}(d)\, x_j} {\sum_{j=1}^{N} x_j} \]
Standarisasi ke Z-skor \[ z(G_i^*) = \frac{\sum_{j} w_{ij} x_j - \bar{x} \sum_{j} w_{ij}} {\sqrt{\frac{\sum_{k} (x_k - \bar{x})^2}{N}} \; \sqrt{\frac{N \sum_{j} w_{ij}^2 - \left(\sum_{j} w_{ij}\right)^2}{N - 1}} } \]
Global mengukur autokorelasi spasial pada seluruh area yang ingin diteliti kemudian menunjukkan apa masing-masing area memiliki kesamaan karakteristik.
Lokal mengukur dimana terjadinya klustering seperti adanya hotspot atau coldspot pada area lokal dan tidak melihat semua area penelitian sekaligus.
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')`
## Loading required package: 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(sf)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
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
setwd("C:\\Users\\HP\\Documents\\KULIAH\\SEMESTER 5\\Spatial")
Indo_Kec <- readRDS('gadm36_IDN_3_sp.rds')
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung",]
plot(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.
set.seed(123)
n_kec <- length(unique(Bandung$NAME_3))
z <- rnorm(n_kec, mean = 50, sd = 10)
lz <- rnorm(n_kec, mean = 5, sd = 2)
beta0 <- 1.2
beta1 <- 0.6
beta2 <- -0.8
lambda <- exp(beta0 + beta1*scale(z) + beta2*scale(lz))
pop <- round(runif(n_kec, 3000, 8000))
kasus <- rpois(n_kec, lambda = lambda * (pop/10000))
data_kasus2 <- data.frame(
NAME_3 = unique(Bandung$NAME_3),
pop = pop,
kepadatan= z,
sanitasi = lz,
lambda = lambda,
kasus = kasus,
rate10k = (kasus/pop)*10000
)
print(data_kasus2)
## NAME_3 pop kepadatan sanitasi lambda kasus rate10k
## 1 Andir 6239 44.39524 5.852928 1.9123482 2 3.205642
## 2 Antapani 4599 47.69823 4.409857 4.6717261 2 4.348771
## 3 Arcamanik 4539 65.58708 6.790251 4.4614782 1 2.203128
## 4 Astanaanyar 4099 50.70508 6.756267 1.8249855 0 0.000000
## 5 Babakan Ciparay 4847 51.29288 6.643162 1.9970982 0 0.000000
## 6 Bandung Kidul 7921 67.15065 6.377281 5.9829166 4 5.049867
## 7 Bandung Kulon 3771 54.60916 6.107835 3.1611140 1 2.651816
## 8 Bandung Wetan 3455 37.34939 4.876177 1.9842473 0 0.000000
## 9 Batununggal 3710 43.13147 4.388075 3.5703309 1 2.695418
## 10 Bojongloa Kaler 6450 45.54338 4.239058 4.4439576 1 1.550388
## 11 Bojongloa Kidul 6096 62.24082 3.610586 16.6728504 10 16.404199
## 12 Buahbatu 7457 53.59814 4.584165 6.1649501 5 6.705109
## 13 Cibeunying Kaler 6365 54.00771 2.469207 17.4079713 10 15.710919
## 14 Cibeunying Kidul 6685 51.10683 9.337912 0.5431423 0 0.000000
## 15 Cibiru 5606 44.44159 7.415924 0.9071392 1 1.783803
## 16 Cicendo 6299 67.86913 2.753783 35.4595296 18 28.575964
## 17 Cidadap 7109 54.97850 4.194230 8.0854823 6 8.440006
## 18 Cinambo 6931 30.33383 4.066689 1.9038646 0 0.000000
## 19 Coblong 7899 57.01356 6.559930 2.9489187 4 5.063932
## 20 Gedebage 5197 45.27209 4.833262 3.2882252 2 3.848374
## 21 Kiaracondong 4559 39.32176 5.506637 1.6551633 1 2.193463
## 22 Lengkong 5047 47.82025 4.942906 3.6461719 2 3.962750
## 23 Mandalajati 3052 39.73996 4.914259 2.2551282 0 0.000000
## 24 Panyileukan 3919 42.71109 7.737205 0.6996463 0 0.000000
## 25 Rancasari 7214 43.74961 4.548458 3.4337227 4 5.544774
## 26 Regol 4156 33.13307 8.032941 0.3380295 0 0.000000
## 27 Sukajadi 4195 58.37787 1.902494 29.8339001 16 38.140644
## 28 Sukasari 3383 51.53373 6.169227 2.5432065 1 2.955956
## 29 Sumur Bandung 4229 38.61863 5.247708 1.7948378 0 0.000000
## 30 Ujung Berung 6661 62.53815 5.431883 7.0966474 5 7.506380
Bandung_sf <- st_as_sf(Bandung)
Bandung_sf$NAME_3 <- as.character(Bandung_sf$NAME_3)
data_kasus2$NAME_3 <- as.character(data_kasus2$NAME_3)
Bandung_merged <- Bandung_sf %>%
left_join(data_kasus2, by = "NAME_3")
names(Bandung_merged)
## [1] "GID_0" "NAME_0" "GID_1" "NAME_1" "NL_NAME_1" "GID_2"
## [7] "NAME_2" "NL_NAME_2" "GID_3" "NAME_3" "VARNAME_3" "NL_NAME_3"
## [13] "TYPE_3" "ENGTYPE_3" "CC_3" "HASC_3" "pop" "kepadatan"
## [19] "sanitasi" "lambda" "kasus" "rate10k" "geometry"
ggplot(Bandung_merged) +
geom_sf(aes(fill = rate10k), color = "black") +
scale_fill_viridis_c(option = "plasma") +
labs(
title = "Simulasi Kasus Diare per 10.000 Penduduk\nKota Bandung",
fill = "Kasus/10k"
) +
theme_minimal()
Bisa dilihat di peta, bahwa kasus diare di Kota Bandung terpusat pada kecamatan Sukajadi dan Cicendo, dengan ada peningkatan pada Cibeunyi Kaler dan Bojongloa Kidul. Menariknya kecamatan sekitar Cicendo dan Sukajadi tidak memiliki kasus diare yang begitu signifikan yang mengindikasi kurangnya pengaruh spasial.
row.names(Bandung) <- as.character(1:30)
W <- poly2nb(Bandung, row.names(Bandung), queen=FALSE)
WB <- nb2listw(W, style='B', zero.policy = TRUE)
WL <- nb2listw(W, style='W')
CoordK <- coordinates(Bandung)
plot(Bandung, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Bandung), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7, col="blue")
plot.nb(W, CoordK, add=TRUE, col="red", lwd=2)
moran_res <- spdep::moran.test(Bandung_merged$rate10k, WL,
randomisation = TRUE, alternative = "two.sided")
moran_res
##
## Moran I test under randomisation
##
## data: Bandung_merged$rate10k
## weights: WL
##
## Moran I statistic standard deviate = 1.5937, p-value = 0.111
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.12830192 -0.03448276 0.01043287
Moran I mendapatkan nilai sebesar 0.12 dimana jika memang dibulatkan mendapatkan 0. Dan dari hasil analisis yang dilakukan, menunjukkan bahwa Moran I tidak signifikan secara statistik. Dari hal tersebut bisa disimoulkan bahwa kasus diare terutama di Kota Bandung tidak mengikuti sebuah pola tertentu melainkan mengikuti pola acak.
geary_res <- spdep::geary.test(Bandung_merged$rate10k, WL,
randomisation = TRUE, alternative = "two.sided")
geary_res
##
## Geary C test under randomisation
##
## data: Bandung_merged$rate10k
## weights: WL
##
## Geary C statistic standard deviate = 0.89647, p-value = 0.37
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic Expectation Variance
## 0.88299530 1.00000000 0.01703457
Didapatkan Geary C sebesar 0.8 dimana nilai sudah mendekati 1, serta juga dengan o-value 0.37 dibawah taraf signifikan 5% kita bisa menerima hipotesis null sehingga menunjukkan bahwa data penyakit diare tidak mengandung autokorelasi spasial. Dibandingkan dengan Moran I bisa dilihat bahwa keduanya mendapatkan kesimpulan yang sama, dengan kedua statistik mendekati taraf yang menunjukkan tidak adanya autokorelasi dan mengindikasi data acak. Keduanya memiliki sensitivitas yang berbeda Geary C lebih sensitif pada perbedaan nilai antar tetangga walaupun Moran I lebih sensitif terhadap kovarians antar lokasi.
x <- scale(Bandung_merged$rate10k)[,1]
lagx <- spdep::lag.listw(WL, x)
lisa <- spdep::localmoran(x, WL, 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()
Didapatkan Kecamatan: a. High-High: Kecamatan Cicadap, Kecamatan
Sukajadi b. Low-High: Kecamatan Coblong, Kecamatan Sukasari
Dilihat bahwa Kecamatan Cicadap dan Sukajadi merupakan kecamatan yang memiliki rawan diare yang tinggi, mungkin karena fasilitas sanitasi pada kedua kecamatan tersebut ataupunf faktor-faktor lainnya, karena itu bisa disebutkan sebuah hotspot. Namun disekitar kedua kecamatan tersebut ada yang memiliki kasus yang rendah yaitu Coblong dan Sukasari walaupun berdekatan dengan kecamatan yang tinggi, mengindikasi perlindungan yang bagus seperti fasilitas sanitasi yang signifikan lebih bagus dari tetanggany ataupun ada kasus-kasus yang tidak dicatat. Walaupun itu, di kecamatan lain tidak ada indikasi kasus diare yang signifikan.
x_raw <- Bandung_merged$rate10k
sum_x <- sum(x_raw)
Wb <- spdep::listw2mat(WB)
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
G_star_raw <- num_Gs / den_Gs
Gz <- spdep::localG(x_raw, listw = WL)
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.01932 Min. :0.03475 Min. :-1.321986 MULTIPOLYGON :30
## 1st Qu.:0.05701 1st Qu.:0.08325 1st Qu.:-0.748370 epsg:NA : 0
## Median :0.11271 Median :0.12441 Median :-0.533315 +proj=long...: 0
## Mean :0.14180 Mean :0.16861 Mean :-0.049237
## 3rd Qu.:0.19421 3rd Qu.:0.21177 3rd Qu.: 0.007445
## Max. :0.38104 Max. :0.49351 Max. : 2.916574
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()
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()
Kecamatan Hot-Spot bisa dilihat pada keempat kecamatan di Kota Bandung
utara yaitu Kecamatan Cicadap, Kecamatan Sukajadi, Kecamatan Coblong,
dan Kecamatan Sukasari menunjukkan besarnya atau signifikannya kasus
diare pada keempat kecamatan tersebut yang memusat pada Kecamatan
Sukajadi. Dibanding dengan peta LISA bisa dilihat bahwa keempat
kecamatan yang sama yang dihighlight dimana keempatnya adalah hotspot
kecamatan dengan kasus diare yang tinggi dan dikelilingi kecamatan yang
tinggi, namun pada LISA dilihat bahwa Sukasari dan Coblong keduanya
dikatakan Low-High menunjukkan bahwa keduanya memiliki tingkat yang
rendah namun dikelilingi kecamatan yang tinggi, menunjukkan bahwa yang
dikatakan rendah di LISA bisa dikatakan tinggi di Getis.
dengan:
MAUP (Modifiable Areal Unit Problem) Hasil analisis bergantung pada pembagian wilayah seperti kecamatan sehingga akan ada masalah dalam merepresentasikan seluruh wilayah tersebut terutama jika berada di perbatasan area. Dan hasil autokorelasi spasial sangatlah bergantung pada hal tersebut sehingga perbedaan area akan mendapatkan hasil yang sangat berbeda.
Ukuran Bobot Spasial Hasil autokorelasi sangat bergantung pada pembuatan matriks bobot spasial (rook, queen, k-nearest neighbors) sehingga jika ada pemilihan matriks bobot yang berbeda akan memberikan hasil serta kesimpulan yang berbeda. Terutana saat pemilihan bobot matriks sangatlah subyektif.
Masalah Multiple Testing Melakukan test yang berbeda serta estimasi berbeda pada masing-masing area bisa menyebabkan overestimation menyebabkan errornya estimasi serta terjadi adanya type-1 atau type-2.