使用R進行空間自相關檢驗

 

「全局溢出」當一個區域的特征變化影響到所有區域的結果時，就會產生全局溢出效應。這甚至適用於區域本身，因為影響可以傳遞到鄰居並返回到自己的區域（反饋）。具體來說，全球溢出效應影響到鄰居、鄰居到鄰居、鄰居到鄰居等等。
「局部溢出」是指影響只落在附近或近鄰的情況，在它們影響鄰鄰區域之前就消失了。

對應全局與局部溢出，存在全局與局部自相關檢驗。全局自相關檢驗指標主要有 moran'I 指數、Geary 指數 C 統計量以及 Getis-Ord global G 統計量；局部自相檢驗指標主要有局部moran指數、Local Getis-Ord Gi 和 Gi* 統計量。此外，也可以通過圖示的方法來顯示空間的依賴關系。

下面將對相關指標的計算以及可視化的方法進行演示(所用數據均可通過公共號后台回復「地圖可視化R」獲取)。

全局空間自相關檢驗

首先導入數據和矩陣，並進行適當轉換

library(readxl)
library("spdep")
# 設置工作路徑
setwd('E:\\空間計量專題\\R-空間計量')
# 導入經濟變量數據
cdata <- read_xlsx("E:/空間計量專題/R-空間計量/cdata.xlsx")
# 導入自定義矩陣並做適當格式轉換
w1 <- read_xlsx("w1.xlsx", col_names = FALSE)
w2 <- as.matrix(w1)
w2[1:5, 1:5]
# 轉換格式並標准化
w <- mat2listw(w2, style="W")

導入shp文件，合並數據

library(rgdal)
# 導入shp文件
shpt <- readOGR("廣東地級市.shp")
# 合並數據
cdatashpt <- merge(shpt, cdata, by = "city")

moran'I 指數

1.Moran’s I 統計量的計算

>>> moran(cdatashpt$gdp2017, listw=w, n=length(cdatashpt$gdp2017), S0=Szero(w))

$I
0.065546870128173
$K
2.99333822437931

I Moran’s I 統計量, K 變量的峰度

當數據中出現孤島或者缺失值時，可通過下面的子選項進行調整：

zero.policy默認值為空，使用全局選項值；如果為真，則將0分配給沒有鄰居的區域的滯后值，如果為假，則分配為NA

NAOK如果 TRUE，則 x 中的任何 NA 或 NaN 或 Inf 值都將傳遞給外部函數。如果 FALSE ，則存在 NA 或 NaN 或 Inf 值將被視為錯誤。

2.Monte-Carlo simulation of Moran's I

>>> # Monte-Carlo simulation of Moran's I
>>> set.seed(12345)
>>> moran.mc(cdatashpt$gdp2017, listw = w, nsim = 999, alternative = 'greater')

 Monte-Carlo simulation of Moran I

data: cdatashpt$gdp2017 
weights: w 
number of simulations + 1: 1000 

statistic = 0.065547, observed rank = 790, p-value = 0.21
alternative hypothesis: greater

3.moran散點圖

>>> # Moran 散點圖
>>> moran.plot(cdatashpt$gdp2017, w, zero.policy=NULL, spChk=NULL, labels=TRUE, xlab=NULL, ylab=NULL, quiet=NULL)

labels給具有高影響度量值的點添加字符標簽，如果設置為FALSE，則不會為具有較大影響的點繪制標簽

4.Moran’s I 空間自相關檢驗

>>> # Moran’s I test for spatial autocorrelation
>>> moran.test(cdatashpt$gdp2017, w)

 Moran I test under randomisation

data: cdatashpt$gdp2017 
weights: w 

Moran I statistic standard deviate = 0.76045, p-value = 0.2235
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance 
 0.06554687 -0.05000000 0.02308712 

Geary 檢驗

>>> geary.test(cdatashpt$gdp2017, listw=w, randomisation=TRUE, alternative="greater")

 Geary C test under randomisation

