Spatial Autocorrelation

Excercise 9

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

  1. Jabar Data (gabung).xlsx

  2. 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.

Import Data

library(raster)
library(sp)
library(spdep)
library(readxl)
library(rgdal)

#import data kependudukan
Jabardata<-read_excel("D:/MATERI KULIAH S2 IPB/SEMESTER 2/SPASIAL/Jabar Data (gabung).xlsx", sheet = "data")
head(Jabardata)
## # A tibble: 6 x 32
##   PROVNO KABKOTNO KODE2010 PROVINSI   KABKOT IDSP2~1  Long   Lat p.mis~2 p.mis~3
##    <dbl>    <dbl>    <dbl> <chr>      <chr>    <dbl> <dbl> <dbl>   <dbl>   <dbl>
## 1     32        1     3201 JAWA BARAT BOGOR     3201  107. -6.56    8.96    8.83
## 2     32        2     3202 JAWA BARAT SUKAB~    3202  107. -7.07    8.96    8.13
## 3     32        3     3203 JAWA BARAT CIANJ~    3203  107. -7.13   12.2    11.6 
## 4     32        4     3204 JAWA BARAT BANDU~    3204  108. -7.10    8.00    7.61
## 5     32        5     3205 JAWA BARAT GARUT     3205  108. -7.36   12.8    11.6 
## 6     32        6     3206 JAWA BARAT TASIK~    3206  108. -7.50   12.0    11.2 
## # ... with 22 more variables: j.miskin15 <dbl>, j.miskin16 <dbl>,
## #   AHH2015 <dbl>, AHH2016 <dbl>, EYS2015 <dbl>, EYS2016 <dbl>, MYS2015 <dbl>,
## #   MYS2016 <dbl>, EXP2015 <dbl>, EXP2016 <dbl>, APM.SD15 <dbl>,
## #   APM.SMP15 <dbl>, APM.SMA15 <dbl>, APM.PT15 <dbl>, APK.SD15 <dbl>,
## #   APK.SMP15 <dbl>, APK.SMA15 <dbl>, APK.PT15 <dbl>, APS.USIA15 <dbl>,
## #   APS.USIA2 <dbl>, APS.USIA3 <dbl>, APS.USIA4 <dbl>, and abbreviated variable
## #   names 1: IDSP2010, 2: p.miskin15, 3: p.miskin16
#import shapefile
petajabar<-readOGR(dsn="D:/MATERI KULIAH S2 IPB/SEMESTER 2/SPASIAL/petaJabar2", layer="Jabar2") #dsn diisi nama folder #layer diisi nama file dalam folder
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\MATERI KULIAH S2 IPB\SEMESTER 2\SPASIAL\petaJabar2", layer: "Jabar2"
## with 26 features
## It has 7 fields
#plot peta jabar 
plot(petajabar) 
text(petajabar,'KABKOT',cex=0.5) #menambahkan nama wilayah pada peta

Eksplorasi Data

library(raster)
colfunc<-colorRampPalette(c("yellow", "orange","red")) 
petajabar$miskin<-Jabardata$p.miskin15
spplot(petajabar, "miskin", col.regions=colfunc(16),
       main="Peta Persentase Penduduk Miskin di Jawa Barat Tahun 2015")

Matriks Bobot

Distance Matrix

# Distance Matrix
longlat<-cbind(Jabardata$Long ,Jabardata$Lat)
plot(longlat)

gdist<-pointDistance(longlat,lonlat=TRUE) #Distance for longitude/latitude coordinates
m.gdist<-as.matrix(gdist)


djarak<-dist(longlat) #Euclide
m.djarak<-as.matrix(djarak)

K-Nearest Neighbour Weight

#k=3
koord <- coordinates(petajabar)
W1<-knn2nb(knearneigh(longlat,k=3,longlat=TRUE)) #matriks bobot dengan knn k=3 
W1
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 78 
## Percentage nonzero weights: 11.53846 
## Average number of links: 3 
## Non-symmetric neighbours list

Standarisasi bobot:

W1<- nb2listw(W1,style='W') #W is row standardised (sums over all links to n). Standardisasi Baris
W1
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 78 
## Percentage nonzero weights: 11.53846 
## Average number of links: 3 
## Non-symmetric neighbours list
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 26 676 26 15.55556 107.7778
plot(petajabar, col='white', border='navy', main ="KNN (k=3)")
plot(W1, longlat, col='blue', lwd=2, add=TRUE)
text(petajabar,'KABKOT',cex=0.5, col='red')

