计算全局莫兰指数(Moran's I)

简介

Moran's I是澳大利亚统计学家莫兰在1950年提出的一种测量空间自相关的测量方法。空间自相关指信号在空间的邻近位置之间呈现相关性，因为是多维和多方向的，它比一维自相关更复杂。莫兰全局指数可以反映数据离散或者聚集的程度。例如，这里白方和黑方是完全分散的，此时全局莫兰指数将是-1。

计算

计算公式如下(Wiki)：

其中x为数据，w为空间权重，W为所有空间权重的和。

权重矩阵w的计算使用较多的是距离矩阵，其中对角线为0，单元格元素为距离的倒数。计算距离的方式有多种比如欧式距离，测地距离等，此外，这里的距离的概念不限于空间，还可以是广义的距离，比如人和人之间的心理距离，区域经济学中的经济距离，世界上最遥远的距离等。

如果将这个简化过莫兰指数公式转换一下会更有助于理解其计算的逻辑：将分子的N移到分母成为1/N，此时分母的公式正是计算方差的公式，而蓝色部分可以看作是用空间权重加权计算的方差，全局莫兰指数反映的便是空间加权方差和数据方差的比。

实际计算可以根据简化的公式进行，从公式看稍微复杂的地方便是分子公式的计算，向量化的计算就是对向量求外积。

s <- sum(weight)
m <- mean(x) #计算x平均
y <- x - m #计算离差
wv <- sum(weight * y %o% y) #对y做外积再加入空间权重信息
v <- sum(y^2) #离差平方和
obs <- (n/s) * (wv/v) #计算莫兰指数

实际计算中还要提前对权重矩阵做了一个标准化。

ROWSUM <- rowSums(weight) 
ROWSUM[ROWSUM == 0] <- 1
weight <- weight/ROWSUM

对于莫兰指数可以计算z分数获取显著性，公式如下：

其中E[I]是期望，V[I]是方差。

V[I]公式中复杂的地方就是前面部分的计算，

带入I的计算公式，化简后如下

根据公式写代码就可以计算出E[I]和SD[I]：

ei <- -1/(n - 1)
S1 <- 0.5 * sum((weight + t(weight))^2)
S2 <- sum((apply(weight, 1, sum) + apply(weight, 2, sum))^2)
s.sq <- s^2
k <- (sum(y^4)/n)/(v/n)^2
vi <- (n * ((n^2 - 3 * n + 3) * S1 - n * S2 + 3 * s.sq) - k * (n * (n - 1) * S1 - 2 * n * S2 + 6 * s.sq))/((n -  1) * (n - 2) * (n - 3) * s.sq) - 1/((n - 1)^2)
sdi <- sqrt(vi)

知道平均数和标准差，就可以根据观测到的莫兰指数(obs)在对应的正态分布中的位置确定其显著性。

pv <- pnorm(obs, mean = ei, sd = sdi) # less
pv <- 1 - pv # greater
pv <- ifelse (obs <= ei, 2 * pv, 2 * (1 - pv)) # two.sided

解释

空间自相关（Global Moran's I）其实是一种统计推断工具，这意味着分析结果需要在无效假设的范围内进行解释。

对于Global Moran's I统计，无效假设是：特征在研究区域中是随机分布的。当返回的p值显著时，可以拒绝无效假设。下表总结了对结果的解释。

软件实现

随便搜索一下，就发现三个软件可以实现全局莫兰指数的计算：

• Moran.I(x,w){ape}
  

• moran(x, listw,n,S0) {spdep}

• Moran {raster}
  

本文中的代码参考了ape的Moran.I(x,w)，这里其他两个函数都没有做参数检验，只返回了全局莫兰指数。但spdep工具包提供了moran.test有更完善的功能，它还检测了参数检验的假设--正态性，另外spdep还提供了非参数检验的方法moran.mc。spdep中listw是w矩阵的另外一种稀疏表达形式，该软件的计算大部分都基于listw，可用mat2listw构建。如果你想要用其他软件(这里说的是没有做标准化的算法，比如spdep)得到和ape相同的结果，需要对w矩阵做标准化(!)。

ROWSUM <- rowSums(weight) 
ROWSUM[ROWSUM == 0] <- 1
weight <- weight/ROWSUM

至于到底需不需要做标准化，一个很简单的验证方式就是生成随机数据，由于莫兰指数的期望只和n相关，可以看随机数据算出的平均莫兰指数是否和期望相匹配。比如，随机生成1000份Desikan模板34个区域的数据，对w做标准化后，平均莫兰指数为-0.03，如果不做标准化，平均莫兰指数为-0.09。前者更符合莫兰指数的期望值-0.03，因此我觉得是需要做标准化的。由于莫兰指数是一个相对的概念，做不做标准化对z分数和p值应该影响不大(?存疑)。

参考

Paradis E, Schliep K (2019). “ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R.” Bioinformatics, 35, 526-528.

https://pro.arcgis.com/en/pro-app/2.8/tool-reference/spatial-statistics/h-how-spatial-autocorrelation-moran-s-i-spatial-st.htm

CSDN: 空间相关分析(一/二/三) 

2022年8月