data: cdatashpt$gdp2017 
weights: w 

Geary C statistic standard deviate = 1.2898, p-value = 0.09856
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance 
 0.79614300 1.00000000 0.02498145 

Geary’s C 檢驗

1.計算 Geary’s C 統計量

>>> geary(cdatashpt$gdp2017, listw=w, n=length(w), n1=length(w)-1, S0=Szero(w))

$C
0.0796142995607693
$K
0.427619746339901

2.Monte-Carlo simulation of Geary's C

>>> geary.mc(cdatashpt$gdp2017, listw=w, nsim=999, alternative="greater")

 Monte-Carlo simulation of Geary C

data: cdatashpt$gdp2017 
weights: w 
number of simulations + 1: 1000 

statistic = 0.79614, observed rank = 116, p-value = 0.116
alternative hypothesis: greater

3.模擬分布圖

>>> set.seed(12345)
>>> gdpgeary <- geary.mc(cdatashpt$gdp2017, listw=w, nsim=999, alternative="greater")
>>> plot(gdpgeary, type='l', col='orange')
>>> gdpgeary.dens = density(gdpgeary$res)
>>> polygon(gdpgeary.dens, col="gray")
>>> abline(v=gdpgeary$statistic, col='orange',lwd=2)

Getis-Ord global G 檢驗

>>> globalG.test(cdatashpt$gdp2017, listw=w, alternative="greater")

 Getis-Ord global G statistic

data: cdatashpt$gdp2017 
weights: w 

standard deviate = -1.1248, p-value = 0.8697
alternative hypothesis: greater
sample estimates:
Global G statistic Expectation Variance 
 4.948903e-02 5.000000e-02 2.063625e-07 

空間相關圖

>>> w.nb <- w$neighbours
>>> spcorrI = sp.correlogram(w.nb, cdatashpt$gdp2017, order = 2, method = "I", style = "W", randomisation = TRUE)
>>> spcorrI

Spatial correlogram for cdatashpt$gdp2017 
method: Moran's I
 estimate expectation variance standard deviate Pr(I) two sided
1 (21) 0.065547 -0.050000 0.023087 0.7605 0.4470
2 (21) -0.130212 -0.050000 0.017101 -0.6134 0.5396

>>> plot(spcorrI, main="Spatial correlogram of gdp2017")

局部空間自相關檢驗

局部moran檢驗

>>> localmoran(cdatashpt$gdp2017, listw=w, alternative = "greater")

A localmoran: 21 × 5 of type dbl
 Ii E.Ii Var.I Z.Ii Pr(z > 0)
1 -0.01370494 -0.05 0.1146316 0.107200142 0.4573151
2 0.65060843 -0.05 0.4279122 1.071021124 0.1420800
3 -0.32033200 -0.05 0.4279122 -0.413256930 0.6602908
4 0.01710619 -0.05 0.4279122 0.102585331 0.4591460
5 0.05648861 -0.05 0.1146316 0.314522027 0.3765623
6 0.48390574 -0.05 0.1929517 1.215458873 0.1120956
7 -0.48582426 -0.05 0.1929517 -0.992172269 0.8394433
8 0.10421133 -0.05 0.1929517 0.351068574 0.3627685
9 -0.26267734 -0.05 0.1146316 -0.628158331 0.7350499
10 -0.44651083 -0.05 0.1929517 -0.902673594 0.8166504
11 -0.38816462 -0.05 0.2712719 -0.649270674 0.7419183
12 0.02013525 -0.05 0.1929517 0.159665854 0.4365722
13 -0.03877614 -0.05 0.1459596 0.029378254 0.4882815
14 0.41014395 -0.05 0.2712719 0.883469050 0.1884914
15 0.02782689 -0.05 0.8978331 0.082135687 0.4672694
16 0.11878306 -0.05 0.2712719 0.324060781 0.3729460
17 0.56506605 -0.05 0.2712719 1.180917005 0.1188178
18 0.47054422 -0.05 0.1929517 1.185040809 0.1180007
19 -0.05142908 -0.05 0.2712719 -0.002743819 0.5010946
20 0.06940159 -0.05 0.1929517 0.271822740 0.3928792
21 0.38968220 -0.05 0.1459596 1.150860065 0.1248949

