knitr::opts_chunk$set(warning = FALSE, message = FALSE)Praktikum 6 : Bivariate Spatial Autocorrelation
Bivariate Spatial Autocorrelation adalah metode statistik yang digunakan untuk mengukur hubungan spasial antara dua peubah di area geografis tertentu. Ini mengevaluasi apakah pola sebaran satu peubah secara spasial berkaitan dengan pola sebaran peubah lainnya.
Beberapa metode yang umum digunakan untuk mengukur bivariate spatial autocorrelation adalah Moran’s I bivariate, Geary’s C bivariate dan Lee’s Spatial Correlation.
Berikut adalah beberapa paket yang perlu di install terlebih dahulu.
# Memuat paket yang diperlukan
library(spdep)
library(sf)
library(bispdep)
library(tmap)
library(ggplot2)
library(dplyr)Kerangka data Columbus memiliki 49 baris dan 22 kolom. Unit analisis: 49 wilayah di Columbus, OH, data tahun 1980. Kerangka data ini berisi kolom-kolom berikut:
X : koordinat x
Y : koordinat y
HOVAL : nilai rumah (dalam “$1.000”)
INC : pendapatan rumah tangga (dalam “$ 1.000”)
CRIME : pencurian rumah dan pencurian kendaraan per seribu rumah tangga di lingkungan tersebut
OPEN : ruang terbuka di wilayah sekitar
PLUMB : persentase unit rumah tanpa pipa ledeng
DISCBD : jarak ke pusat bisnis (CBD)
columbus: kerangka data, berisi 22 variabel
col.gal.nb: adalah objek dari kelas “nb”, daftar vektor, satu untuk setiap unit spasial, dan berisi nomor urut tetangga (file contiguity ini menggunakan definisi queen untuk Columbus)
coords: koordinat centroid yang dapat digunakan untuk membangun bobot berbasis jarak
# Memuat data columbus
columbus <- st_read(system.file("shapes/columbus.gpkg", package = "spData"), quiet = TRUE)plot(columbus)
tm_shape(columbus) + tm_polygons("INC") + tm_layout(title = "Map of INCOME")
tm_shape(columbus) + tm_polygons("HOVAL") + tm_layout(title = "Map of HOVAL")
Matriks pembobot spasial yang akan digunakan dalam praktikum ini adalah Queen Contiguity (Silahkan jika ingin mencoba menggunakan matriks pembobot yang lain).
queen <- poly2nb(columbus, queen = TRUE)
W.queen <- nb2listw(queen, style='W')
summary(queen)Neighbour list object:
Number of regions: 49
Number of nonzero links: 236
Percentage nonzero weights: 9.829238
Average number of links: 4.816327
Link number distribution:
2 3 4 5 6 7 8 9 10
5 9 12 5 9 3 4 1 1
5 least connected regions:
1 6 42 46 47 with 2 links
1 most connected region:
20 with 10 linksBivariate Geary’s Cxy
Fungsi sederhana untuk menghitung Bivariate Geary’s Cxy adalah sebagai berikut:
Fungsi yang dapat digunakan adalah geary.bi dari paket bispdep, fungsi ini digunakan untuk menghitung koefisien Geary bivariate (Cxy), yang merupakan ukuran ketergantungan spasial.
Sebagai contoh, kita ingin melihat apakah terdapat autokorelasi spasial antara peubah crime dan income.
geary_bi <- geary.bi(columbus$HOVAL, columbus$INC, W.queen, zero.policy=TRUE, alternative = "greater")
geary_bi$C
[1] 1.303785
$Kx
[1] 4.312181
$Ky
[1] 3.771052Geary’s C memiliki skala yang berbeda dibandingkan Moran’s I. Nilai Geary’s C berkisar dari 0 hingga 2:
Jika C = 1, ini menunjukkan tidak ada autokorelasi spasial (random).
Jika C < 1, ini menunjukkan autokorelasi positif (wilayah yang berdekatan cenderung memiliki nilai yang mirip).
Jika C > 1, ini menunjukkan autokorelasi negatif (wilayah yang berdekatan cenderung memiliki nilai yang berbeda atau berlawanan).
Geary’s Cxy > 1 menunjukkan autokorelasi spasial negatif antara variabel CRIME dan INC. Dengan kata lain, wilayah yang berdekatan cenderung memiliki nilai CRIME dan INC yang berlawanan (jika CRIMEtinggi, INC rendah, atau sebaliknya)
Kurtosis menunjukkan bahwa sebaran CRIME mendekati distribusi normal, sedangkan INC memiliki sebaran yang lebih ekstrem.
Uji Bivariate Geary’s C untuk Autokorelasi Spasial
Fungsi yang dapat digunakan adalah gearybi.test.
# Uji Geary's C Bivariate
gearybi_test <- gearybi.test(columbus$HOVAL, columbus$INC, W.queen, randomisation = TRUE, alternative = "greater")
gearybi_test
Geary C_{xy} test under randomisation
data: columbus$HOVAL
weights: W.queen
Zc_{xy} Statistics of Geary's C_{xy} = -2.926, p-value = 0.9983
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C_{xy} statistic Expectation Variance
1.30378453 1.00000000 0.01077903 alternative="greater", uji ini bertujuan untuk menguji apakah ada autokorelasi spasial positif, alternative="less": menguji apakah ada autokorelasi spasial negatif (lokasi yang bertetangga memiliki nilai yang sangat berbeda),
Berdasarkan output di atas, diperoleh p-value 0.9983 > 0.05, artinya gagal Tolak H0 yang menyatakan bahwa tidak terdapat autokorelasi spasial positif antar peubah HOVAL dan INCOME pada taraf nyata 5%.
# Uji Geary's C Bivariate
gearybi_test <- gearybi.test(columbus$HOVAL, columbus$INC, W.queen, randomisation = TRUE, alternative = "less")
gearybi_test
Geary C_{xy} test under randomisation
data: columbus$HOVAL
weights: W.queen
Zc_{xy} Statistics of Geary's C_{xy} = -2.926, p-value = 0.001717
alternative hypothesis: Expectation less than statistic
sample estimates:
Geary C_{xy} statistic Expectation Variance
1.30378453 1.00000000 0.01077903 Berdasarkan output di atas, diperoleh p-value 0.001717 < 0.05, artinya Tolak H0 yang menyatakan bahwa terdapat autokorelasi spasial negatif antar peubah HOVAL dan INCOME pada taraf nyata 5%. Dengan kata lain, wilayah yang berdekatan cenderung memiliki nilai CRIME dan INC yang berlawanan (jika CRIMEtinggi, INC rendah, atau sebaliknya).
Bivariate Moran’s Ixy
Menghitung Bivariate Moran’s Ixy
Fungsi moran.bi digunakan untuk menghitung koefisien Moran bivariate (Ixy), yang merupakan ukuran autokorelasi spasial antara dua peubah.
moran_bi <- moran.bi(columbus$HOVAL, columbus$INC, W.queen, zero.policy=TRUE)
moran_bi$I
[1] 0.2208649
$Kx
[1] 4.312181
$Ky
[1] 3.771052Moran’s I biasanya berkisar antara -1 dan 1:
Nilai positif menunjukkan autokorelasi spasial positif (wilayah yang berdekatan cenderung memiliki nilai variabel yang serupa).
Nilai negatif menunjukkan autokorelasi spasial negatif (wilayah yang berdekatan cenderung memiliki nilai variabel yang berlawanan).
Nilai mendekati 0 menunjukkan tidak ada autokorelasi
Nilai 0.2208649 menunjukkan autokorelasi spasial positif yang lemah antara nilai properti rumah (HOVAL) dan pendapatan (INC). Ini berarti wilayah dengan nilai properti rumah yang tinggi cenderung berdekatan dengan wilayah yang juga memiliki pendapatan tinggi, meskipun hubungannya tidak terlalu kuat.
Kurtosis: Kedua peubah (HOVAL dan INC) memiliki sebaran yang lebih ekstrim (kurtosis lebih tinggi dari 3).
Uji Bivariate Moran’s Ixy untuk Autokorelasi Spasial
Fungsi moranbi.test digunakan untuk melakukan uji Moran bivariate untuk mendeteksi autokorelasi spasial antara dua peubah menggunakan matriks bobot spasial.
moranbi_test <- moranbi.test(columbus$HOVAL, columbus$INC, listw=W.queen, randomisation = TRUE, zero.policy=TRUE)
moranbi_test
Bivariate Moran I_{xy} test under randomisation
data: columbus$HOVAL
weights: W.queen
Bivariate Moran Z(I_{xy}) statistic = 3.0912, p-value = 0.0009969
alternative hypothesis: greater
sample estimates:
Bivariate Moran I_{xy} statistic Expectation
0.220864911 -0.010414137
Variance
0.005597975 Berdasarkan output diatas diperoleh p-value 0.009969 < 0.05 artinya Tolak H0 yang menyatakan bahwa terdapat autokorelasi spasial positif antara peubah HOVAL dan INC. Artinya, wilayah-wilayah dengan nilai properti yang tinggi cenderung berdekatan dengan wilayah-wilayah dengan pendapatan yang tinggi.
Bivariate Moran’s Ixy (Simulasi Monte Carlo)
Fungsi moranbi.mc melakukan uji Moran Bivariate (Ixy) menggunakan simulasi Monte Carlo untuk menghitung autokorelasi spasial antara dua variabel.
set.seed(123)
moran_mc <- moranbi.mc(columbus$HOVAL, columbus$INC, W.queen, nsim=999, zero.policy=TRUE, alternative = "greater")
moran_mc
Monte-Carlo simulation of Bivariate Moran I
data: columbus$HOVAL
weights: W.queen
number of simulations + 1: 1000
statistic = 0.22086, observed rank = 999, p-value = 0.001
alternative hypothesis: greaterBerdasarkan output simulasi monte carlo diatas diperoleh p-value 0.001 < 0.05 artinya Tolak H0 yang menyatakan bahwa terdapat autokorelasi spasial positif antara peubah HOVAL dan INC. Artinya, wilayah-wilayah dengan nilai properti yang tinggi cenderung berdekatan dengan wilayah-wilayah dengan pendapatan yang tinggi.
Bivariate Moran Scatterplot
Fungsi moranbi.plot membuat scatterplot Moran bivariate, yang menampilkan hubungan antara variabel Y dan variabel X yang di-lag secara spasial.
HOVAL <- as.vector(scale(columbus$HOVAL))
INCOME <- as.vector(scale(columbus$INC))
moranbi.plot(HOVAL, INCOME, W.queen, label=columbus$POLYID, zero.policy=TRUE, plot=TRUE)
Secara umum, tren garis diagonal ini mengindikasikan adanya hubungan spasial positif antara HOVAL dan INCOME. Artinya, wilayah dengan nilai properti rumah yang lebih tinggi cenderung berada di dekat wilayah tetangga dengan pendapatan lebih tinggi, yang sejalan dengan hasil Bivariate Moran’s I yang menunjukkan autokorelasi spasial positif.
BiLISA- Bivariate Local Indicators of Spatial Association Moran’s Ixy statistic
Bivariate local spatial statistic Moran’s I dihitung untuk setiap zona berdasarkan objek pembobot spasial yang digunakan. Nilai yang dikembalikan termasuk nilai Z, dan dapat digunakan sebagai alat diagnostik. Statistik tersebut adalah
localmoran_bi <- localmoran.bi(columbus$HOVAL, columbus$INC, W.queen, zero.policy=TRUE, alternative = "two.sided")
localmoran_bi Ixyi E.Ixyi Var.Ixyi Z.Ixyi Pr(z != 0)
1 0.5308712143 -5.507621e-02 2.3508675675 0.382159419 0.7023431
2 -0.0193374585 -1.171818e-03 0.0365734440 -0.094987791 0.9243245
3 0.0091157320 -4.554183e-03 0.1035290160 0.042484917 0.9661121
4 -0.0089242936 -8.548020e-04 0.0195771127 -0.057672908 0.9540092
5 0.1403099306 -7.213691e-03 0.0741392753 0.541797999 0.5879577
6 -0.0012776656 -2.925083e-03 0.1394925505 0.004410910 0.9964806
7 -0.4543359635 -4.168034e-02 0.8764891380 -0.440772508 0.6593777
8 0.0229103060 -5.360228e-05 0.0007824714 0.820940301 0.4116803
9 -0.0260904517 -6.254421e-03 0.0644054663 -0.078161593 0.9376995
10 0.5493480801 -1.047472e-01 1.8995288790 0.474589455 0.6350796
11 0.2438773057 -1.094443e-02 0.1920044546 0.581541649 0.5608755
12 0.2947586952 -1.071203e-02 0.1530367458 0.780857406 0.4348864
13 -0.0432690112 -3.321009e-04 0.0076139109 -0.492070009 0.6226699
14 -0.0574317217 -6.212040e-04 0.0090578679 -0.596919484 0.5505611
15 0.2628033875 -1.302059e-02 0.1851399081 0.641035362 0.5214997
16 0.2242230919 -1.202113e-02 0.1223424963 0.675417831 0.4994103
17 0.0750583254 -3.423540e-04 0.0107029326 0.728826011 0.4661081
18 -0.2268750273 -1.449702e-02 0.3229415542 -0.373721061 0.7086119
19 0.0448916092 -1.914445e-03 0.0596624975 0.191624625 0.8480363
20 0.1593234253 -5.719276e-02 0.4014111351 0.341739870 0.7325467
21 0.0006413976 -1.062552e-02 0.3253450363 0.019752998 0.9842404
22 -0.0110271416 -1.988438e-03 0.0289143114 -0.053155684 0.9576079
23 0.1547918404 -2.694598e-03 0.0838439314 0.543885331 0.5865204
24 -0.1672807011 -6.795542e-03 0.0818840237 -0.560835233 0.5749099
25 0.2086558265 -1.314832e-02 0.1335051205 0.607044749 0.5438212
26 0.1998352419 -1.025470e-02 0.1466401218 0.548629122 0.5832600
27 -0.0099087650 -5.862098e-04 0.0134329036 -0.080435916 0.9358906
28 0.1594766108 -7.573751e-03 0.0674118764 0.643397269 0.5199664
29 0.0506645413 -1.098626e-03 0.0133910246 0.447315490 0.6546473
30 0.1223681377 -7.917719e-03 0.1397649428 0.348496280 0.7274675
31 -0.0598597581 -1.373003e-03 0.0428353102 -0.282589857 0.7774913
32 0.0652322752 -1.082969e-04 0.0024839808 1.311018493 0.1898515
33 -0.0562098563 -6.862399e-03 0.1552740994 -0.125231986 0.9003399
34 -0.1099204578 -3.109079e-03 0.0708840528 -0.401183674 0.6882849
35 -0.0953750334 -4.077505e-03 0.0494034175 -0.410752791 0.6812538
36 -0.0290178941 -1.422733e-04 0.0025511260 -0.571696298 0.5675277
37 -0.0417031007 -7.375238e-04 0.0107514384 -0.395080676 0.6927833
38 0.0749087147 -7.720230e-03 0.1109691595 0.248045330 0.8040993
39 0.0126330880 -4.222466e-05 0.0013208528 0.348763857 0.7272666
40 0.4520387910 -1.723749e-02 0.2985148147 0.858906040 0.3903924
41 0.0899397404 -4.184913e-04 0.0130812023 0.790030563 0.4295099
42 0.0435777111 -1.084073e-03 0.0518892102 0.196063853 0.8445602
43 0.0152680188 -5.057205e-03 0.0730846553 0.075183446 0.9400688
44 0.0117176483 -7.596595e-04 0.0136047528 0.106973341 0.9148101
45 -0.0280868356 -3.571564e-03 0.0813524688 -0.085951080 0.9315053
46 0.4772920969 -4.422592e-02 1.9337831731 0.375029495 0.7076385
47 0.1416698452 -5.148594e-04 0.0246719011 0.905214647 0.3653517
48 -0.0502692961 -4.221370e-03 0.0960277467 -0.148597569 0.8818712
49 0.0005677425 -2.166678e-04 0.0067753382 0.009529668 0.9923965
attr(,"call")
localmoran.bi(varX = columbus$HOVAL, varY = columbus$INC, listw = W.queen,
zero.policy = TRUE, alternative = "two.sided")
attr(,"class")
[1] "localmoran.bi" "matrix" "array" st_crs(columbus)Coordinate Reference System:
User input: Undefined Cartesian SRS with unknown unit
wkt:
ENGCRS["Undefined Cartesian SRS with unknown unit",
EDATUM["Unknown engineering datum"],
CS[Cartesian,2],
AXIS["x",unspecified,
ORDER[1],
LENGTHUNIT["unknown",0]],
AXIS["y",unspecified,
ORDER[2],
LENGTHUNIT["unknown",0]]]columbus <- st_set_crs(columbus, 4326)columbus$lmI <- localmoran_bi[, "Ixyi"] # local Moran's I
columbus$lmZ <- localmoran_bi[, "Z.Ixyi"] # z-scores
# p-values corresponding to alternative greater
columbus$lmp <- localmoran_bi[, "Pr(z != 0)"]p1 <- tm_shape(columbus) +
tm_polygons(col = "lmI", title = "Local Moran's I",
style = "quantile") +
tm_layout(legend.outside = TRUE)
p1
p2 <- tm_shape(columbus) +
tm_polygons(col = "lmp", title = "p-value",
breaks = c(-Inf, 0.05, Inf)) +
tm_layout(legend.outside = TRUE)
p2
tmap_mode("plot")tm_shape(columbus) + tm_polygons(col = "lmZ",
title = "Local Moran's I", style = "fixed",
breaks = c(-Inf, -1.96, 1.96, Inf),
labels = c("Negative SAC", "No SAC", "Positive SAC"),
palette = c("blue", "white", "red")) +
tm_layout(legend.outside = TRUE)
Dalam uji dua sisi ini, nilai z-score yang lebih rendah dari –1,96 menunjukkan autokorelasi spasial negatif, dan nilai z-score yang lebih besar dari 1,96 menunjukkan autokorelasi spasial positif.
Berdasarkan output di atas, dapat dilihat bahwa tidak terdapat autokorelasi spasial lokal antar peubah HOVAL dan INCOME di semua lokasi.
Lee’s Spatial Correlation (Lee’s L)
Lee’s L adalah koefisien korelasi spasial bivariat yang digunakan untuk mengukur hubungan antara dua set observasi yang dilakukan di lokasi spasial yang sama. Berbeda dengan ukuran asosiasi standar seperti koefisien korelasi Pearson, yang tidak mempertimbangkan dimensi spasial data, Lee’s L dirancang untuk mengatasi masalah tersebut dengan mengurangi potensi pembesaran asosiasi yang disebabkan oleh autokorelasi spasial. Lee’s L dapat ditemukan di berbagai pustaka perangkat lunak analisis spasial seperti spdep di R. Metode ini telah diterapkan dalam berbagai bidang, seperti penelitian polusi udara, vitikultura, dan analisis harga sewa rumah.
Untuk data spasial x_i dan y_i yang diukur di N lokasi yang terhubung dengan matriks bobot spasial w_ij, pertama tentukan vektor lag spasial:
Untuk data spasial x_i dan y_i yang diukur di N lokasi yang terhubung dengan matriks bobot spasial w_ij, pertama tentukan vektor lag spasial:
dengan definisi serupa untuk y_i.
Kemudian, Lee’s L didefinisikan sebagai:
di mana x̄ dan ȳ adalah nilai rata-rata dari x_i dan y_i. Ketika matriks bobot spasial dinormalisasi dengan cara baris (sehingga ∑_j w_ij = 1), faktor pertama adalah 1. Selengkapnya ada di “https://en.wikipedia.org/wiki/Lee%27s_L”
library(spData)
library(spdep)
data(oldcol)
col.W <- nb2listw(COL.nb, style="W")
crime <- COL.OLD$CRIME
lee.test(crime, crime, col.W, zero.policy=TRUE)
Lee's L statistic randomisation
data: crime , crime
weights: col.W
Lee's L statistic standard deviate = 5.2343, p-value = 8.279e-08
alternative hypothesis: greater
sample estimates:
Lee's L statistic Expectation Variance
0.547064219 0.239417989 0.003454459 Hasil uji Lee’s L statistic menunjukkan bahwa terdapat korelasi spasial positif yang signifikan antara variabel crime dan crime itu sendiri di wilayah yang dianalisis. Nilai Lee’s L statistic yang diperoleh adalah 0.547064219, yang mengindikasikan asosiasi spasial positif antara data yang dianalisis. Artinya, wilayah dengan tingkat kriminalitas tinggi cenderung memiliki tingkat kriminalitas tinggi di wilayah sekitarnya, dan begitu juga sebaliknya.
Dengan p-value yang sangat kecil (8.279e-08), kita dapat menolak hipotesis nol yang menyatakan tidak ada hubungan spasial, dan menyimpulkan bahwa hasil ini sangat signifikan secara statistik. Ini berarti bahwa ada bukti yang kuat bahwa data kriminalitas di wilayah yang berdekatan menunjukkan asosiasi spasial positif yang konsisten.
Secara keseluruhan, hasil ini menunjukkan bahwa ada hubungan spasial yang signifikan dan kuatantara variabel kriminalitas di wilayah yang dianalisis, dengan pola yang menunjukkan bahwa wilayah dengan tingkat kejahatan tinggi cenderung terhubung dengan wilayah lainnya yang memiliki tingkat kejahatan tinggi juga.
# Menggunakan Data Lain
data(boston, package="spData")lw<-nb2listw(boston.soi)
x<-boston.c$CMEDV
y<-boston.c$CRIM
lee.test(x, y, lw, zero.policy=TRUE, alternative="less")
Lee's L statistic randomisation
data: x , y
weights: lw
Lee's L statistic standard deviate = -11.54, p-value < 2.2e-16
alternative hypothesis: less
sample estimates:
Lee's L statistic Expectation Variance
-0.326297206 -0.105040316 0.000367637 Hasil uji Lee’s L statistic menunjukkan adanya korelasi spasial negatif yang signifikan antara cMEDV(nilai median rumah) dan cCRIM (tingkat kriminalitas) di wilayah yang dianalisis. Nilai Lee’s L statisticyang diperoleh adalah -0.3262972060, yang mengindikasikan adanya asosiasi spasial negatif antara kedua variabel. Hal ini berarti bahwa wilayah dengan tingkat kriminalitas tinggi cenderung memiliki nilai median rumah yang lebih rendah, dan sebaliknya.
Dengan p-value yang sangat kecil (lebih kecil dari 2.2e-16), kita dapat menolak hipotesis nol yang menyatakan bahwa tidak ada hubungan spasial, dan menyimpulkan bahwa hasil ini sangat signifikan secara statistik. Ini menunjukkan bukti yang kuat bahwa ada hubungan spasial negatif yang signifikan antara kedua variabel yang dianalisis. Secara keseluruhan, hasil ini menunjukkan bahwa kriminalitas dan nilai rumah memiliki hubungan spasial negatif yang signifikan, dengan korelasi yang lebih kuat di wilayah yang lebih dekat satu sama lain. berikut merupakan Plot yang dihasilkan:
ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) +
geom_point(color = "blue") +
geom_smooth(method = "lm", se = FALSE, color = "red", linetype = "dashed") +
labs(title = "Lee's Scatter plot (CRIM vs MEDV)",
x = "CRIM (Kriminalitas)",
y = "MEDV (Nilai Median Rumah)") +
theme_minimal()
Dari Lee’s Scatter Plot yang ditampilkan, menggambarkan hubungan antara CRIM (kriminalitas) dan MEDV (nilai median rumah). Garis tren merah putus-putus (garis regresi linier) menunjukkan hubungan linier negatif antara kedua variabel. Garis ini memperlihatkan bahwa kejahatan yang lebih tinggi di suatu area cenderung berhubungan dengan nilai rumah yang lebih rendah, yang konsisten dengan hasil Lee’s L statistic yang menunjukkan asosiasi spasial negatif.
Secara keseluruhan, plot ini memperkuat temuan bahwa ada asosiasi spasial negatif yang kuat antara kriminalitas dan nilai rumah, di mana area dengan tingkat kriminalitas lebih tinggi cenderung memiliki nilai rumah yang lebih rendah.
# Menggunakan Simulasi Monte Carlo
lee.mc(x, y, nsim=99, lw, zero.policy=TRUE, alternative="two.sided")
Monte-Carlo simulation of Lee's L
data: x , y
weights: lw
number of simulations + 1: 100
statistic = -0.3263, observed rank = 1, p-value = 0.02
alternative hypothesis: two.sidedHasil Monte-Carlo simulation ini memperkuat temuan bahwa ada korelasi spasial negatif yang signifikan antara kedua variabel yang dianalisis, dan ini menunjukkan bahwa nilai rendah pada satu variabelcenderung terkait dengan nilai rendah pada variabel lain di wilayah yang berdekatan.
Local Lee’s L
X <- boston.c$CMEDV
Y <- boston.c$CRIM
simula_lee <- function(x, y, listw, nsim = nsim, zero.policy = NULL, na.action = na.fail) {
if (deparse(substitute(na.action)) == "na.pass")
stop ("na.pass not permitted")
na.act <- attr(na.action(cbind(x, y)), "na.action")
x[na.act] <- NA
y[na.act] <- NA
x <- na.action(x)
y <- na.action(y)
if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
}
n <- length(listw$neighbours)
if ((n != length(x)) | (n != length(y)))
stop ("objects of different length")
gamres <- suppressWarnings(nsim > gamma(n + 1))
if (gamres)
stop ("nsim too large for this number of observations")
if (nsim < 1)
stop ("nsim too small")
xy <- data.frame(x, y)
S2 <- sum((unlist(lapply(listw$weights, sum)))^2)
lee_boot <- function(var, i, ...) {
return(spdep::lee(x = var[i, 1], y = var[i, 2], ...)$localL)
}
res <- boot::boot(xy, statistic = lee_boot, R = nsim, sim = "permutation",
listw = listw, n = n, S2 = S2, zero.policy = zero.policy)
}
lw <- nb2listw(boston.soi, style = "B", zero.policy = T)
W <- as(lw, "symmetricMatrix")
W <- as.matrix(W / Matrix::rowSums(W))
W[which(is.na(W))] <- 0
# Local Lee's L values
m <- lee(x = X,
y = Y,
listw = lw,
n = length(X),
zero.policy = TRUE,
NAOK = TRUE)
# Local Lee's L simulations
local_sims <- simula_lee(x = X,
y = Y,
listw = lw,
nsim = 10000,
zero.policy = TRUE,
na.action = na.omit)
m_i <- m[[2]] # local values
# Identify the significant values
alpha <- 0.05 # for a 95% confidence interval
probs <- c(alpha/2, 1-alpha/2)
intervals <- t(apply(t(local_sims[[2]]), 1, function(x) quantile(x, probs = probs)))
sig <- ( m_i < intervals[ , 1] ) | ( m_i > intervals[ , 2] )# Prepare for plotting
map_sf <- sf::st_as_sf(data.frame(boston.utm), coords = c("x", "y"))
map_sf$sig <- sig
# Identify the Lee's L clusters
Xp <- scale(X)[ , 1]
Yp <- scale(Y)[ , 1]
patterns <- as.character(interaction(Xp > 0, W %*% Yp > 0))
patterns <- patterns %>%
stringr::str_replace_all("TRUE","High") %>%
stringr::str_replace_all("FALSE","Low")
patterns[map_sf$sig == 0] <- "Not significant"
map_sf$patterns <- patterns
# Rename Lee's L clusters
map_sf$patterns2 <- factor(map_sf$patterns,
levels = c("High.High", "High.Low", "Low.High", "Low.Low", "Not significant"),
labels = c("High value - High crime", "High value - Low crime", "Low value - High crime","Low value - Low crime", "Not significant"))
# Plot
ggplot() +
geom_sf(data = map_sf, ggplot2::aes(color = patterns2)) +
scale_color_manual(values = c("red", "pink", "light blue", "dark blue", "grey80")) +
guides(color = ggplot2::guide_legend(title = "Lee's L clusters")) +
theme_minimal()
Referensi:
Melo, C., Melo, O., & Melo, S. (2024). Statistical Tools for Bivariate Spatial Dependence Analysis (bispdep) (Version 1.0-1). CRAN. https://github.com/carlosm77/bispdep.
See Lee (2004) for details on how the asymptotic expectation and variance of Lee’s L is computed. In particular, check Lee (2004), table 1, page 1690. https://en.wikipedia.org/wiki/Lee%27s_L