moran指数 r语言_使用R进行空间自相关检验

「全局溢出」当一个区域的特征变化影响到所有区域的结果时，就会产生全局溢出效应。这甚至适用于区域本身，因为影响可以传递到邻居并返回到自己的区域(反馈)。具体来说，全球溢出效应影响到邻居、邻居到邻居、邻居到邻居等等。

「局部溢出」是指影响只落在附近或近邻的情况，在它们影响邻邻区域之前就消失了。

对应全局与局部溢出，存在全局与局部自相关检验。全局自相关检验指标主要有 moran'I 指数、Geary 指数 C 统计量以及 Getis-Ord global G 统计量；局部自相检验指标主要有局部moran指数、Local Getis-Ord Gi 和 Gi* 统计量。此外，也可以通过图示的方法来显示空间的依赖关系。

下面将对相关指标的计算以及可视化的方法进行演示(所用数据均可通过公共号后台回复「地图可视化R」获取)。

全局空间自相关检验

首先导入数据和矩阵，并进行适当转换

library(readxl)

library("spdep")

# 设置工作路径

setwd('E:\\空间计量专题\\R-空间计量')

# 导入经济变量数据

cdata

# 导入自定义矩阵并做适当格式转换

w1

w2

w2[1:5, 1:5]

# 转换格式并标准化

w

导入shp文件，合并数据

library(rgdal)

# 导入shp文件

shpt

# 合并数据

cdatashpt

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

>>> 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

>>> 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.050.11463160.1072001420.4573151

20.65060843-0.050.42791221.0710211240.1420800

3-0.32033200-0.050.4279122-0.4132569300.6602908

40.01710619-0.050.42791220.1025853310.4591460

50.05648861-0.050.11463160.3145220270.3765623

60.48390574-0.050.19295171.2154588730.1120956

7-0.48582426-0.050.1929517-0.9921722690.8394433

80.10421133-0.050.19295170.3510685740.3627685

9-0.26267734-0.050.1146316-0.6281583310.7350499

10-0.44651083-0.050.1929517-0.9026735940.8166504

11-0.38816462-0.050.2712719-0.6492706740.7419183

120.02013525-0.050.19295170.1596658540.4365722

13-0.03877614-0.050.14595960.0293782540.4882815

140.41014395-0.050.27127190.8834690500.1884914

150.02782689-0.050.89783310.0821356870.4672694

160.11878306-0.050.27127190.3240607810.3729460

170.56506605-0.050.27127191.1809170050.1188178

180.47054422-0.050.19295171.1850408090.1180007

19-0.05142908-0.050.2712719-0.0027438190.5010946

200.06940159-0.050.19295170.2718227400.3928792

210.38968220-0.050.14595961.1508600650.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.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.050.11684940.0400069120.48404381

20.8370573121-0.050.44047981.3365607000.09068304

3-1.1222105809-0.050.4404798-1.6155376940.94690285

40.2373930487-0.050.44047980.4330252950.33249820

5-0.1380790378-0.050.1168494-0.2576673350.60166817

6-0.1081663716-0.050.1977570-0.1307994940.55203304

70.1159985869-0.050.19775700.3732832320.35446883

80.7082098711-0.050.19775701.7049966320.04409753

90.0612709718-0.050.11684940.3255132600.37239632

100.4204181464-0.050.19775701.0578355490.14506521

110.3774028474-0.050.27866460.8096485160.20907111

12-0.0868802622-0.050.1977570-0.0829331370.53304765

130.1244368071-0.050.14921240.4515809860.32578544

140.1874982314-0.050.27866460.4499036260.32638997

15-0.5517955762-0.050.9259254-0.5214814050.69898427

16-0.1211737778-0.050.2786646-0.1348277020.55362595

17-0.0030386117-0.050.27866460.0889610790.46455642

18-0.4168336708-0.050.1977570-0.8249037600.79528688

19-0.0539797929-0.050.2786646-0.0075391020.50300764

20-0.3916838150-0.050.1977570-0.7683489440.77886005

21-0.0001822581-0.050.14921240.1289678790.44869153

其他相关检验类似，只需将变量替换为残差项即可，不再演示。此外，这里仅演示了各命令的常见用法，需要具体了解的话，可以点击「阅读原文」，这里给出了各命令的详细说明。

更多精彩内容，扫码关注微信公众号！