Spatial Autocorrelation

Joint Count Statistics

Joint count merupakan metode paling dasar dalam menentukan autokorelasi spasial antar area. Penjelasan lebih lengkap mengenai metode ini dapat dilihat pada Lizazaro (2016). Ilustrasi berikut ini diadaptasi dari Lizazaro (2016). Untuk melakukan analisis joint count, terlebih dulu kita load beberapa package pada R.

library(raster)
library(sp)
library(spdep)

Sebagai ilustrasi, kita akan membuat data contoh dengan menggunakan syntax berikut.

pri <- rep(1,12)
seg <- rep(0,4)
ter <- rep(1,2)
cua <- rep(0,4)
qui <- rep(1,2)
sex <-  rep(0,12)

A <- matrix(c(pri, seg, ter, cua ,qui, sex), nrow=6, byrow=FALSE)

A

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    0    0    0    0
## [2,]    1    1    0    0    0    0
## [3,]    1    1    0    0    0    0
## [4,]    1    1    0    0    0    0
## [5,]    1    1    1    1    0    0
## [6,]    1    1    1    1    0    0

Selanjutnya matriks A dikonversi menjadi raster menggunakan fungsi raster().

rA <- raster(A)

rA

## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

Berikut ini adalah plot dari area A.

plot(rA)
text(coordinates(rA), labels=rA[], cex=1.5)

Kriteria ketetanggaan selanjutnya digunakan untuk mengkuantifikasi pola kedekatan data tersebut.

pA <- rasterToPolygons(rA, dissolve=FALSE)

pA

## class       : SpatialPolygonsDataFrame 
## features    : 36 
## extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 1
## names       : layer 
## min values  :     0 
## max values  :     1

Seandainya kita akan menggunakan kriteria queen contiguity, maka dapat dilakukan dengan syntax berikut.

spA <- SpatialPolygons(pA@polygons)

nb1 <- poly2nb(spA, queen = T)

nb1

## Neighbour list object:
## Number of regions: 36 
## Number of nonzero links: 220 
## Percentage nonzero weights: 16.97531 
## Average number of links: 6.111111

Selanjutnya kita dapat memvisualisasikan link dari ketetanggannya.

par(mai=c(0,0,0,0))
plot(spA, col='gray', border='blue')
xy <- coordinates(spA)
plot(nb1, xy, col='red', lwd=2, add=TRUE)

Seandainya kriteria yang digunakan adalah rook contiguity.

nb2 <- poly2nb(spA, queen = F)

nb2

## Neighbour list object:
## Number of regions: 36 
## Number of nonzero links: 120 
## Percentage nonzero weights: 9.259259 
## Average number of links: 3.333333

par(mai=c(0,0,0,0))
plot(spA, col='gray', border='blue')
xy <- coordinates(spA)
plot(nb2, xy, col='green', lwd=2, add=TRUE)

Agar dapat mengidentifikasi pola autokorelasi pada data ini, maka matriks bobot perlu disusun dengan tipe biner.

wl1 <- nb2listw(nb1, style='B')

wl2 <- nb2listw(nb2, style='B')

jc_test1 <- joincount.test(as.factor(pA$layer), wl1)

jc_test1

## 
##  Join count test under nonfree sampling
## 
## data:  as.factor(pA$layer) 
## weights: wl1 
## 
## Std. deviate for 0 = 5.1529, p-value = 1.282e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##              53.00000              33.17460              14.80263 
## 
## 
##  Join count test under nonfree sampling
## 
## data:  as.factor(pA$layer) 
## weights: wl1 
## 
## Std. deviate for 1 = 4.7634, p-value = 9.52e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##              37.00000              20.95238              11.34999

jc_test2 <- joincount.test(as.factor(pA$layer), wl2)

jc_test2

