揭秘Moran‘s I与Geary‘s C：如何用R语言精准识别空间自相关

第一章：揭秘Moran's I与Geary's C：空间自相关的R语言探索

在空间数据分析中，理解地理现象的空间依赖性是核心任务之一。Moran's I 与 Geary's C 是衡量空间自相关性的两个经典统计量，它们帮助我们判断邻近区域的观测值是否呈现出相似性或差异性。Moran's I 更关注全局趋势，正值表示正向空间自相关，负值则意味着负向自相关；而 Geary's C 对局部差异更为敏感，其值接近 1 表示无自相关，小于 1 表示正相关，大于 1 则为负相关。

数据准备与空间权重矩阵构建

使用 R 语言进行分析时，首先需要加载必要的包并构建空间权重矩阵。常用 `spdep` 包来定义邻接关系。

# 加载所需库
library(spdep)
library(sf)

# 读取空间数据（例如：区县多边形）
nc <- st_read(system.file("shape/nc.shp", package="sf"))

# 构建邻接列表
nb <- poly2nb(nc)

# 创建空间权重矩阵（行标准化）
listw <- nb2listw(nb, style = "W")

Moran's I 与 Geary's C 的计算

接下来可利用观测变量（如人均收入）计算两个指标：

# 假设变量为 'income'
moran_result <- moran.test(nc$BIR74, listw)
geary_result <- geary.test(nc$BIR74, listw)

print(moran_result)
print(geary_result)

• Moran's I 统计量基于协方差思想，强调相似值的空间聚集
• Geary's C 使用差值平方，对局部变化更敏感
• 两者均需通过显著性检验判断结果可靠性

统计量	范围	解释
Moran's I	约 [-1, 1]	>0: 正自相关；<0: 负自相关
Geary's C	约 [0, 2]	<1: 正自相关；>1: 负自相关

graph TD A[加载空间数据] --> B[构建邻接关系] B --> C[生成空间权重] C --> D[计算Moran's I / Geary's C] D --> E[检验显著性] E --> F[解释空间模式]

第二章：空间自相关理论基础与R实现

2.1 Moran's I指数的数学原理与空间权重构建

Moran's I 是衡量空间自相关性的核心统计量，其数学表达式为：

I = (n / S₀) * ΣᵢΣⱼ wᵢⱼ (xᵢ - x̄) (xⱼ - x̄) / Σᵢ (xᵢ - x̄)²

其中，n 为样本数量，wᵢⱼ 构成空间权重矩阵，S₀ = ΣᵢΣⱼ wᵢⱼ 为权重总和，x̄ 为变量均值。该公式通过协方差结构量化邻近区域属性值的相似性。

空间权重矩阵的构建策略

常用的空间权重定义包括邻接关系（Rook或Queen）和距离衰减函数（如反距离）。例如，基于欧氏距离的权重可表示为：

import numpy as np
def inverse_distance_weight(coords):
    n = len(coords)
    w = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            if i != j:
                d = np.linalg.norm(coords[i] - coords[j])
                w[i][j] = 1 / d if d > 0 else 0
    return w / np.sum(w)  # 标准化

该代码实现反距离权重矩阵的生成，距离越近的单元对空间相关性影响越大。

权重矩阵的标准化处理

为消除规模偏差，通常对原始权重进行行标准化，使每行权重之和为1，从而提升Moran’s I的解释一致性。

2.2 Geary's C统计量的差异性测度机制解析

Geary's C是一种用于探测空间自相关的统计量，侧重于衡量相邻区域之间的属性值差异程度。其核心思想是通过比较各空间单元与其邻居的偏差来识别聚集或离散模式。

数学表达与参数解析

Geary's C的公式定义如下：

C = \frac{(n - 1)}{2 \sum_{i,j} w_{ij}} \cdot \frac{\sum_{i,j} w_{ij} (x_i - x_j)^2}{\sum_{i} (x_i - \bar{x})^2}

其中，\(n\) 为区域总数，\(w_{ij}\) 表示空间权重矩阵元素，\(x_i\) 和 \(x_j\) 为区域属性值，\(\bar{x}\) 为均值。当 \(C < 1\) 时，呈现正空间自相关；\(C > 1\) 则为负相关。

