Autokorelasi Spasial (1)

Review: Autokorelasi Temporal

set.seed(0)
d<-(1:20)+rnorm(20)

a <- d[-length(d)]
b <- d[-1]
plot(a, b, xlab='t', ylab='t-1')

acf(d)

Pemodelan Proses Spasial

Menurut Cressie dan Moores (2021), suatu pemodelan proses spasial dapat dibedakan menjadi tiga jenis, yaitu:

• point process (titik)

• lattice process (bersifat diskret)

• geostatistical process (bersifat kontinu)

Autokorelasi Spasial pada Lattice Process

Ilustrasi ini diambil dari Georgia Tech City Planning (2019), yang mengadopsi materi dari Gimond (2019).

Ilustrasi: Data Pendapatan di Maine

Sebelum menjalankan fungsi-fungsi pada bagian ini, kita akan mempersiapkan dulu beberapa package yang akan digunakan untuk analisis dan perhitungan autkorelasi spasial.

library(sp)
library(spdep)

## Loading required package: spData

## 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

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

library(tmap)

Data yang digunakan berikut ini diperoleh dari Gimond (2021) pada link ini.

load(url("https://github.com/raoy/data/raw/master/moransI.RData"))

class(s1)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

str(s1, max.level = 2)

## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  16 obs. of  5 variables:
##   ..@ polygons   :List of 16
##   ..@ plotOrder  : int [1:16] 1 3 2 4 5 7 6 11 16 9 ...
##   ..@ bbox       : num [1:2, 1:2] 336338 4772272 660529 5255569
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   ..$ comment: chr "TRUE"

s1@data

##            NAME Income   NoSchool  NoSchoolSE IncomeSE
## 0     Aroostook  21024 0.01338720 0.001406960  250.909
## 1      Somerset  21025 0.00521153 0.001150020  390.909
## 2   Piscataquis  21292 0.00633830 0.002128960  724.242
## 3     Penobscot  23307 0.00684534 0.001025450  242.424
## 4    Washington  20015 0.00478188 0.000966036  327.273
## 5      Franklin  21744 0.00508507 0.001641740  530.909
## 6        Oxford  21885 0.00700822 0.001318160  536.970
## 7         Waldo  23020 0.00498141 0.000918837  450.909
## 8      Kennebec  25652 0.00570358 0.000917087  360.000
## 9  Androscoggin  24268 0.00830953 0.001178660  460.606
## 10      Hancock  28071 0.00238996 0.000784584  585.455
## 11         Knox  27141 0.00652269 0.001863920  684.849
## 12      Lincoln  27839 0.00278315 0.001030800  571.515
## 13   Cumberland  32549 0.00494917 0.000683236  346.061
## 14    Sagadahoc  28122 0.00285524 0.000900782  544.849
## 15         York  28496 0.00529228 0.000737195  332.121

tmap_mode("view")
tm_basemap("Stamen.TonerLite") + 
  tm_shape(s1) +
  tm_polygons(style = "quantile", col = "Income", alpha = 0.5) +
  tm_legend(outside = TRUE, text.size = .8)

Menentukan Ketetanggaan

Ketetanggaan dapat mengacu pada beberapa kriteria, di antaranya adalah:

• contiguity

• jarak

• n-tetangga terdekat

• general weights

Terdapat beberapa kriteria contiguity seperti terlihat pada ilustrasi berikut.

Penyusunan matriks pembobot spasial dengan kriteria queen dan rook contiguity dapat dilakukan menggunakan fungsi poly2nb()

nb <- poly2nb(s1, queen=TRUE)

class(nb) 

## [1] "nb"

str(nb)

## List of 16
##  $ : int [1:4] 2 3 4 5
##  $ : int [1:6] 1 3 4 6 8 9
##  $ : int [1:3] 1 2 4
##  $ : int [1:6] 1 2 3 5 8 11
##  $ : int [1:3] 1 4 11
##  $ : int [1:4] 2 7 9 10
##  $ : int [1:4] 6 10 14 16
##  $ : int [1:6] 2 4 9 11 12 13
##  $ : int [1:6] 2 6 8 10 13 15
##  $ : int [1:5] 6 7 9 14 15
##  $ : int [1:3] 4 5 8
##  $ : int [1:2] 8 13
##  $ : int [1:4] 8 9 12 15
##  $ : int [1:4] 7 10 15 16
##  $ : int [1:4] 9 10 13 14
##  $ : int [1:2] 7 14
##  - attr(*, "class")= chr "nb"
##  - attr(*, "region.id")= chr [1:16] "0" "1" "2" "3" ...
##  - attr(*, "call")= language poly2nb(pl = s1, queen = TRUE)
##  - attr(*, "type")= chr "queen"
##  - attr(*, "sym")= logi TRUE