Radial Distance Weight

#d=60
W2<-dnearneigh(koord,0,60,longlat=TRUE) #dnearneigh(x, d1, d2, row.names = NULL, longlat = NULL, bounds=c("GE", "LE"), use_kd_tree=TRUE, symtest=FALSE)

W2
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 142 
## Percentage nonzero weights: 21.00592 
## Average number of links: 5.461538

Standarisasi Bobot:

W2 <- nb2listw(W2,style='W') #W is row standardised (sums over all links to n). Standardisasi Baris
W2
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 142 
## Percentage nonzero weights: 21.00592 
## Average number of links: 5.461538 
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 26 676 26 9.871522 104.9136
plot(petajabar, col='white', border='navy', main ="Radial Distance Weight d=60")
plot(W2, longlat, col='blue', lwd=2, add=TRUE)
text(petajabar,'KABKOT',cex=0.5, col='red')

Power Distance Weight

Power distance weight dengan \(\alpha = 1\):

# Jarak euclidean
W3a<-1/(m.djarak^1) 

#dinormalisasi 
diag(W3a) <-0
rtot<-rowSums(W3a,na.rm=TRUE)

W3a = mat2listw(W3a,style='W')  
summary(W3a)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 650 
## Percentage nonzero weights: 96.15385 
## Average number of links: 25 
## Link number distribution:
## 
## 25 
## 26 
## 26 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 with 25 links
## 26 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 with 25 links
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 26 676 26 6.638168 105.5718
plot(petajabar, col='white', border='navy', main ="Power Distance Weight alpha=1")
plot(W3a, longlat, col='blue', lwd=2, add=TRUE)
text(petajabar,'KABKOT',cex=0.5, col='red')

Power distance weight dengan \(\alpha = 2\):

# Jarak Euclidean
W3b<-1/(m.djarak^2)

#dinormalisasi 
diag(W3b) <-0
rtot<-rowSums(W3b,na.rm=TRUE)

W3b<-W3b/rtot #row-normalized
rowSums(W3b,na.rm=TRUE) 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
W3b = mat2listw(W3b,style='W')#matriks bobot power distance dengan alpha=2 
summary(W3b)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 650 
## Percentage nonzero weights: 96.15385 
## Average number of links: 25 
## Link number distribution:
## 
## 25 
## 26 
## 26 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 with 25 links
## 26 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 with 25 links
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1      S2
## W 26 676 26 14.35111 108.236
plot(petajabar, col='white', border='navy', main ="Power Distance Weigth alpha=2, Euclidean")
plot(W3b, longlat, col='blue', lwd=2, add=TRUE)
text(petajabar,'KABKOT',cex=0.5, col='red')

Exponential Distance Weight

alpha=2
W4<-exp((-alpha)*m.djarak)

#dinormalisasi 
diag(W4) <-0
rtot<-rowSums(W4,na.rm=TRUE)
rtot
##        1        2        3        4        5        6        7        8 
## 5.385850 3.943443 5.449284 6.769776 5.421621 4.889553 4.646033 5.043533 
##        9       10       11       12       13       14       15       16 
## 5.099263 6.038772 6.611018 4.833578 5.986258 6.742714 5.272713 5.140898 
##       17       18       19       20       21       22       23       24 
## 7.169330 5.547803 5.403421 7.394625 5.076805 5.258575 5.189097 7.467575 
##       25       26 
## 5.487780 4.420439
W4<-W4/rtot #row-normalized
rowSums(W4,na.rm=TRUE)
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
W4 = mat2listw(W4,style='W')  
summary(W4)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 650 
## Percentage nonzero weights: 96.15385 
## Average number of links: 25 
## Link number distribution:
## 
## 25 
## 26 
## 26 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 with 25 links
## 26 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 with 25 links
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 26 676 26 3.497596 104.4985
plot(petajabar, col='white', border='navy', main ="Exponential Distance Weight Alpha=2")
plot(W4, longlat, col='blue', lwd=2, add=TRUE)

Spatial Contiguity Weight

plot(petajabar) 
text(petajabar,'KABKOT',cex=0.5) #menambahkan nama wilayah pada peta

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

Rook

#Rook
W5<-poly2nb(petajabar,queen=FALSE)
W5 #matriks bobot Rook
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 102 
## Percentage nonzero weights: 15.08876 
## Average number of links: 3.923077
#memetakan jabar dengan matriks bobot Rook
plot(petajabar,col='skyblue',border='white',main="Peta Persentase Kemiskinan Jabar \n Tahun 2015 dengan Rook")
xy<-coordinates(petajabar) 
plot(W5,xy,col='red',lwd=2,add=TRUE)

