Spatial Analysis: 09

Spatial Autocorrelation

0. Loading R packages and data

rm(list = ls())
library(spdep)
library(sf)
library(tmap)

# setwd("C:/Data")
TWPOP_sf <- st_read("./Data/Popn_TWN2.shp",options = "encoding=Big5" )

options:        encoding=Big5 
Reading layer `Popn_TWN2' from data source `C:\Users\wenth\Desktop\WK09\Data\Popn_TWN2.shp' using driver `ESRI Shapefile'
Simple feature collection with 368 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -197572.5 ymin: 2295201 xmax: 606875.6 ymax: 2919551
Projected CRS: TWD97 / TM2 zone 121

ID<- TWPOP_sf$COUNTY_ID
Sel<- ID == "65000" | ID == "63000" 
NorthTW_sf<- TWPOP_sf[Sel,]
length(NorthTW_sf)

[1] 15

NorthTW_lyr <- tm_shape(NorthTW_sf)+tm_polygons(col="grey")

1. Define Spatial Neighbors

a. Contiguity: QUEEN vs. ROOK

#?poly2nb: Construct neighbours list 
TWN_nb<-poly2nb(NorthTW_sf) #QUEEN = TRUE
summary(TWN_nb)

Neighbour list object:
Number of regions: 41 
Number of nonzero links: 202 
Percentage nonzero weights: 12.01666 
Average number of links: 4.926829 
Link number distribution:

 2  3  4  5  6  7  8  9 11 
 3  7  6 12  7  3  1  1  1 
3 least connected regions:
278 292 296 with 2 links
1 most connected region:
231 with 11 links

TWN_nb2<-poly2nb(NorthTW_sf,queen=FALSE ) # ROOK

# 1.1. Finding neighbors of a district
TWN_nb[1]

[[1]]
[1]  2  3  4  9 10

# 1-2 Buiding Neighborhood Matrix
TWN_nb_w.mat <- nb2mat(TWN_nb, style="B")
# style = B is the basic binary coding, 
# W is row standardised, C is globally standardised

# 1.3. Plot the Neighborhood Matrix
coords_sf<-st_centroid(NorthTW_sf)
coords_lyr <- tm_shape(coords_sf)+tm_dots(col="red", size=0.3)
NorthTW_lyr + coords_lyr

TW.centroid <- st_centroid(NorthTW_sf)
TW.coords <- st_coordinates(TW.centroid)
TW.net <- nb2lines(TWN_nb,coords=TW.coords)
TW.net_sf <- st_as_sf(TW.net)
st_crs(TW.net_sf) <-  st_crs(TWPOP_sf)

NorthTW_lyr + coords_lyr + 
  tm_shape(TW.net_sf) + tm_lines(col='red')

b. K-nearest Neighbors (KNN)

IDs <-NorthTW_sf$TOWN_ID
TWN_kn1<-knn2nb(knearneigh(TW.coords, k=2), row.names=IDs)
TWN_kn1[1]

[[1]]
[1] 2 4

TW.net <- nb2lines(TWN_kn1,coords=TW.coords)
TW.net_sf <- st_as_sf(TW.net)
st_crs(TW.net_sf) <-  st_crs(TWPOP_sf)

NorthTW_lyr + coords_lyr + 
  tm_shape(TW.net_sf) + tm_lines(col='red')

c. Fixed distance band

TWN_ran1<-dnearneigh(TW.coords, d1=0, d2=20000, row.names=IDs)
TW.net <- nb2lines(TWN_ran1,coords=TW.coords)
TW.net_sf <- st_as_sf(TW.net)
st_crs(TW.net_sf) <-  st_crs(TWPOP_sf)

NorthTW_lyr + coords_lyr + 
  tm_shape(TW.net_sf) + tm_lines(col='red')

2. Weighting matrix

# From Spatial Neighbors to ListW 

TWN_nb_w<- nb2listw(TWN_nb, zero.policy=T) # default: style = "W" (row standardised)
TWN_nb_w$weight[1]

[[1]]
[1] 0.2 0.2 0.2 0.2 0.2

3. Moran’s I & General G Statistics

Popn<-NorthTW_sf$A65UP_CNT
NorthTW_sf$Density<- Popn * 10^6 / st_area(NorthTW_sf)


# 1. Mapping the attribute
tm_shape(NorthTW_sf) + 
  tm_polygons("Density", palette = "OrRd", style = "jenks", title = "Elder Density")

# 2.Moran I Statistic
Density <- NorthTW_sf$Density
Density <- as.vector(Density)
M <- moran.test(Density, listw=TWN_nb_w, zero.policy=T)

# 3.Monte-Carlo simulation of Moran I
bperm <- moran.mc(Density,listw=TWN_nb_w,nsim=999)
hist(bperm$res, freq=TRUE, breaks=20, xlab="Simulated Moran's I")
abline(v=0.486, col="red")

# 4.Moran Spatial Correlogram
cor<-sp.correlogram(TWN_nb, Density, order=3, method="I", style="W")
print(cor); plot(cor)

Spatial correlogram for Density 
method: Moran's I
         estimate expectation   variance standard deviate Pr(I) two sided    
1 (41)  0.6327816  -0.0250000  0.0094261           6.7751       1.243e-11 ***
2 (41)  0.1797066  -0.0250000  0.0044722           3.0611        0.002206 ** 
3 (41) -0.1876583  -0.0250000  0.0031154          -2.9142        0.003566 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# 5.Moran Scatter Plot
nci <- moran.plot (Density, TWN_nb_w, labels=IDs , xlab="Popn", ylab="SL Popn")

# 6.Getis-Ord General G Statistic
TWN_ran1_wb <- nb2listw(TWN_ran1, style="B", zero.policy=T)
G <- globalG.test(Density, listw=TWN_ran1_wb)
G

    Getis-Ord global G statistic

data:  Density 
weights: TWN_ran1_wb   

standard deviate = 4.813, p-value = 7.434e-07
alternative hypothesis: greater
sample estimates:
Global G statistic        Expectation           Variance 
       0.969818340        0.540243902        0.007966092