## 
##  Join count test under nonfree sampling
## 
## data:  as.factor(pA$layer) 
## weights: wl2 
## 
## Std. deviate for 0 = 5.4677, p-value = 2.28e-08
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             30.000000             18.095238              4.740611 
## 
## 
##  Join count test under nonfree sampling
## 
## data:  as.factor(pA$layer) 
## weights: wl2 
## 
## Std. deviate for 1 = 5.1203, p-value = 1.525e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             22.000000             11.428571              4.262554

Berdasarkan kedua output di atas, dengan p-value yang sangat kecil, artinya kita dapat menolak hipotesis nol yang menyatakan bahwa tidak terdapat autokorelasi. Sesuai dengan output di atas, alternative hypothesis: greater, artinya kita dapat menyimpulkan bahwa terdapat cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf nyata 5%.

Pengujian hipotesis dapat pula dilakukan dengan melibatkan algoritma monte carlo seperti di bawah ini.

set.seed(123)
jc_test3 <- joincount.mc(as.factor(pA$layer), wl1, nsim=99)

jc_test3

## 
##  Monte-Carlo simulation of join-count statistic
## 
## data:  as.factor(pA$layer) 
## weights: wl1 
## number of simulations + 1: 100 
## 
## Join-count statistic for 0 = 53, rank of observed statistic = 100,
## p-value = 0.01
## alternative hypothesis: greater
## sample estimates:
##     mean of simulation variance of simulation 
##               33.01010               15.15296 
## 
## 
##  Monte-Carlo simulation of join-count statistic
## 
## data:  as.factor(pA$layer) 
## weights: wl1 
## number of simulations + 1: 100 
## 
## Join-count statistic for 1 = 37, rank of observed statistic = 100,
## p-value = 0.01
## alternative hypothesis: greater
## sample estimates:
##     mean of simulation variance of simulation 
##               20.37374               12.37930

Global Autocorrelation

Moran’s I

Sebagai ilustrasi, kita masih akan menggunakan data yang sama. Seandainya kita ingin menguji autokorelasi menggunakan pendekatan indeks moran, maka kita dapat menggunakan fungsi moran.test().

I1 <- moran.test(pA$layer,wl1)

I1

## 
##  Moran I test under randomisation
## 
## data:  pA$layer  
## weights: wl1    
## 
## Moran I statistic standard deviate = 7.4654, p-value = 4.151e-14
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.624090909      -0.028571429       0.007643049

Sama halnya seperti yang telah diperlihatkan pada uji dengan metode joint count, uji moran juga dapat dilakukan dengan melibatkan simulasi monte carlo.

set.seed(123)
MC<- moran.mc(pA$layer, wl1, nsim=599)

# View results (including p-value)
MC

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  pA$layer 
## weights: wl1  
## number of simulations + 1: 600 
## 
## statistic = 0.62409, observed rank = 600, p-value = 0.001667
## alternative hypothesis: greater

Global Geary’s C

Geary’s C merupakan alternatif dari indeks Moran, yang memiliki nilai antara 0 s.d 2. Nilai 0 menunjukkan autokorelasi positif, 1 menunjukkan tidak ada autokorelasi, dan 2 menunjukkan autokorelasi negatif.

C1 <- geary.test(pA$layer,wl1)

C1

## 
##  Geary C test under randomisation
## 
## data:  pA$layer 
## weights: wl1 
## 
## Geary C statistic standard deviate = 7.4869, p-value = 3.526e-14
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.357954545       1.000000000       0.007354046

GS1 <- geary.mc(pA$layer, wl1, nsim=599)

GS1

## 
##  Monte-Carlo simulation of Geary C
## 
## data:  pA$layer 
## weights: wl1 
## number of simulations + 1: 600 
## 
## statistic = 0.35795, observed rank = 1, p-value = 0.001667
## alternative hypothesis: greater

Local Autocorrelation

Local Moran’s I

Pendekatan ini termasuk ke dalam Local Indicators for Spatial Association (LISA), yang mengindentifikasi autokorelasi pada tingkat lokal.

