Spatial Neighbors and Dependency

Loading libraries and GIS data

rm(list=ls())
setwd("D:/_Lab_Data/R/R_Labs")
library(rgdal)
library(spdep)  

# Loading Shapefiles
TWN <- readOGR(dsn = "shapefiles", layer = "Popn_TWN", encoding="big5")

Study Area (Taipei City) for this session

Code<-TWN@data$Code
sel2<- Code < 42 & Code >=30
TWN_North<- TWN[sel2,]
plot(TWN_North); head(TWN_North@data)

##    UNI_ID   Town County Area Code towncode Male Female Age25 Age34 Age44
## 29     30 松山區 臺北市    9   30     6301 3148   2758   170   660  1109
## 30     31 信義區 臺北市   11   31     6302 4200   3402   242   934  1246
## 31     32 大安區 臺北市   11   32     6303 4036   3540   417   959  1313
## 32     33 中山區 臺北市   14   33     6304 4618   4553   439  1361  1808
## 33     34 中正區 臺北市    8   34     6305 4054   3243   417  1005  1174
## 34     35 大同區 臺北市    6   35     6306 3611   3020   189   856  1242
##    Age54 Age64 Age_L65 Popn
## 29  1392  1252    1323 5906
## 30  1795  1617    1768 7602
## 31  1604  1528    1755 7576
## 32  2123  1867    1573 9171
## 33  1577  1429    1695 7297
## 34  1506  1351    1487 6631

Spatial Neighbors

Contiguity: QUEEN vs. ROOK
K-nearest Neighbors (KNN)
Distance-based (fixed distance band)

1. Contiguity: QUEEN vs. ROOK

#?poly2nb: Construct neighbours list 
TWN_nbq<-poly2nb(TWN_North) #QUEEN = TRUE
summary(TWN_nbq)

## Neighbour list object:
## Number of regions: 12 
## Number of nonzero links: 48 
## Percentage nonzero weights: 33.33333 
## Average number of links: 4 
## Link number distribution:
## 
## 1 2 4 5 6 
## 1 1 6 3 1 
## 1 least connected region:
## 40 with 1 link
## 1 most connected region:
## 32 with 6 links

TWN_nbq2<-poly2nb(TWN_North,queen=FALSE )

1.1. Buiding Neighborhood Matrix

TWN_nbq_w.mat <-nb2mat(TWN_nbq, style="B"); TWN_nbq_w.mat

##    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## 29    0    1    1    1    0    0    0    0    1     1     0     0
## 30    1    0    1    0    0    0    0    1    1     0     0     0
## 31    1    1    0    1    1    0    0    1    0     0     0     0
## 32    1    0    1    0    1    1    0    0    0     1     1     0
## 33    0    0    1    1    0    1    1    1    0     0     0     0
## 34    0    0    0    1    1    0    1    0    0     0     1     0
## 35    0    0    0    0    1    1    0    0    0     0     0     0
## 36    0    1    1    0    1    0    0    0    1     0     0     0
## 37    1    1    0    0    0    0    0    1    0     1     0     0
## 38    1    0    0    1    0    0    0    0    1     0     1     0
## 39    0    0    0    1    0    1    0    0    0     1     0     1
## 40    0    0    0    0    0    0    0    0    0     0     1     0
## attr(,"call")
## nb2mat(neighbours = TWN_nbq, style = "B")

1.2. Finding neighbors of a district

TWN.region.id <- attr(TWN_nbq, "region.id")
TWN.neighbors.index = TWN_nbq[[match("33", TWN.region.id)]]
TWN.neighbors = rownames(TWN_nbq_w.mat[TWN.neighbors.index,]); TWN.neighbors

## [1] "31" "32" "34" "35" "36"

1.3. Plot the Neighborhood Matrix

coords<-coordinates(TWN_North)
plot(TWN_North)
plot(TWN_nbq, coords, add=T)

2. K-nearest Neighbors (KNN)

