Spatial Analysis: 10

Spatial Autocorrelation (categorical variables)

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\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,]  #研究區：台北＋新北

1. Defining Neighborhood

# 定義鄰近關係與鄰近矩陣
TWN_nb <- poly2nb(NorthTW_sf) #QUEEN = TRUE
TWN_nb_w <- nb2listw(TWN_nb, zero.policy=T)

2. Creating categorical variables

# 計算切點：連續變數 轉換 類別變數分類
ave <- mean(NorthTW_sf$A0A14_CNT)
max <- max(NorthTW_sf$A0A14_CNT)

HIPOP <- cut(NorthTW_sf$A0A14_CNT, breaks=c(0,ave,max), labels=c("LOW","HIGH"))
names(HIPOP) <- rownames(NorthTW_sf)

3-1. Join Counts Statistics: two categories

joincount.test(HIPOP, TWN_nb_w)

    Join count test under nonfree sampling

data:  HIPOP 
weights: TWN_nb_w 

Std. deviate for LOW = 3.2383, p-value = 0.0006013
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            7.7702381             5.7750000             0.3796304 


    Join count test under nonfree sampling

data:  HIPOP 
weights: TWN_nb_w 

Std. deviate for HIGH = 4.4985, p-value = 3.421e-06
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            6.9260823             4.2750000             0.3473002 

joincount.multi(HIPOP, TWN_nb_w)

          Joincount Expected Variance z-value
LOW:LOW     7.77024  5.77500  0.37963  3.2383
HIGH:HIGH   6.92608  4.27500  0.34730  4.4985
HIGH:LOW    5.80368 10.45000  1.03357 -4.5702
Jtot        5.80368 10.45000  1.03357 -4.5702

3-2. Join Counts Statistics: more categories

# 兩個以上的類別
cut1 <- quantile(NorthTW_sf$A0A14_CNT,0.33)
cut2 <- quantile(NorthTW_sf$A0A14_CNT,0.66)
max <- max(NorthTW_sf$A0A14_CNT)

HIPOP <- cut(NorthTW_sf$A0A14_CNT, breaks=c(0, cut1, cut2, max), 
             labels=c("LOW","MID","HIGH"))
names(HIPOP) <- rownames(NorthTW_sf)

joincount.test(HIPOP, TWN_nb_w)

    Join count test under nonfree sampling

data:  HIPOP 
weights: TWN_nb_w 

Std. deviate for LOW = 4.3537, p-value = 6.692e-06
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            4.4547619             2.2750000             0.2506667 


    Join count test under nonfree sampling

data:  HIPOP 
weights: TWN_nb_w 

Std. deviate for MID = 0.81447, p-value = 0.2077
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            2.3380952             1.9500000             0.2270534 


    Join count test under nonfree sampling

data:  HIPOP 
weights: TWN_nb_w 

Std. deviate for HIGH = 3.3322, p-value = 0.0004308
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            3.9433261             2.2750000             0.2506667 

joincount.multi(HIPOP, TWN_nb_w)

          Joincount Expected Variance z-value
LOW:LOW     4.45476  2.27500  0.25067  4.3537
MID:MID     2.33810  1.95000  0.22705  0.8145
HIGH:HIGH   3.94333  2.27500  0.25067  3.3322
MID:LOW     2.71548  4.55000  0.54680 -2.4809
HIGH:LOW    2.25364  4.90000  0.57924 -3.4771
HIGH:MID    4.79470  4.55000  0.54680  0.3309
Jtot        9.76382 14.00000  0.90287 -4.4582