r语言自相关图和偏自相关图怎么看 r语言自相关检验代码

空间自相关指数又称莫兰指数，是空间分析常采用的指标，但是使用不同软件计算出的莫兰指数有时会不一致，这是因为不同软件设定的默认选项不一样。本篇介绍如何在R语言中计算莫兰指数和局部莫兰指数，使用的工具包为spdep。

该包名称是“Spatial Dependence”的缩写，是R语言中专门做空间相关性分析的工具包。在spdep中，计算莫兰指数的过程分为三个步骤，即根据矢量对象创建空间邻接矩阵、根据邻接矩阵创建空间权重矩阵和根据权重矩阵计算空间自相关指数。

前面已经介绍，在R中空间矢量对象有sp和sf两种不同的存储格式，spdep对两种格式都适用，且使用方法基本一致。这里以sf格式的nc数据集为例，并使用as函数将其转换为sp格式作为备用，以在必要的情况下说明两种格式在使用方法上的某些细微差别。

library(sf)
library(sp)
nc <- st_read(system.file("shape/nc.shp", package = "sf"))
nc.sp <- as(nc, "Spatial")

1 创建邻接矩阵

spdep工具包可以分别以面要素和点要素为对象建立空间邻接矩阵，前者主要以邻接关系为依据，后者主要以距离为依据。

1.1 基于面要素的空间邻接矩阵

一阶邻接矩阵

判断两个面要素是否具有邻接关系有两种规则：一是两者之间必须具有共同的边界，即ROOK规则，另一种是两者至少有一个交点，即QUEEN规则。相应的函数及其语法规则为：

poly2nb(pl, row.names = NULL, snap=sqrt(.Machine$double.eps),
        queen=TRUE, useC=TRUE, foundInBox=NULL)

• pl：面要素对象；
• queen：是否选择QUEEN规则判断邻接关系，默认为TRUE，FALSE表示使用ROOK规则。

使用plot函数绘制空间邻接矩阵的语法规则：

plot(x, coords, col="black", points=TRUE, add=FALSE,
     arrows=FALSE, length=0.1, xlim=NULL, ylim=NULL, ...)

• x：空间邻接矩阵；
• coords：绘图的坐标信息；sp对象或它的坐标；sf对象的sfc信息；
• points：在绘制时是否用点代替面要素的位置，默认为TRUE；
• add：图层叠加参数；
• arrows：在邻接关系非对称时是否使用箭头表示连接方向，默认为FALSE；
• 该函数的输出结果实际上并不是矩阵，而是与x等长度的list；每个元素存储与对应面要素存在邻接关系的要素编码。

library(spdep)

nb <- poly2nb(nc)
nb2 <- poly2nb(nc, queen = F)

# 绘制邻接关系图
plot(st_geometry(nc), main = "QUEEN")
plot(nb, st_geometry(nc), add = T, col = "red")

plot(st_geometry(nc), main = "ROOK")
plot(nb2, st_geometry(nc), arrows = T, add = T, col = "blue")

# sp对象
nb.sp <- poly2nb(nc.sp)
plot(nc.sp, main = "sp")
plot(nb.sp, nc.sp, add = T, col = "red")

高阶邻接矩阵

poly2nb函数生成的是一阶邻接矩阵，在此基础上使用nblag和nblag_cumul函数可以生成更高阶的邻接矩阵。语法结构如下：

nblag(neighbours, maxlag)
nblag_cumul(nblags)

• neighbours：poly2nb函数输出的邻接矩阵；
• maxlag：邻接阶数；
• nblag用于输出邻接矩阵的滞后值，nblag_cumul用于将其转换为新的空间邻接矩阵。

nb.lag <- nblag(nb, maxlag = 3)
nb.lag2 <- nblag_cumul(nb.lag)

class(nb.lag)

## [1] "list"

class(nb.lag2)

## [1] "nb"

1.2 基于点要素的空间邻接矩阵

提取nc的质心作为点要素对象：

ncpoint <- st_centroid(st_geometry(nc))

# sp对象
nc.coord <- coordinates(nc.sp)

k邻接方法

k邻接方法就是基于点之间距离的排序，将距离点最近的k个点作为该点的邻接要素。使用的函数为knearneigh和knn2nb函数：