summary(nb)

## Neighbour list object:
## Number of regions: 16 
## Number of nonzero links: 66 
## Percentage nonzero weights: 25.78125 
## Average number of links: 4.125 
## Link number distribution:
## 
## 2 3 4 5 6 
## 2 3 6 1 4 
## 2 least connected regions:
## 11 15 with 2 links
## 4 most connected regions:
## 1 3 7 8 with 6 links

nb[1]

## [[1]]
## [1] 2 3 4 5

class(nb[1])

## [1] "list"

nb[[1]]

## [1] 2 3 4 5

class(nb[[1]])

## [1] "integer"

Selanjutnya kita dapat pula memeriksa link antar lokasi yang bertetangga dengan cara berikut ini.

xy <- sp::coordinates(s1)
plot(s1, border='grey', lwd=1)
plot(nb, xy, col='red', lwd=1, add=TRUE)

Kita dapat memperoleh nama lokasi yang bertetangga dengan perintah berikut.

s1$NAME

##  [1] Aroostook    Somerset     Piscataquis  Penobscot    Washington  
##  [6] Franklin     Oxford       Waldo        Kennebec     Androscoggin
## [11] Hancock      Knox         Lincoln      Cumberland   Sagadahoc   
## [16] York        
## 16 Levels: Androscoggin Aroostook Cumberland Franklin Hancock Kennebec ... York

s1$NAME[1]

## [1] Aroostook
## 16 Levels: Androscoggin Aroostook Cumberland Franklin Hancock Kennebec ... York

class(s1$NAME[1])

## [1] "factor"

s1$NAME[nb[[1]]]

## [1] Somerset    Piscataquis Penobscot   Washington 
## 16 Levels: Androscoggin Aroostook Cumberland Franklin Hancock Kennebec ... York

Selanjutnya, kita perlu mempersiapkan object listw agar dapat menggunakan fungsi moran() untuk menghitung autokorelasi dengan pendekatan indeks Moran.

lw <- nb2listw(nb, style="W", zero.policy=TRUE)

class(lw)

## [1] "listw" "nb"

lw$weights

## [[1]]
## [1] 0.25 0.25 0.25 0.25
## 
## [[2]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## 
## [[3]]
## [1] 0.3333333 0.3333333 0.3333333
## 
## [[4]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## 
## [[5]]
## [1] 0.3333333 0.3333333 0.3333333
## 
## [[6]]
## [1] 0.25 0.25 0.25 0.25
## 
## [[7]]
## [1] 0.25 0.25 0.25 0.25
## 
## [[8]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## 
## [[9]]
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## 
## [[10]]
## [1] 0.2 0.2 0.2 0.2 0.2
## 
## [[11]]
## [1] 0.3333333 0.3333333 0.3333333
## 
## [[12]]
## [1] 0.5 0.5
## 
## [[13]]
## [1] 0.25 0.25 0.25 0.25
## 
## [[14]]
## [1] 0.25 0.25 0.25 0.25
## 
## [[15]]
## [1] 0.25 0.25 0.25 0.25
## 
## [[16]]
## [1] 0.5 0.5
## 
## attr(,"mode")
## [1] "binary"
## attr(,"W")
## [1] TRUE
## attr(,"comp")
## attr(,"comp")$d
##  [1] 4 6 3 6 3 4 4 6 6 5 3 2 4 4 4 2

lw$weights[[1]]

## [1] 0.25 0.25 0.25 0.25

Matriks pembobot spasial yang digunakan dapat berupa matriks biner (berisi nilai nol atau 1), row standardized, dsb.

• Contoh binary weight matrix

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

• Contoh row-standardized weight matrix

##      [,1] [,2] [,3] [,4]
## [1,]  0.0 1.00 0.00 0.00
## [2,]  0.0 0.00 0.50 0.50
## [3,]  0.5 0.50 0.00 0.00
## [4,]  0.0 0.33 0.33 0.33

Berikut adalah style yang tersedia pada fungsi nb2listw().