IDs <-TWN_North@data$UNI_ID
TWN_kn1<-knn2nb(knearneigh(coords, k=1), row.names=IDs)
plot(TWN_North)
plot(TWN_kn1, coords, add=T)
TWN_PTS<-SpatialPointsDataFrame(coords, TWN_North@data)
sel<- TWN_PTS@data$Area > 30
plot(TWN_PTS[sel,], col="blue", add=TRUE)

3. Distance-based (fixed distance band)

#?nbdists; ?dnearneigh
dist<-unlist(nbdists(TWN_kn1, coords)); summary(dist)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2249    2261    2752    3188    4167    4660

TWN_ran1<-dnearneigh(coords, d1=0, d2=7000, row.names=IDs)
plot(TWN_North)
plot(TWN_ran1, coords, add=T)

From Spatial Neighbors to ListW (W matrix)

# Row-standardized weights matrix
TWN_nbq_w<- nb2listw(TWN_nbq, zero.policy=T) # default: style = "W" 
TWN_nbq_w

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 12 
## Number of nonzero links: 48 
## Percentage nonzero weights: 33.33333 
## Average number of links: 4 
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 12 144 12 6.776667 49.63833

# Binary matrix
TWN_nbq_wb2<-nb2listw(TWN_nbq, style="B", zero.policy=T); TWN_nbq_wb2

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 12 
## Number of nonzero links: 48 
## Percentage nonzero weights: 33.33333 
## Average number of links: 4 
## 
## Weights style: B 
## Weights constants summary:
##    n  nn S0 S1  S2
## B 12 144 48 96 848

Global Analysis

Moran’s I
Getis-Ord General G

Moran’s I

1. Moran’s I Statistic

Popn<-TWN_North@data$Popn; Area<-TWN_North@data$Area
Density<-Popn/Area
TWN_North@data$Density<- Density
M<-moran.test(Density, listw=TWN_nbq_w, zero.policy=T); M

## 
##  Moran's I test under randomisation
## 
## data:  Density  
## weights: TWN_nbq_w  
## 
## Moran I statistic standard deviate = 2.4554, p-value = 0.007037
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.36409091       -0.09090909        0.03433907

2. Monte-Carlo simulation of Moran’s I

set.seed(123456)
bperm<-moran.mc(Density,listw=TWN_nbq_w,nsim=999)
hist(bperm$res, freq=TRUE, breaks=20, xlab="Simulated Moran's I")
abline(v=0.39519718, col="red")

3. Spatial correlogram

cor2<-sp.correlogram(TWN_nbq, Density, order=2, method="I", style="W")
print(cor2); plot(cor2)

## Spatial correlogram for Density 
## method: Moran's I
##         estimate expectation  variance standard deviate Pr(I) two sided  
## 1 (12)  0.364091   -0.090909  0.034339           2.4554         0.01407 *
## 2 (12) -0.305631   -0.090909  0.024241          -1.3791         0.16786  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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 = 3.264, p-value = 0.0005492
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##        0.747403137        0.500000000        0.005745246

Local Analysis

Local Moran’s I (LISA)
Local G: Gi(d)

Local Moran’s I (LISA)

1. LISA statistic

LISA.Density <- localmoran(Density, TWN_nbq_w, zero.policy=T)
LISA.Density2 <- as.data.frame(localmoran.exact(lm(Density ~1, TWN_North@data), nb=TWN_nbq, style="C"))
class(LISA.Density)

## [1] "localmoran" "matrix"

printCoefmat(data.frame(LISA.Density))