Queen

#Queen
W6<-poly2nb(petajabar,queen=TRUE)
W6 #matriks bobot Queen
## Neighbour list object:
## Number of regions: 26 
## Number of nonzero links: 102 
## Percentage nonzero weights: 15.08876 
## Average number of links: 3.923077
#memetakan jabar dengan matriks bobot Queen
par(mai=c(0,0,0,0))

plot(petajabar,col='skyblue',border='white',main="Peta Persentase Kemiskinan Jabar Tahun 2015 dengan Queen")
xy<-coordinates(petajabar) 
plot(W6,xy,col='orange',lwd=2,add=TRUE)

Global Autocorrelation

Moran’s I

K-Nearest Neighbor Weight dengan k=3:

MI1 <- moran.test(petajabar$miskin,W1)  
MI1  
## 
##  Moran I test under randomisation
## 
## data:  petajabar$miskin  
## weights: W1    
## 
## Moran I statistic standard deviate = 3.2398, p-value = 0.000598
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.42074017       -0.04000000        0.02022419

Radial Distance Weight

MI2 <- moran.test(petajabar$miskin,W2)  
MI2
## 
##  Moran I test under randomisation
## 
## data:  petajabar$miskin  
## weights: W2    
## 
## Moran I statistic standard deviate = 3.2719, p-value = 0.0005341
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.31496539       -0.04000000        0.01176964

Power Distance Weight

  1. Power distance weight dengan alpha = 1
MI3a <- moran.test(petajabar$miskin,W3a)  
MI3a
## 
##  Moran I test under randomisation
## 
## data:  petajabar$miskin  
## weights: W3a    
## 
## Moran I statistic standard deviate = 2.1932, p-value = 0.01415
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.141165530      -0.040000000       0.006823324
  1. Power distance weight dengan alpha = 2
MI3b <- moran.test(petajabar$miskin,W3b)  
MI3b
## 
##  Moran I test under randomisation
## 
## data:  petajabar$miskin  
## weights: W3b    
## 
## Moran I statistic standard deviate = 2.5972, p-value = 0.004699
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##         0.3120068        -0.0400000         0.0183689

Exponential distance weight

MI4 <- moran.test(petajabar$miskin,W4)  
MI4
## 
##  Moran I test under randomisation
## 
## data:  petajabar$miskin  
## weights: W4    
## 
## Moran I statistic standard deviate = 4.4975, p-value = 3.439e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.167150691      -0.040000000       0.002121481

Spatial Contiguity Weight

  1. Rook
MI5 <- moran.test(petajabar$miskin,W4)  
MI5
## 
##  Moran I test under randomisation
## 
## data:  petajabar$miskin  
## weights: W4    
## 
## Moran I statistic standard deviate = 4.4975, p-value = 3.439e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.167150691      -0.040000000       0.002121481
  1. Queen
MI6 <- moran.test(petajabar$miskin,W4)  
MI6
## 
##  Moran I test under randomisation
## 
## data:  petajabar$miskin  
## weights: W4    
## 
## Moran I statistic standard deviate = 4.4975, p-value = 3.439e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.167150691      -0.040000000       0.002121481

Simulasi Monte Carlo

Uji moran juga dapat dilakukan dengan melibatkan simulasi monte carlo.

set.seed(123)
MI7<- moran.mc(petajabar$miskin,W1, nsim=399)

# View results (including p-value)
MI7
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  petajabar$miskin 
## weights: W1  
## number of simulations + 1: 400 
## 
## statistic = 0.42074, observed rank = 397, p-value = 0.0075
## alternative hypothesis: greater

Hasil Global Moran’s I Autocorrelation

##                   Bobot.Spatial   Moran.I   p.value
## 1                    KNN Weight 0.4207402 5.980e-04
## 2        Radial Distance Weight 0.3149654 5.341e-04
## 3 Power Distance Weight alpha=1 0.1411655 1.415e-02
## 4 Power Distance Weight alpha=2 0.3120068 4.699e-03
## 5   Exponential Distance Weight 0.1671507 3.439e-06
## 6     Spatial Contiguity (Rook) 0.1671507 3.439e-06
## 7    Spatial Contiguity (Queen) 0.1671507 3.439e-06
## 8          Simulasi Monte Carlo 0.4207400 7.500e-03

Berdasarkan hasil di atas, diperoleh p-value < 0.05 berarti H0 ditolak, dengan alternative hypothesis: greater berarti pada taraf nyata 5%, terdapat autokorelasi spasial positif.