输出对比表

C值范围	空间模式解释
< 1	相似值聚集（高-高或低-低）
≈ 1	无显著空间自相关
> 1	相异值相邻（异质性高）

2.3 空间邻接关系与空间权重矩阵的R语言建模

空间邻接关系的基本概念

在空间计量分析中，空间邻接关系用于描述地理单元之间的相邻结构。常见的定义方式包括Rook邻接（共享边界）和Queen邻接（共享边界或顶点）。这些关系将被转化为二元邻接矩阵。

构建空间权重矩阵的R实现

使用R语言中的spdep包可高效构建空间权重矩阵：

# 加载必要库
library(spdep)
library(sf)

# 读取空间多边形数据
nc <- st_read(system.file("shape/nc.shp", package="sf"))

# 构建Queen邻接关系
nb_q <- poly2nb(nc, queen = TRUE)

# 生成行标准化的空间权重矩阵
w_mat <- nb2listw(nb_q, style = "W", zero.policy = TRUE)

上述代码中，poly2nb识别相邻多边形，nb2listw将邻接列表转换为可用于回归分析的空间权重矩阵，参数style = "W"表示行标准化。

权重矩阵的可视化

可通过plot函数结合lines展示邻接结构。

2.4 全局与局部空间自相关的概念辨析与可视化

空间自相关衡量地理空间中邻近位置属性值的相似性程度。全局空间自相关（如Moran's I）评估整个研究区域内是否存在聚集模式，而局部空间自相关（如LISA）则识别特定位置周围的聚类或异常。

全局与局部指标对比

• 全局Moran's I：反映整体空间依赖性，值域[-1,1]，接近1表示强正相关。
• LISA：分解全局指标至每个空间单元，可识别热点（高-高）、冷点（低-低）及异常点（高-低、低-高）。

可视化代码示例

from esda.moran import Moran, Moran_Local
import matplotlib.pyplot as plt

# 计算全局Moran's I
moran = Moran(y=values, w=w_matrix)
print(f"Global Moran's I: {moran.I:.3f}, p-value: {moran.p_sim:.4f}")

# 计算局部LISA
lisa = Moran_Local(y=values, w=w_matrix)
plt.scatter(values, lisa.Is); plt.xlabel('Value'); plt.ylabel('Local Moran\'s I')

上述代码首先计算全局自相关以判断整体聚集趋势，随后通过局部指标生成散点图，揭示各区域与其邻居的关系类型，为后续空间聚类分析提供基础。

2.5 假设检验与显著性评估在R中的实现方法

基本假设检验流程

在R中执行假设检验通常包括设定原假设、选择合适统计检验方法、计算p值并做出决策。常用函数如t.test()、chisq.test()和prop.test()可快速完成分析。

独立样本t检验示例

# 生成两组样本数据
group_a <- rnorm(30, mean = 50, sd = 10)
group_b <- rnorm(30, mean = 55, sd = 10)

# 执行双样本t检验
result <- t.test(group_a, group_b)
print(result)

该代码比较两组正态分布样本的均值差异。t.test()默认采用Welch校正，不假定方差齐性；输出包含t统计量、自由度、p值及置信区间，用于判断是否拒绝原假设。

常见检验方法对照表

数据类型	检验目的	R函数
连续变量	均值比较	t.test()
分类变量	独立性检验	chisq.test()
比例数据	比例差异	prop.test()

第三章：R语言空间数据分析环境搭建

3.1 sp、sf与spdep包的核心功能与数据结构

空间数据的表示与操作

R语言中的sp包是早期处理空间数据的基础工具，提供SpatialPoints、SpatialPolygons等类来封装地理信息。其核心在于将几何对象与属性数据绑定，支持基本的空间查询与投影变换。

现代空间数据框架：sf包

sf包采用简单要素（Simple Features）标准，使用sfc和sfg类型统一表达点、线、面。数据以data.frame扩展形式存储，兼容dplyr操作。

library(sf)
nc <- st_read("shapefile.shp")
st_crs(nc) # 查看坐标系

上述代码加载Shapefile并查看其CRS，st_read()自动解析几何与属性，st_crs()提取空间参考系统。

空间依赖建模：spdep包