##           Ii      E.Ii    Var.Ii      Z.Ii Pr.z...0.
## 29 -0.032517 -0.090909  0.116304  0.171221    0.4320
## 30 -0.107933 -0.090909  0.166793 -0.041683    0.5166
## 31  0.122681 -0.090909  0.116304  0.626303    0.2656
## 32  0.084420 -0.090909  0.082645  0.609884    0.2710
## 33  0.560185 -0.090909  0.116304  1.909179    0.0281
## 34  0.460745 -0.090909  0.166793  1.350760    0.0884
## 35  1.228949 -0.090909  0.419236  2.038439    0.0208
## 36 -0.347543 -0.090909  0.166793 -0.628384    0.7351
## 37  0.217835 -0.090909  0.166793  0.755980    0.2248
## 38  0.335753 -0.090909  0.166793  1.044709    0.1481
## 39  0.076730 -0.090909  0.166793  0.410475    0.3407
## 40  1.769785 -0.090909  0.924124  1.935574    0.0265

LISA.Density[,5] # p-value

##         29         30         31         32         33         34 
## 0.43202490 0.51662448 0.26555822 0.27096920 0.02811951 0.08838617 
##         35         36         37         38         39         40 
## 0.02075302 0.73512385 0.22483046 0.14807871 0.34072870 0.02645992

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

spplot(TWN_North, zcol="z.li", col.regions=lm.palette(20), main="Local Moran")

head(TWN_North@data, n=12)

##    UNI_ID        Town County Area Code towncode Male Female Age25 Age34
## 29     30      松山區 臺北市    9   30     6301 3148   2758   170   660
## 30     31      信義區 臺北市   11   31     6302 4200   3402   242   934
## 31     32      大安區 臺北市   11   32     6303 4036   3540   417   959
## 32     33      中山區 臺北市   14   33     6304 4618   4553   439  1361
## 33     34      中正區 臺北市    8   34     6305 4054   3243   417  1005
## 34     35      大同區 臺北市    6   35     6306 3611   3020   189   856
## 35     36 萬\xb5?\xcf 臺北市    9   36     6307 4090   3188   164   737
## 36     37      文山區 臺北市   32   37     6308 3818   3173   326   766
## 37     38      南港區 臺北市   22   38     6309 4001   2844   136   849
## 38     39      內湖區 臺北市   32   39     6310 3933   3093   202   981
## 39     40      士林區 臺北市   62   40     6311 4394   3253   442   777
## 40     41      北投區 臺北市   57   41     6312 3805   2675   210   632
##    Age44 Age54 Age64 Age_L65 Popn   Density        z.li     pvalue
## 29  1109  1392  1252    1323 5906  656.2222  0.17122125 0.43202490
## 30  1246  1795  1617    1768 7602  691.0909 -0.04168347 0.51662448
## 31  1313  1604  1528    1755 7576  688.7273  0.62630266 0.26555822
## 32  1808  2123  1867    1573 9171  655.0714  0.60988439 0.27096920
## 33  1174  1577  1429    1695 7297  912.1250  1.90917894 0.02811951
## 34  1242  1506  1351    1487 6631 1105.1667  1.35075995 0.08838617
## 35  1238  1759  1586    1794 7278  808.6667  2.03843907 0.02075302
## 36  1485  1784  1365    1265 6991  218.4688 -0.62838418 0.73512385
## 37  1385  1725  1413    1337 6845  311.1364  0.75598043 0.22483046
## 38  1550  2061  1396     836 7026  219.5625  1.04470912 0.14807871
## 39  1415  1868  1660    1485 7647  123.3387  0.41047519 0.34072870
## 40  1187  1705  1526    1220 6480  113.6842  1.93557450 0.02645992

writeOGR(TWN_North, dsn="shapefiles", layer="Popn_LISA", driver="ESRI Shapefile")

2. Moran Scatter Plot

nci <- moran.plot (Density, TWN_nbq_w, labels=IDs , xlab="Popn Density", ylab="SL Popn Density")

3. LISA map

chk<-Density-mean(Density)
zi<- LISA.Density[,4]
quadrant <- vector(mode="numeric",length=nrow(LISA.Density))
quadrant[chk>0 & zi>0] <- 1 # H-H
quadrant[chk<0 & zi>0] <- 2 # L-L
quadrant[chk>0 & zi<0] <- 3 # H-L
quadrant[chk<0 & zi<0] <- 4 # L-H