knearneigh(x, k = 1, longlat = NULL, use_kd_tree=TRUE)
knn2nb(knn, row.names = NULL, sym = FALSE)

• x：输入点要素；
• k：k邻接方法的k参数；
• knn：knearneigh的输出对象；
• knearneigh函数用于储存邻接对象，knn2nb函数用于将其转换为空间邻接矩阵。

nb.k <- knn2nb(knearneigh(ncpoint))
nb.k2 <- knn2nb(knearneigh(ncpoint, k = 2))
nb.k4 <- knn2nb(knearneigh(ncpoint, k =4))

plot(st_geometry(nc), main = "k = 1")
plot(nb.k, ncpoint, arrows = T, add = T, col = "red")

plot(st_geometry(nc), main = "k = 2")
plot(nb.k2, ncpoint, add = T, col = "red")

plot(st_geometry(nc), main = "k = 4")
plot(nb.k4, ncpoint, add = T, col = "red")

使用nbdists函数可以输出具有邻接关系的要素之间具体的距离：

nbdists(nb, coords, longlat = NULL)

• nb：空间邻接矩阵对象；
• coords：点坐标信息，可以是点要素，也可以是距离矩阵。

nbdists(nb.k, ncpoint)
nbdists(nb.k2, ncpoint)
nbdists(nb.k4, ncpoint)

基于距离范围

该方法是将与点要素的距离在一定范围内的要素作为其邻接要素。使用函数为dnearneigh：

dnearneigh(x, d1, d2, row.names = NULL, longlat = NULL, bounds=c("GT", "LE"))

• d1和d2分别为距离范围的下限和上限。

max = max(unlist(nbdists(nb.k, ncpoint)))
nb.dist <- dnearneigh(ncpoint, 0, 0.5*max)
nb.dist2 <- dnearneigh(ncpoint, 0, max)
nb.dist3 <- dnearneigh(ncpoint, 0, 1.5*max)

plot(st_geometry(nc), main = "小范围")
plot(nb.dist, ncpoint, add = T, col = "red")

plot(st_geometry(nc), main = "中范围")
plot(nb.dist2, ncpoint, add = T, col = "red")

plot(st_geometry(nc), main = "大范围")
plot(nb.dist3, ncpoint, add = T, col = "red")

基于图形关系

除了根据距离判断点要素之间的邻接关系外，还可以根据点要素分布的图形关系建立邻接关系，优点是邻接要素之间的连线不会出现交叉，更符合人们对“邻接”概念的认知。

根据Delaunay三角剖分算法建立邻接矩阵：

tri2nb(coords, row.names = NULL)

nb.tri <- tri2nb(ncpoint)

plot(st_geometry(nc), main = "Delaunay triangulation")
plot(nb.tri, ncpoint, add = T, col = "red")

除Delaunay算法外，还有以下函数可以建立点要素的图形关系：

• gabrielneigh函数；
• relativeneigh函数；
• soi.graph函数：输入对象为tri2nb函数的输出对象，使用前需加载dbscan工具包。

以上函数只是建立图形关系，输出对象还需要使用graph2nb函数才能生成空间邻接矩阵。

gabrielneigh(coords, nnmult=3)
relativeneigh(coords, nnmult=3)
soi.graph(tri.nb, coords, quadsegs=10)
graph2nb(gob, row.names=NULL,sym=FALSE)

• 具体算法可查看相关资料。

graph <-  gabrielneigh(ncpoint)
nb.graph <- graph2nb(graph)
plot(st_geometry(nc), main = "Gabriel neighbor graph")
plot(nb.graph, ncpoint, add = T, col = "red")

graph2 <- relativeneigh(ncpoint)
nb.graph2 <- graph2nb(graph2)
plot(st_geometry(nc), main = "Relative neighbor graph")
plot(nb.graph2, ncpoint, add = T, col = "red")

library(dbscan)
graph3 <- soi.graph(nb.tri, ncpoint)
nb.graph3 <- graph2nb(graph3)
plot(st_geometry(nc), main = "SOI neighbor graph")
plot(nb.graph3, ncpoint, add = T, col = "red")

2 创建空间权重矩阵

2.1 二分式的空间权重矩阵