Geary

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.

KNN Weights

C1 <- geary.test(petajabar$miskin,W1) 
C1
## 
##  Geary C test under randomisation
## 
## data:  petajabar$miskin 
## weights: W1 
## 
## Geary C statistic standard deviate = 2.5241, p-value = 0.005799
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.63708305        1.00000000        0.02067246

Radial Distance Weight

C2 <- geary.test(petajabar$miskin,W2)  
C2
## 
##  Geary C test under randomisation
## 
## data:  petajabar$miskin 
## weights: W2 
## 
## Geary C statistic standard deviate = 3.1399, p-value = 0.000845
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.66535312        1.00000000        0.01135907

Power Distance Weight

  1. Power distance weight dengan alpha = 1
C3a <- geary.test(petajabar$miskin,W3a)  
C3a
## 
##  Geary C test under randomisation
## 
## data:  petajabar$miskin 
## weights: W3a 
## 
## Geary C statistic standard deviate = 2.0638, p-value = 0.01952
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.825735951       1.000000000       0.007129903
  1. Power distance weight dengan alpha = 2
C3b <- geary.test(petajabar$miskin,W3b)  
C3b
## 
##  Geary C test under randomisation
## 
## data:  petajabar$miskin 
## weights: W3b 
## 
## Geary C statistic standard deviate = 2.448, p-value = 0.007182
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##         0.6608234         1.0000000         0.0191966

Exponential distance weight

C4 <- geary.test(petajabar$miskin,W4)  
C4
## 
##  Geary C test under randomisation
## 
## data:  petajabar$miskin 
## weights: W4 
## 
## Geary C statistic standard deviate = 4.1302, p-value = 1.812e-05
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.805313190       1.000000000       0.002221941

Spatial Contiguity Weight

  1. Rook
C5 <- geary.test(petajabar$miskin,W4)  
C5
## 
##  Geary C test under randomisation
## 
## data:  petajabar$miskin 
## weights: W4 
## 
## Geary C statistic standard deviate = 4.1302, p-value = 1.812e-05
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.805313190       1.000000000       0.002221941
  1. Queen
C6 <- geary.test(petajabar$miskin,W4)  
C6
## 
##  Geary C test under randomisation
## 
## data:  petajabar$miskin 
## weights: W4 
## 
## Geary C statistic standard deviate = 4.1302, p-value = 1.812e-05
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.805313190       1.000000000       0.002221941

Simulasi Monte Carlo

#Dengan monte carlo:
set.seed(018)
C7 <- geary.mc(petajabar$miskin,W1, nsim=399)
C7
## 
##  Monte-Carlo simulation of Geary C
## 
## data:  petajabar$miskin 
## weights: W1 
## number of simulations + 1: 400 
## 
## statistic = 0.63708, observed rank = 2, p-value = 0.005
## alternative hypothesis: greater

Hasil Global Geary’s C Autocorrelation

##                   Bobot.Spatial Geary.s.C   p.value
## 1                    KNN Weight 0.6370830 5.799e-03
## 2        Radial Distance Weight 0.6653531 8.450e-04
## 3 Power Distance Weight alpha=1 0.8257360 1.952e-02
## 4 Power Distance Weight alpha=2 0.6608234 7.182e-03
## 5   Exponential Distance Weight 0.8053132 1.812e-05
## 6     Spatial Contiguity (Rook) 0.8053132 1.812e-05
## 7    Spatial Contiguity (Queen) 0.8053132 1.812e-05
## 8          Simulasi Monte Carlo 0.6370800 5.000e-03