Code	Description
B	basic binary coding
W	row standardised (sums over all links to n)
C	globally standardised (sums over all links to n)
U	equal to C divided by the number of neighbours (sums over all links to unity)
S	the variance-stabilizing coding scheme

moran(s1$Income, lw, n=length(lw$neighbours), S0=Szero(lw))

## $I
## [1] 0.2828111
## 
## $K
## [1] 2.249902

Output di atas menunjukkan bahwa pendapatan di Maine memiliki indeks Moran sebesar 0.2828111 dan kurtosis sebesar 2.249902.

Pengujian hipotesis dapat menggunakan fungsi moran.test() dan moran.mc(). Secara default, hipotesis yang digunakan adalah \(H_A\): \(I > 0\). Kita dapat menambahkan argumen alternative untuk mengganti hipotesis tersebut.

moran.test(s1$Income,lw)

## 
##  Moran I test under randomisation
## 
## data:  s1$Income  
## weights: lw    
## 
## Moran I statistic standard deviate = 2.2472, p-value = 0.01231
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.28281108       -0.06666667        0.02418480

Output di atas menunjukkan bahwa p-value = 0.01231, artinya kita dapat menolak hipotesis nol pada taraf 0.05 dan menyimpulkan bahwa terdapat autokorelasi spasial pada data pendapatan tersebut.

Pengujian dapat pula dilakukan secara permutasi menggunakan pendekatan Monte Carlo seperti berikut ini.

moran.mc(s1$Income,lw, nsim=99)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  s1$Income 
## weights: lw  
## number of simulations + 1: 100 
## 
## statistic = 0.28281, observed rank = 100, p-value = 0.01
## alternative hypothesis: greater

• Moran Scatterplot

Moran plot dapat dibuat menggunakan fungsi moran.plot() seperti berikut ini.

moran.plot(s1$Income, lw, labels=s1$NAME)

dengan \(x_i\) merupakan nilai dari peubah \(x\) yang diamati di lokasi \(i\). Sementara \(\overline(x)\) merupakan nilai rataan dari peubah \(x\) pada semua lokasi dan \(s_x\) adalah simpangan baku dari peubah \(x\). Morans Scatterplot disajikan berbasis pada data z-score suatu lokasi pada sumbu (x), dan nilai z-score rata-rata tetangganya pada sumbu (y). Secara visual Morans Scatterplot terbagi 4 kuadran.

• Local Moran

Autokorelasi lokal dapat dihitung menggunakan beberapa pendekatan, salah satunya adalah dengan Indeks Moran, yang dikenal pula sebagai Local Indicators of Spatial Association (LISA) fungsi localmoran(). Fungsi tersebut akan menghasilkan nilai-nilai berikut:

• Ii: statistik local Moran

• E.Ii: nilai harapan dari statistik local Moran

• Var.Ii: ragam dari statistik local Moran

• Z.Ii: z-score yang dapat dihitung dengan formula \(Z_i=\frac{I_i-E_i}{\sqrt{Var_i}}\)

• Pr(z!=E(Ii)): p-value

local <- localmoran(s1$Income, listw = lw)
head(local)

##          Ii        E.Ii    Var.Ii      Z.Ii Pr(z != E(Ii))
## 0 0.9861946 -0.07341799 0.2138016 2.2916147     0.02192789
## 1 0.6091324 -0.07337822 0.1165609 1.9990930     0.04559830
## 2 0.8106915 -0.06314511 0.2704357 1.6803442     0.09289036
## 3 0.2626943 -0.01068952 0.0181290 2.0304195     0.04231392
## 4 0.2209313 -0.11903606 0.4793896 0.4910129     0.62341730
## 5 0.3621512 -0.04757243 0.1424006 1.0857630     0.27758385

Pada umumnya local Moran disajikan secara visual seperti berikut ini (Mendez, 2020).

moran.map<-cbind(s1,local)
tm_shape(moran.map) +
  tm_fill(col = "Ii",
          style = "quantile",
          title = "local moran statistic") 

• Plot LISA Clusters

quadrant <- vector(mode="numeric",length=nrow(local))

# centers the variable of interest around its mean
m.qualification <- s1$Income - mean(s1$Income)     

# centers the local Moran's around the mean
m.local <- local[,1] - mean(local[,1])    

# significance threshold
signif <- 0.05

