「全局溢出」当一个区域的特征变化影响到所有区域的结果时, 就会产生全局溢出效应. 这甚至适用于区域本身, 因为影响可以传递到邻居并返回到自己的区域 (反馈). 具体来说, 全球溢出效应影响到邻居, 邻居到邻居, 邻居到邻居等等.
「局部溢出」是指影响只落在附近或近邻的情况, 在它们影响邻邻区域之前就消失了.
对应全局与局部溢出, 存在全局与局部自相关检验. 全局自相关检验指标主要有 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
其他相关检验类似, 只需将变量替换为残差项即可, 不再演示. 此外, 这里仅演示了各命令的常见用法, 需要具体了解的话, 可以点击「阅读原文」, 这里给出了各命令的详细说明.
来源: https://www.cnblogs.com/xgzhu/p/12449965.html