signif <- 0.05
quadrant[LISA.Density[, 5]> signif] <- 5
colors <- c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
par(mar=c(0,0,1,0)) # sets margin parameters for plot space; 
# A numeric vector of length 4, which sets the margin sizes in the following order: bottom, left, top, and right. The default is c(5.1, 4.1, 4.1, 2.1).
plot(TWN_North, border="grey", col=colors[quadrant], main = "LISA Cluster Map, Population Density")
legend("bottomright",legend=c("High-High","Low-Low","High-Low","Low-High"), fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)

#dev.off()

Adjustment for p-value

TWN_North$p.adj <- p.adjust(LISA.Density[,5], method="bonferroni")
head(TWN_North@data[,16:19], n=12)

##      Density        z.li     pvalue     p.adj
## 29  656.2222  0.17122125 0.43202490 1.0000000
## 30  691.0909 -0.04168347 0.51662448 1.0000000
## 31  688.7273  0.62630266 0.26555822 1.0000000
## 32  655.0714  0.60988439 0.27096920 1.0000000
## 33  912.1250  1.90917894 0.02811951 0.3374341
## 34 1105.1667  1.35075995 0.08838617 1.0000000
## 35  808.6667  2.03843907 0.02075302 0.2490363
## 36  218.4688 -0.62838418 0.73512385 1.0000000
## 37  311.1364  0.75598043 0.22483046 1.0000000
## 38  219.5625  1.04470912 0.14807871 1.0000000
## 39  123.3387  0.41047519 0.34072870 1.0000000
## 40  113.6842  1.93557450 0.02645992 0.3175190

Local G

TWN_nbq_in<-include.self(TWN_nbq); summary(TWN_nbq_in)

## Neighbour list object:
## Number of regions: 12 
## Number of nonzero links: 60 
## Percentage nonzero weights: 41.66667 
## Average number of links: 5 
## Link number distribution:
## 
## 2 3 5 6 7 
## 1 1 6 3 1 
## 1 least connected region:
## 40 with 2 links
## 1 most connected region:
## 32 with 7 links

TWN_nbq_in_w<- nb2listw(TWN_nbq_in, zero.policy=T)
LG<-localG(Density, TWN_nbq_in_w); LG

##  [1] -0.05179141 -0.25373181  0.99011884  0.99812167  1.97406323
##  [6]  1.57593542  2.40689738  0.19703028 -1.08014601 -1.31116043
## [11] -0.86816724 -1.97333072
## attr(,"gstari")
## [1] TRUE
## attr(,"call")
## localG(x = Density, listw = TWN_nbq_in_w)
## attr(,"class")
## [1] "localG"

class(LG)

## [1] "localG"

LG1<-0
for (i in 1:12){LG1[i]<-LG[i]}
TWN_North$LG<-LG1
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(TWN_North, zcol="LG", col.regions=lm.palette(20), main="Local G")

Spatial Neighbors and Dependency

Loading libraries and GIS data

Study Area (Taipei City) for this session

Spatial Neighbors

Contiguity: QUEEN vs. ROOK

K-nearest Neighbors (KNN)

Distance-based (fixed distance band)

1. Contiguity: QUEEN vs. ROOK

1.1. Buiding Neighborhood Matrix

1.2. Finding neighbors of a district

1.3. Plot the Neighborhood Matrix

2. K-nearest Neighbors (KNN)

3. Distance-based (fixed distance band)

From Spatial Neighbors to ListW (W matrix)

Global Analysis

Moran’s I

Getis-Ord General G

Moran’s I

1. Moran’s I Statistic

2. Monte-Carlo simulation of Moran’s I

3. Spatial correlogram

Getis-Ord General G Statistic

Local Analysis

Local Moran’s I (LISA)

Local G: Gi(d)

Local Moran’s I (LISA)

1. LISA statistic

2. Moran Scatter Plot

3. LISA map

Adjustment for p-value

Local G