Local Getis-Ord Gi 和 Gi*統計量

>>> localG(cdatashpt$gdp2017, listw=w)

 [1] -0.05909956 1.09552796 -0.05815482 -0.09536069 -1.69368469 -1.94012537
 [7] 0.69295805 0.34995778 0.96964875 -0.34538087 -0.15784709 0.51802708
[13] 0.20826225 -0.87609799 -0.19862319 -1.09970094 -1.06316301 -1.42853393
[19] 0.38484915 1.17379925 -1.43159536
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = cdatashpt$gdp2017, listw = w)
attr(,"class")
[1] "localG"

基於殘差項的moran檢驗

>>> ols <- lm(gdp2017 ~ kj2017 + l2017 + ks2017 + pe2017 + inex2017 + new_inc2017 + pri_en2017 + high_stu2017, data = cdatashpt)
>>> lm.morantest(ols, listw = w, alternative = "two.sided")

 Global Moran I for regression residuals

data: 
model: lm(formula = gdp2017 ~ kj2017 + l2017 + ks2017 + pe2017 +
inex2017 + new_inc2017 + pri_en2017 + high_stu2017, data = cdatashpt)
weights: w

Moran I statistic standard deviate = 0.95884, p-value = 0.3376
alternative hypothesis: two.sided
sample estimates:
Observed Moran I Expectation Variance 
 0.001873225 -0.123283278 0.017037907 

散點圖的繪制

>>> moran.plot(ols$residuals, w)

局部moran檢驗

localmoran(ols$residuals, w)

A localmoran: 21 × 5 of type dbl
 Ii E.Ii Var.I Z.Ii Pr(z > 0)
1 -0.0363243420 -0.05 0.1168494 0.040006912 0.48404381
2 0.8370573121 -0.05 0.4404798 1.336560700 0.09068304
3 -1.1222105809 -0.05 0.4404798 -1.615537694 0.94690285
4 0.2373930487 -0.05 0.4404798 0.433025295 0.33249820
5 -0.1380790378 -0.05 0.1168494 -0.257667335 0.60166817
6 -0.1081663716 -0.05 0.1977570 -0.130799494 0.55203304
7 0.1159985869 -0.05 0.1977570 0.373283232 0.35446883
8 0.7082098711 -0.05 0.1977570 1.704996632 0.04409753
9 0.0612709718 -0.05 0.1168494 0.325513260 0.37239632
10 0.4204181464 -0.05 0.1977570 1.057835549 0.14506521
11 0.3774028474 -0.05 0.2786646 0.809648516 0.20907111
12 -0.0868802622 -0.05 0.1977570 -0.082933137 0.53304765
13 0.1244368071 -0.05 0.1492124 0.451580986 0.32578544
14 0.1874982314 -0.05 0.2786646 0.449903626 0.32638997
15 -0.5517955762 -0.05 0.9259254 -0.521481405 0.69898427
16 -0.1211737778 -0.05 0.2786646 -0.134827702 0.55362595
17 -0.0030386117 -0.05 0.2786646 0.088961079 0.46455642
18 -0.4168336708 -0.05 0.1977570 -0.824903760 0.79528688
19 -0.0539797929 -0.05 0.2786646 -0.007539102 0.50300764
20 -0.3916838150 -0.05 0.1977570 -0.768348944 0.77886005
21 -0.0001822581 -0.05 0.1492124 0.128967879 0.44869153

其他相關檢驗類似，只需將變量替換為殘差項即可，不再演示。此外，這里僅演示了各命令的常見用法，需要具體了解的話，可以點擊「閱讀原文」，這里給出了各命令的詳細說明。

 

 

更多精彩內容，掃碼關注微信公眾號！