oid <- order(pA$layer)
resI <- localmoran(pA$layer, wl1)
head(resI)

##     Ii        E.Ii   Var.Ii      Z.Ii  Pr(z > 0)
## 1 3.75 -0.08571429 2.817443 2.2851709 0.01115140
## 2 1.75 -0.14285714 4.402701 0.9021074 0.18349991
## 3 0.40 -0.14285714 4.402701 0.2587176 0.39792657
## 4 4.00 -0.14285714 4.402701 1.9744237 0.02416679
## 5 4.00 -0.14285714 4.402701 1.9744237 0.02416679
## 6 2.40 -0.08571429 2.817443 1.4808929 0.06931756

pA$z.li <- resI[,4]
pA$pvalue <- resI[,5]
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(pA, zcol="z.li", col.regions=lm.palette(20), main="Local Moran")

moran.plot(pA$layer,wl1)

Getis-Ord Gi

Menurut Mendez (2020), pendekatan Getis-ord Gi dapat membantu mengidentifikasi pola penggerombolan berdasarkan ukuran autokorelasi pada level lokal.

local_g <- localG(pA$layer, wl1)

local_g

##  [1]  2.0615528  0.8246211 -0.2730593 -2.1844747 -2.1844747 -1.6383560
##  [7]  2.7487371  1.2598378 -0.5233637 -2.9126330 -2.9126330 -2.1844747
## [13]  2.7487371  1.2598378 -0.5233637 -2.9126330 -2.9126330 -2.1844747
## [19]  2.7487371  2.0615528  1.0694824 -1.3197868 -2.1162099 -2.1844747
## [25]  2.7487371  2.8632678  2.0615528 -0.3435921 -1.3197868 -2.1844747
## [31]  2.0615528  2.7487371  2.7487371  0.8246211 -0.2730593 -1.6383560
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = pA$layer, listw = wl1)
## attr(,"class")
## [1] "localG"

Output di atas menghasilkan z-score, yang biasanya disajikan secara visual untuk mengidentifikasi cluster maupun hotspot.

pA$localg <- as.numeric(local_g)
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(pA, zcol="localg", col.regions=lm.palette(20), main="Local Gi")

Exercise

Sebagai latihan, Anda dipersilahkan menggunakan data yang tersedia pada: https://github.com/raoy/SpatialReg . Terdapat dua data yang harus Anda download, yaitu:

Jabar Data (gabung).xlsx
petaJabar2.zip

Data pertama (dengan format Excel) menyimpan data kependudukan yang diperoleh dari BPS. Sedangkan data kedua merupakan data shapefile berisi peta Provinsi Jawa Barat. Silahkan manfaatkan kedua data tersebut untuk mengeksplorasi pola depedensi spasial untuk peubah kemiskinan antar kota/kabupaten di Jawa Barat pada tahun 2015. Data tersebut terdapat pada kolom I dengan nama kolom p.miskin15 pada file Excel.

References

Lizarazo, I. (2016, October 31). Areal pattern analysis. https://rstudio-pubs-static.s3.amazonaws.com/223305_944ddc517306448f8fb0d60ca29dd94b.html

Mendez, C. (2020). Spatial autocorrelation analysis in R. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/spatial-autocorrelation

rahmaanisa@apps.ipb.ac.id ↩︎

Spatial Autocorrelation

Rahma Anisa¹

November 12, 2020

Joint Count Statistics

Global Autocorrelation

Moran’s I

Global Geary’s C

Local Autocorrelation

Local Moran’s I

Getis-Ord Gi

Exercise

References

Spatial Autocorrelation

Rahma Anisa1

November 12, 2020

Joint Count Statistics

Global Autocorrelation

Moran’s I

Global Geary’s C

Local Autocorrelation

Local Moran’s I

Getis-Ord Gi

Exercise

References

Rahma Anisa¹