使用R进行空间自相关检验

您所在的位置:网站首页 莫兰指数在r中的实现 使用R进行空间自相关检验

使用R进行空间自相关检验

2023-11-29 06:53| 来源: 网络整理| 查看: 265

 

「全局溢出」当一个区域的特征变化影响到所有区域的结果时,就会产生全局溢出效应。这甚至适用于区域本身,因为影响可以传递到邻居并返回到自己的区域(反馈)。具体来说,全球溢出效应影响到邻居、邻居到邻居、邻居到邻居等等。 「局部溢出」是指影响只落在附近或近邻的情况,在它们影响邻邻区域之前就消失了。

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

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

全局空间自相关检验

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

library(readxl)library("spdep")# 设置工作路径setwd('E:\\空间计量专题\\R-空间计量')# 导入经济变量数据cdata >> # 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.21alternative 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.2235alternative hypothesis: greatersample 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.09856alternative hypothesis: Expectation greater than statisticsample 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))

$C0.0796142995607693$K0.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.116alternative 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.8697alternative hypothesis: greatersample 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 sided1 (21) 0.065547 -0.050000 0.023087 0.7605 0.44702 (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.45731512 0.65060843 -0.05 0.4279122 1.071021124 0.14208003 -0.32033200 -0.05 0.4279122 -0.413256930 0.66029084 0.01710619 -0.05 0.4279122 0.102585331 0.45914605 0.05648861 -0.05 0.1146316 0.314522027 0.37656236 0.48390574 -0.05 0.1929517 1.215458873 0.11209567 -0.48582426 -0.05 0.1929517 -0.992172269 0.83944338 0.10421133 -0.05 0.1929517 0.351068574 0.36276859 -0.26267734 -0.05 0.1146316 -0.628158331 0.735049910 -0.44651083 -0.05 0.1929517 -0.902673594 0.816650411 -0.38816462 -0.05 0.2712719 -0.649270674 0.741918312 0.02013525 -0.05 0.1929517 0.159665854 0.436572213 -0.03877614 -0.05 0.1459596 0.029378254 0.488281514 0.41014395 -0.05 0.2712719 0.883469050 0.188491415 0.02782689 -0.05 0.8978331 0.082135687 0.467269416 0.11878306 -0.05 0.2712719 0.324060781 0.372946017 0.56506605 -0.05 0.2712719 1.180917005 0.118817818 0.47054422 -0.05 0.1929517 1.185040809 0.118000719 -0.05142908 -0.05 0.2712719 -0.002743819 0.501094620 0.06940159 -0.05 0.1929517 0.271822740 0.392879221 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.43159536attr(,"gstari")[1] FALSEattr(,"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.3376alternative hypothesis: two.sidedsample 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.484043812 0.8370573121 -0.05 0.4404798 1.336560700 0.090683043 -1.1222105809 -0.05 0.4404798 -1.615537694 0.946902854 0.2373930487 -0.05 0.4404798 0.433025295 0.332498205 -0.1380790378 -0.05 0.1168494 -0.257667335 0.601668176 -0.1081663716 -0.05 0.1977570 -0.130799494 0.552033047 0.1159985869 -0.05 0.1977570 0.373283232 0.354468838 0.7082098711 -0.05 0.1977570 1.704996632 0.044097539 0.0612709718 -0.05 0.1168494 0.325513260 0.3723963210 0.4204181464 -0.05 0.1977570 1.057835549 0.1450652111 0.3774028474 -0.05 0.2786646 0.809648516 0.2090711112 -0.0868802622 -0.05 0.1977570 -0.082933137 0.5330476513 0.1244368071 -0.05 0.1492124 0.451580986 0.3257854414 0.1874982314 -0.05 0.2786646 0.449903626 0.3263899715 -0.5517955762 -0.05 0.9259254 -0.521481405 0.6989842716 -0.1211737778 -0.05 0.2786646 -0.134827702 0.5536259517 -0.0030386117 -0.05 0.2786646 0.088961079 0.4645564218 -0.4168336708 -0.05 0.1977570 -0.824903760 0.7952868819 -0.0539797929 -0.05 0.2786646 -0.007539102 0.5030076420 -0.3916838150 -0.05 0.1977570 -0.768348944 0.7788600521 -0.0001822581 -0.05 0.1492124 0.128967879 0.44869153

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

 

 

更多精彩内容,扫码关注微信公众号!



【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3