二分式的空间权重矩阵即对与某要素邻接和不邻接的要素分别赋予不同的值，spdep工具包中的赋值方法有以下几种：

• B型：具有邻接关系的赋值为1，否则赋值为0；
• W型：对B型空间权重矩阵进行行标准化，即使每行元素之和都为1；
• U型：对B型空间权重矩阵进行全局标准化，即使所有元素之和为1；
• C型：U型空间权重矩阵乘以矢量对象所包含的空间单元个数；
• minmax型：对B型空间权重矩阵分别求出最大行和和最大列和，然后让B型空间权重矩阵除以其中的较小者；
• S型：先对B型空间权重矩阵的每个元素进行开方得到矩阵Q，然后对Q进行类似W型的行标准化得到P，再对P进行类似C型的操作。

• B、U、C和minmax型空间权重矩阵相互之间具有倍数关系。

生成空间权重矩阵的函数及其语法结构如下：

nb2listw(neighbours, glist=NULL, style="W", zero.policy=NULL)
nb2mat(neighbours, glist=NULL, style="W", zero.policy=NULL)
listw2mat(listw)

• neighbours：前文任意一种类型的空间邻接矩阵；
• style：权重赋值方法，默认为W型；ArcGIS中默认为B型，GeoDa中默认也为W型；
• nb2listw函数输出的空间权重矩阵的数据类型为list，nb2mat函数输出的为矩阵；listw2mat函数可以将前者转换为后者。

w1 <- nb2listw(nb)
w2 <- nb2mat(nb, style = "B")
w3 <- listw2mat(w1)

class(w1)

## [1] "listw" "nb"

class(w2)

## [1] "matrix"

class(w3)

## [1] "matrix"

2.2 距离衰减式的空间权重矩阵

该方法对邻接要素之间的权重按照距离衰减的原则进行赋值，使用的函数及其语法结构如下：

nb2listwdist(neighbours, x, type="idw", style="raw", 
             alpha = 1, dmax = NULL, longlat = NULL,
             zero.policy=NULL)

• neighbours：前文任意一种类型的空间邻接矩阵；
• x：空间矢量对象，可以为sp、sf和sfc对象；
• type：距离衰减函数，默认为idw（反距离衰减），可选的有exp（指数衰减）和dpd；
• style：空间权重矩阵的标准化类型，默认值raw表示不标准化，可选的同nb2listw的style参数，即W、B、C、U、minmax和S；
• alpha、dmax：分别相当于距离衰减函数中的 和 参数。

距离衰减函数的表达式：

idw：

exp：

dpd：

nb2listwdist(nb, nc)
nb2listwdist(nb.k2, nc)

3 计算空间自相关指数

3.1 全局莫兰指数

全局莫兰指数的表达式为：

数学期望只与空间单元数量有关：

上式中， 为空间单元数量， 为空间权重矩阵的元素， 为属性值， 为空间权重矩阵的元素之和。

计算莫兰指数的函数为moran，语法结构如下：

moran(x, listw, n, S0, zero.policy=NULL, NAOK=FALSE)

• x、n和S0参数的含义同表达式；
• spdep中的Szero函数专门用于计算 ；
• listw：空间权重矩阵，必须为list数据类型。

moran.plot函数用于绘制莫兰散点图：

moran.plot(x, listw, zero.policy=NULL, spChk=NULL,
           labels=NULL, xlab=NULL, ylab=NULL, quiet=NULL, 
           plot=TRUE, return_df=TRUE, ...)

moran.test函数用于参数检验，但在计算莫兰指数时实际上比moran函数更方便：

moran.test(x, listw, randomisation=TRUE, zero.policy=NULL,
           alternative="greater", rank = FALSE, na.action=na.fail, 
           spChk=NULL, adjust.n=TRUE, drop.EI2=FALSE)

• alternative：假设检验的选项，greater和less为单边减压，two.sided为双边检验；默认值为greater。

moran.mc函数使用蒙特卡洛试验进行显著性检验：

moran.mc(x, listw, nsim, zero.policy=NULL, 
         alternative="greater", na.action=na.fail,
         spChk=NULL, return_boot=FALSE, adjust.n=TRUE)

# 计算空间权重矩阵
weight.nc <- nb2listw(nb)