Berdasarkan hasil di atas, diperoleh p-value < 0.05 berarti H0 ditolak, dengan alternative hypothesis: greater berarti pada taraf nyata 5%, terdapat autokorelasi spasial positif.

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(petajabar$miskin)
lisa <- localmoran(petajabar$miskin,W1)
lisa
##             Ii          E.Ii      Var.Ii       Z.Ii Pr(z != E(Ii))
## 1   0.40274292 -0.0034206161 0.027081939  2.4680903     0.01358360
## 2   0.03773872 -0.0034166067 0.027050305  0.2502303     0.80240924
## 3   0.02761853 -0.0161007157 0.125851779  0.1232375     0.90191900
## 4   0.36763087 -0.0128274967 0.100600119  1.1995211     0.23032539
## 5   0.48124464 -0.0257890053 0.199595686  1.1349094     0.25641327
## 6   0.44022459 -0.0130813970 0.102564955  1.4154416     0.15693906
## 7  -0.15545147 -0.0032722181 0.025910890 -0.9453971     0.34445619
## 8   1.00842208 -0.0514436416 0.387666591  1.7022434     0.08870975
## 9   1.10819674 -0.0743755166 0.546925737  1.5990554     0.10980830
## 10  0.74264273 -0.0572632690 0.428874375  1.2214450     0.22191759
## 11  0.14192119 -0.0061121635 0.048260950  0.6738475     0.50040826
## 12  1.06134396 -0.0810614425 0.591785520  1.4850398     0.13753329
## 13  0.34047226 -0.0168621909 0.131701867  0.9846433     0.32479932
## 14 -0.01870636 -0.0023289311 0.018458973 -0.1205430     0.90405304
## 15 -0.10441646 -0.0004750149 0.003771937 -1.6924137     0.09056712
## 16  1.49678574 -0.0722084307 0.532233077  2.1506533     0.03150358
## 17 -0.83483233 -0.0233510560 0.181179285 -1.9064448     0.05659252
## 18  0.84930716 -0.0185301894 0.144484193  2.2831170     0.02242348
## 19  0.03849214 -0.0046566449 0.036822186  0.2248608     0.82208755
## 20  0.50160978 -0.0937423249 0.674917906  0.7246836     0.46864618
## 21  0.13245000 -0.0004627854 0.003674871  2.1925308     0.02834120
## 22  1.63244773 -0.0665280966 0.493366755  2.4188128     0.01557125
## 23  1.63020958 -0.1868664910 1.207137722  1.6538446     0.09815909
## 24  0.52468954 -0.0557668757 0.418330065  0.8974496     0.36947905
## 25 -0.26765454 -0.1285172520 0.889782290 -0.1475033     0.88273480
## 26 -0.64588524 -0.0215396278 0.167434507 -1.5258175     0.12705533
## attr(,"call")
## localmoran(x = petajabar$miskin, listw = W1)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##         mean    median     pysal
## 1    Low-Low   Low-Low   Low-Low
## 2    Low-Low   Low-Low   Low-Low
## 3  High-High High-High High-High
## 4    Low-Low   Low-Low   Low-Low
## 5  High-High High-High High-High
## 6  High-High High-High High-High
## 7   Low-High  Low-High  Low-High
## 8  High-High High-High High-High
## 9  High-High High-High High-High
## 10 High-High High-High High-High
## 11 High-High High-High High-High
## 12 High-High High-High High-High
## 13 High-High High-High High-High
## 14  Low-High  Low-High  Low-High
## 15  High-Low  High-Low  High-Low
## 16   Low-Low   Low-Low   Low-Low
## 17  High-Low  High-Low  High-Low
## 18   Low-Low   Low-Low   Low-Low
## 19   Low-Low   Low-Low   Low-Low
## 20   Low-Low   Low-Low   Low-Low
## 21 High-High High-High High-High
## 22   Low-Low   Low-Low   Low-Low
## 23   Low-Low   Low-Low   Low-Low
## 24   Low-Low   Low-Low   Low-Low
## 25  High-Low  High-Low  High-Low
## 26  Low-High  Low-High  Low-High
petajabar$z.li <- lisa[,4]
petajabar$pvalue <- lisa[,5]
lm.palette <- colorRampPalette(c("yellow","orange", "red"), space = "rgb")
spplot(petajabar, zcol="z.li", col.regions=lm.palette(20), main="Local Moran")

Warna yang lebih pekat berarti memiliki nilai Z_Score yang besar, apabila Z-Score > 2 atau Z-score < -2 berarti daerah tersebut memiliki pengaruh yang signifikan terhadap tetangganya.

moran.plot(petajabar$miskin,W1)

Getis-Ord Gi

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

local_gi <- localG(petajabar$miskin,W1)

local_gi
##  [1] -2.4680903 -0.2502303  0.1232375 -1.1995211  1.1349094  1.4154416
##  [7]  0.9453971  1.7022434  1.5990554  1.2214450  0.6738475  1.4850398
## [13]  0.9846433  0.1205430 -1.6924137 -2.1506533 -1.9064448 -2.2831170
## [19] -0.2248608 -0.7246836  2.1925308 -2.4188128 -1.6538446 -0.8974496
## [25] -0.1475033  1.5258175
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = petajabar$miskin, listw = W1)
## attr(,"class")
## [1] "localG"

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

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

Notes

Pilih bobot yang menghasilkan autokorelasi spasial terbesar.