spdep包构建空间权重矩阵，用于莫兰指数等分析。通过邻接关系生成listw对象，支撑空间回归模型构建。

3.2 空间数据的读取、转换与拓扑关系处理

在地理信息系统（GIS）开发中，空间数据的读取是首要步骤。常用工具如GDAL/OGR支持多种格式的矢量与栅格数据读取。以Python为例，可通过以下代码加载Shapefile：

from osgeo import ogr
# 打开数据源
dataSource = ogr.Open("roads.shp")
layer = dataSource.GetLayer()
for feature in layer:
    geom = feature.GetGeometryRef()
    print(geom.ExportToWkt())  # 输出几何对象的WKT表示

该代码段打开Shapefile文件并遍历图层中的每个要素，获取其几何信息并以WKT格式输出，便于后续解析与分析。

坐标系转换

空间数据常需统一坐标参考系统（CRS）。使用OSR模块可实现投影转换：

from osgeo import osr
# 定义源和目标坐标系
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(4326)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(3857)
# 创建转换器
transform = osr.CoordinateTransformation(src_srs, dst_srs)

此转换器可用于点坐标的动态重投影。

拓扑关系处理

基于Shapely库可判断要素间的拓扑关系，如包含、相交等：

• geom1.intersects(geom2)：判断是否相交
• geom1.contains(geom2)：判断包含关系
• geom1.touches(geom2)：判断是否接触

这些谓词为地理分析提供逻辑基础。

3.3 构建可复用的空间权重对象：listw详解

在空间计量分析中，`listw` 对象是构建空间依赖关系的核心数据结构。它封装了空间单元间的邻接关系与权重信息，支持多种空间矩阵形式。

listw 的生成流程

通常由邻接列表（`nb`）和权重类型共同构建：

library(spdep)
# 假设已构建 nb 邻接列表
listw <- nb2listw(nb, style = "W", zero.policy = TRUE)

其中，`style = "W"` 表示行标准化，使每行权重和为1；`zero.policy = TRUE` 允许存在孤立区域。

常用权重风格对比

风格	说明
W	行标准化，最常用
B	二值化权重
C	列标准化

重用 `listw` 可避免重复计算，提升模型效率，适用于空间滞后、莫兰指数等分析场景。

第四章：空间自相关建模实战案例分析

4.1 城市GDP分布的全局Moran's I检验与解读

在空间计量分析中，全局Moran's I是衡量区域间属性值空间自相关性的核心指标。通过对城市GDP数据进行检验，可判断高GDP或低GDP城市是否在地理上呈现集聚分布。

计算流程与代码实现

from esda.moran import Moran
import numpy as np

# 假设gdp_data为各城市的GDP数组，w为标准化的空间权重矩阵
moran = Moran(gdp_data, w)
print(f"Moran's I: {moran.I:.3f}")
print(f"P-value: {moran.p_sim:.4f}")

上述代码利用`esda`库中的`Moran`类计算全局Moran's I。其中，`gdp_data`表示城市GDP观测值，`w`为空间邻接权重矩阵（通常经行标准化处理）。输出的`I`值若显著大于0，表明存在正向空间自相关，即相似经济水平的城市趋于集聚。

结果解读

• 当Moran's I > 0：表示高-高或低-低集聚；
• I ≈ 0：无明显空间模式；
• I < 0：呈现分散或棋盘状分布。

结合p值可判断统计显著性，通常p < 0.05认为空间自相关显著。

4.2 居民收入差异的Geary's C空间模式识别

空间自相关分析原理

Geary's C是一种用于探测空间自相关的统计量，值域通常在0到2之间。C接近1表示无显著空间聚集性，小于1表明存在正向空间自相关（相似值聚集），大于1则为负向自相关。

计算实现与代码解析

from pysal.explore import esda
from pysal.lib import weights

# 构建空间权重矩阵
w = weights.Queen.from_dataframe(income_data)
w.transform = 'r'

# 计算Geary's C
geary = esda.Geary(income_data['income'], w)
print(f"Geary's C: {geary.C}, p-value: {geary.p_val}")

上述代码首先基于邻接关系构建Queen权重矩阵，随后标准化并计算Geary's C统计量。参数transform='r'表示行标准化，确保权重可比性。