moran.plot(nc$PERIMETER, weight.nc)
a <- moran.plot(nc$PERIMETER, weight.nc)
head(a)
##       x       wx is_inf labels      dfb.1_         dfb.x       dffit    cov.r
## 1 1.442 1.501000  FALSE      1 -0.04278069  0.0278853957 -0.06430210 1.026424
## 2 1.231 1.685333  FALSE      2  0.04605112 -0.0365300971  0.05392410 1.036641
## 3 1.630 1.478600  FALSE      3 -0.03264263  0.0080770921 -0.09058791 1.014197
## 4 2.968 2.593500   TRUE      4 -0.43006790  0.5006278266  0.53389501 1.043857
## 5 2.206 1.861500  FALSE      5  0.01243412 -0.0174393665 -0.02346575 1.043528
## 6 1.670 1.880333  FALSE      6  0.01661839 -0.0003639193  0.05900216 1.023748
##         cook.d        hat
## 1 0.0020815773 0.01231623
## 2 0.0014665581 0.01848151
## 3 0.0041112280 0.01008014
## 4 0.1394532266 0.08282389
## 5 0.0002780899 0.02233750
## 6 0.0017523463 0.01000038

moran(nc$PERIMETER, weight.nc, n = nrow(nc), S0 = Szero(weight.nc))

## $I
## [1] 0.3275279
## 
## $K
## [1] 5.947649

• I为莫兰指数，K为样本峰度。

moran.test(nc$PERIMETER, weight.nc)

## 
##  Moran I test under randomisation
## 
## data:  nc$PERIMETER  
## weights: weight.nc    
## 
## Moran I statistic standard deviate = 5.2608, p-value = 7.17e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.327527853      -0.010101010       0.004118771

moran.mc(nc$PERIMETER, weight.nc, nsim = 1000)

# 
##  Monte-Carlo simulation of Moran I
## 
## data:  nc$PERIMETER 
## weights: weight.nc  
## number of simulations + 1: 1001 
## 
## statistic = 0.32753, observed rank = 1001, p-value = 0.000999
## alternative hypothesis: greater

比较使用不同类型的空间权重矩阵的计算结果：

moran.test(nc$PERIMETER, nb2listw(nb, style = "W"))$estimate[1]

## Moran I statistic 
##         0.3275279

moran.test(nc$PERIMETER, nb2listw(nb, style = "B"))$estimate[1]

## Moran I statistic 
##         0.2333919

moran.test(nc$PERIMETER, nb2listw(nb, style = "C"))$estimate[1]

## Moran I statistic 
##         0.2333919

moran.test(nc$PERIMETER, nb2listw(nb, style = "U"))$estimate[1]

## Moran I statistic 
##         0.2333919

moran.test(nc$PERIMETER, nb2listw(nb, style = "minmax"))$estimate[1]

## Moran I statistic 
##         0.2333919

moran.test(nc$PERIMETER, nb2listw(nb, style = "S"))$estimate[1]

## Moran I statistic 
##          0.274386

• 由于B、C、U和minmax型空间权重矩阵相互之间具有倍数关系，因此计算出的莫兰指数相同。

3.2 局部莫兰指数

局部莫兰指数，也称LISA指数，表达式为：

计算局部莫兰指数的函数为localmoran，语法结构如下：

localmoran(x, listw, zero.policy=NULL, na.action=na.fail, 
           alternative = "greater", p.adjust.method="none", mlvar=TRUE,
           spChk=NULL, adjust.x=FALSE)

localmoran_perm通过置换试验对其进行显著性检验：

localmoran_perm(x, listw, nsim=499, zero.policy=NULL, na.action=na.fail, 
                alternative = "greater", p.adjust.method="none", mlvar=TRUE,
                spChk=NULL, adjust.x=FALSE, sample_Ei=TRUE, iseed=NULL)

locaI = localmoran(nc$PERIMETER, weight.nc)
head(locaI)