# builds a data quadrant
quadrant[m.qualification >0 & m.local>0] <- 4  
quadrant[m.qualification <0 & m.local<0] <- 1      
quadrant[m.qualification <0 & m.local>0] <- 2
quadrant[m.qualification >0 & m.local<0] <- 3
quadrant[local[,5]>signif] <- 0   

# plot in r
brks <- c(0,1,2,3,4)
colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")
plot(s1,border="lightgray",col=colors[findInterval(quadrant,brks,all.inside=FALSE)])

## Warning in wkt(obj): CRS object has no comment

box()
legend("bottomleft", legend = c("insignificant","low-low","low-high","high-low","high-high"),
       fill=colors,bty="n")

Terdapat beberapa pendekatan lain yang dapat digunakan untuk menganalisa autokorelasi spasial selain Indeks Moran, di antaranya adalah Getis-Ord dan Geary’s C. Perhitungan Getis-Ord pada R software dapat menggunakan fungsi globalG.test() dan localG(). Sedangkan perhitungan untuk Geary’s C dapat menggunakan fungsi geary.test().

Autokorelasi Spasial pada Geostatistical Process

Spatial Autocorrelogram

Ilustrasi berikut ini mengacu pada Petrkeil (2013) yang dapat diakses pada link ini.

# packages used for the data generation
library(raster)
library(colorRamps) # for some crispy colors
library(vegan) # will be used for PCNM
# empty matrix and spatial coordinates of its cells
side=30
my.mat <- matrix(NA, nrow=side, ncol=side)
x.coord <- rep(1:side, each=side)
y.coord <- rep(1:side, times=side)
xy <- data.frame(x.coord, y.coord)
# all paiwise euclidean distances between the cells
xy.dist <- dist(xy)
# PCNM axes of the dist. matrix (from 'vegan' package)
pcnm.axes <- pcnm(xy.dist)$vectors
# using 8th PCNM axis as my atificial z variable
z.value <- pcnm.axes[,8]*200 + rnorm(side*side, 0, 1)
# plotting the artificial spatial data
my.mat[] <- z.value
r <- raster(my.mat)
plot(r, axes=F, col=matlab.like(20))

set.seed(1)
library(ncf)

## 
## Attaching package: 'ncf'

## The following object is masked from 'package:vegan':
## 
##     mantel.correlog

ncf.cor <- correlog(x.coord, y.coord, z.value,
                    increment=2, resamp=100)

## 10  of  100 
20  of  100 
30  of  100 
40  of  100 
50  of  100 
60  of  100 
70  of  100 
80  of  100 
90  of  100 
100  of  100 

plot(ncf.cor)

# 'nb' - neighbourhood of each cell
r.nb <- dnearneigh(as.matrix(xy), d1=0.5, d2=1.5)
# 'nb' - an alternative way to specify the neighbourhood
# r.nb <- cell2nb(nrow=side, ncol=side, type="queen")
sp.cor <- sp.correlogram(r.nb, z.value, order=15,
                      method="I",randomisation=FALSE)
plot(sp.cor)

Variogram

Ilustrasi berikut ini mengacu pada Guélat (2013) yang dapat diakses pada link ini.

library(gstat)
xy$z.value<-z.value
coordinates(xy)<-c("x.coord", "y.coord")
variocloud <- variogram(z.value ~ 1, data=xy, cloud = TRUE)
plot(variocloud)

vario <- variogram(z.value ~ 1, xy, cutoff = 30)
plot(vario, pch = 16, cex = 1.5)

Referensi

Cressie, N., & Moores, M. T. (2021). Spatial statistics. arXiv preprint arXiv:2105.07216. Retrieved from: https://arxiv.org/pdf/2105.07216.pdf

Georgia Tech City Planning. (2019, April 3). Spatial autocorrelation. RPubs. Retrieved from: https://rpubs.com/spring19cp6521/Week12_Wednesday1

Gimond, M. (2021). Intro to GIS and spatial analysis. Main repos | mgimond.github.io. Retrieved from: https://mgimond.github.io/Spatial/index.html

Guélat, J. (2013). Spatial autocorrelation (introduction). Retrieved from: https://rstudio-pubs-static.s3.amazonaws.com/9688_a49c681fab974bbca889e3eae9fbb837.html

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

Petrkeil. (2013, May 21). Spatial correlograms in R: A mini overview. R-bloggers. https://www.r-bloggers.com/2013/05/spatial-correlograms-in-r-a-mini-overview/