结果解读

Geary's C	p-value	解释
0.78	0.002	显著正向空间聚集

低C值配合显著p值，说明居民收入呈现明显空间集聚特征。

4.3 LISA聚类图绘制与热点区域探测

LISA（Local Indicators of Spatial Association）聚类图是探测空间数据中局部聚集模式的重要工具，广泛应用于地理信息系统与区域经济分析中。

核心计算流程

使用PySAL库实现LISA指标的计算：

from pysal.explore import esda
from pysal.lib import weights

# 构建空间权重矩阵
w = weights.Queen.from_dataframe(gdf)
w.transform = 'r'

# 计算局部莫兰指数
lisa = esda.Moran_Local(gdf['value'], w)

上述代码首先基于邻接关系构建Queen权重矩阵，并进行行标准化。Moran_Local函数逐区域计算局部空间自相关性，输出每个区域的LISA值及其显著性。

显著性分类与可视化

根据LISA结果可划分四类区域：高-高（HH）、低-低（LL）、高-低（HL）、低-高（LH）。通过如下逻辑判定热点区域：

• HH：高值被高值包围，显著正相关，视为热点区
• LL：低值被低值包围，显著负相关，视为冷点区
• HL/LH：空间异质性强，为异常区域

4.4 空间自相关结果的稳健性检验与报告生成

稳健性检验策略

为确保空间自相关分析结果的可靠性，需进行多维度稳健性检验。常用方法包括：调整空间权重矩阵的构建方式（如邻接、距离衰减）、更换全局或局部指标（如Moran’s I与Getis-Ord G*对比），以及进行置换检验（permutation test）以评估显著性。

1. 使用不同空间尺度重新计算指标
2. 执行999次随机置换获取p值分布
3. 比较子区域与全区域结果一致性

自动化报告生成示例

通过Python脚本集成检验流程并输出结构化报告：

from pysal.explore import esda
from pysal.lib import weights

# 构建空间权重矩阵
w = weights.Queen.from_dataframe(geo_data)
w.transform = 'r'

# 计算Moran's I并进行置换检验
moran = esda.Moran(y, w, permutations=999)
print(f"Moran's I: {moran.I:.3f}, p-value: {moran.p_sim:.3f}")

上述代码中，permutations=999启用蒙特卡洛模拟评估统计显著性，w.transform='r'对权重矩阵行标准化，确保结果可比性。

第五章：结论与空间计量分析的进阶方向

空间计量模型不仅揭示了地理单元间的相互依赖性，还为区域经济、城市规划和环境科学提供了强有力的分析工具。随着数据粒度提升与计算能力增强，其应用正从静态横截面分析转向动态面板与高维空间网络建模。

非线性空间杜宾模型的实现

在处理复杂空间溢出效应时，传统线性假设常显不足。以下为基于 R 的非线性空间杜宾模型拟合示例：

library(spdplyr)
library(spdep)

# 构建空间权重矩阵
w <- nb2listw(knn2nb(knearneigh(coordinates(data), k=5)))

# 拟合非线性SDM（含平方项与交互项）
model_nlsdm <- lagsarlm(
  log(GDP) ~ POP + INVEST + EDUC + I(POP^2) + POP:INVEST,
  data = data, 
  listw = w, 
  type = "mixed"
)
summary(model_nlsdm)

多尺度地理加权回归（MGWR）的应用场景

相较于传统GWR，MGWR允许不同变量在不同带宽下进行局部估计，更符合现实机制。例如，在分析中国城市群创新效率时，研发资本的邻近影响范围可达300公里，而人才流动仅限于100公里圈层。

• 使用Python的mGWR库支持并行带宽选择
• 通过AICc准则比较MGWR与GWR拟合优度
• 可视化各变量带宽分布以识别空间作用尺度差异

空间机器学习融合路径

将图神经网络（GNN）与空间滞后变量结合，可捕捉非欧几里得空间依赖。如在房价预测中，以行政区为节点构建图结构，嵌入POI密度与地铁连通性作为边权重：

特征类型	数据来源	空间权重角色
商业设施密度	高德API	节点属性
地铁直达频次	交通卡口数据	边权重
学区评级	教育局公开数据	节点属性