##             Ii        E.Ii    Var.Ii        Z.Ii    Pr(z > 0)
## 1  0.172453378 -0.01010101 0.3105185  0.32760351 3.716057e-01
## 2 -0.023745273 -0.01010101 0.3105185 -0.02448535 5.097673e-01
## 3  0.036255898 -0.01010101 0.1826379  0.10847235 4.568105e-01
## 4  5.176465758 -0.01010101 0.4703693  7.56241468 1.978273e-14
## 5  0.436385992 -0.01010101 0.2305931  0.92979172 1.762395e-01
## 6 -0.002665337 -0.01010101 0.3105185  0.01334371 4.946768e-01

4 与其他软件的比较

ArcGIS和GeoDa平台也可以计算空间自相关指数，但是二者在创建空间权重矩阵时的默认选项是不同的。

4.1 ArcGIS的设置界面

在ArcGIS中，按照ArcToolbox -> Spatial Statistics Tools -> Analyzing Patterns -> Spatial Autocorrelation (Morans I)的操作顺序打开计算莫兰指数的窗口：

• conceptualization of Spatial Relationships：空间邻接矩阵的计算方法，默认选项为INVERSE_DISTANCE，其中CONTIGUITY_EDGES_ONLY选项相当于ROOK，CONTIGUITY_EDGES_CORNERS选项相当于QUEEN；
• Distance Method：空间邻接矩阵与距离有关时，选择计算距离的方法；
• Standardization：空间权重矩阵的计算方法，默认选项为NONE相当于B型，选项ROW相当于W型。

4.2 GeoDa的设置界面

基于面要素的邻接关系的设置界面：

基于点要素的距离关系的设置界面：

有兴趣的读者可以在两个平台上分别操作，然后观察计算结果是否存在区别。

5 手动计算

空间自相关指数并不涉及太复杂的数学知识，完全可以通过自定义程序进行手动计算。

5.1 手动计算空间权重矩阵

B型空间权重矩阵实际上就是空间邻接矩阵，其他类型的空间权重矩阵的计算都是以此为基础的。在这里，将其记为w0。

w0 <- nb2mat(nb, style = "B")

• 这里使用输出数据结构为矩阵的空间权重计算函数。

# W型
ww <- nb2mat(nb, style = "W") # 函数计算
ww2 <- w0/rowSums(w0) # 手动计算

# U型
wu <- nb2mat(nb, style = "U")
wu2 <- w0/sum(w0)

# C型
wc <- nb2mat(nb, style = "C")
wc2 <- ncol(w0)*w0/sum(w0)

# minmax型
wmin <- nb2mat(nb, style = "minmax")
wmin2 <- w0/min(max(colSums(w0)), max(rowSums(w0)))

# S型
ws <- nb2mat(nb, style = "S")
Q <- w0/sqrt(rowSums(w0))
ws2 <- ncol(w0)*Q/sum(Q)

• 手动计算和函数计算结果完全一致。

5.2 手动计算莫兰指数

按照前面莫兰指数的表达式，需要写一个包含两层循环的程序进行计算，虽然理论上可行，但是实际上可能由于循环过程中会不断进行四舍五入，使得计算结果与真实值有所偏差。为了避免循环，可以引入矩阵运算：

先计算 中心化后的数值，记为 ：

然后莫兰指数可以转换成以下形式：

手动计算过程如下：

# 加载工具包和示例数据
library(sf)
library(spdep)
nc <- st_read(system.file("shape/nc.shp", package = "sf"))

nb <- poly2nb(nc) # 计算空间邻接矩阵
X <- nc$PERIMETER # 属性变量
n <- length(X) # 空间单元个数

# 函数计算
w <- nb2listw(nb, style = "W") # 必须是list类型的权重矩阵
moran.test(X, w) # 可以使用未中心化的变量

# 手动计算
w2 <- w0/sum(w0) # 必须是矩阵类型的权重矩阵
Z <- scale(X, scale = F) # 必须对变量进行中心化
ncol(w2)*cov(Z, w2%*%Z)/(var(Z)*sum(w2))

# 部分输出结果

moran.test(X, w) # 可以使用未标准化的变量
## 
##  Moran I test under randomisation
## 
## data:  X  
## weights: w    
## 
## Moran I statistic standard deviate = 3.9951, p-value = 3.234e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.233391907      -0.010101010       0.003714697

ncol(w2)*cov(Z, w2%*%Z)/(var(Z)*sum(w2))
##           [,1]
## [1,] 0.2333919