局部莫兰指数(LISA)原理、计算与可视化:空间热点探测全解析

1. 项目概述:从全局到局部,洞察空间异质性的关键一步

在空间数据分析领域,我们常常需要回答一个核心问题:某个区域的现象,是否与其周边区域存在某种关联?全局莫兰指数(Global Moran‘s I)为我们提供了一个宏观的答案,它像一把尺子,衡量了整个研究区域内空间现象的总体聚集或离散趋势。然而,这把“全局尺”有一个明显的局限:它无法告诉我们,具体是哪些地方在“抱团取暖”(高值聚集),哪些地方在“鹤立鸡群”(高值被低值包围),或者哪些地方表现得“格格不入”(空间异常值)。这就好比看一张全国经济地图,全局指数告诉你“经济存在空间正相关”,但你不知道是“长三角、珠三角在高速聚集”,还是“某些资源型城市形成了低增长洼地”。

局部莫兰指数(Local Moran‘s I),正是为了解决这个问题而生的利器。它是对全局莫兰指数的自然深化与解构,其核心思想是为研究区域内的每一个空间单元(如区县、街道、网格)都计算一个独立的指数值。这个值,结合其统计显著性检验,能够精准地揭示每个单元与其邻近单元之间的局部空间关联模式。它不再满足于一个笼统的“整体相关”结论,而是致力于绘制一幅精细的“空间关联显微图”,将隐藏在全局趋势下的局部热点、冷点、异常点一一标识出来。

理解与计算局部莫兰指数,对于地理学、城市规划、流行病学、环境科学、房地产经济学等任何涉及空间数据的领域都至关重要。它能帮助研究者识别犯罪高发区集群、追踪传染病的扩散源头、发现环境污染的热点区域,或是评估政策实施的局部效应差异。对于数据分析师和研究者而言,掌握这项技术,意味着能从“知道有模式”进阶到“看清模式在哪里”,从而支撑起更具针对性的决策与分析。

2. 局部莫兰指数的核心原理与类型解析

要真正用好局部莫兰指数,不能只停留在调用软件函数,必须深入理解其数学本质和结果解读的逻辑。这不仅是正确计算的前提,更是避免误读误用的关键。

2.1 数学公式的直观理解

局部莫兰指数 ( I_i ) 对于空间单元 ( i ) 的计算公式为:

[ I_i = \frac{(x_i - \bar{x})}{S^2} \sum_{j=1, j\neq i}^{n} w_{ij}(x_j - \bar{x}) ]

这个公式看起来有些复杂,但我们可以将其拆解为三个核心部分来理解:

  1. ( (x_i - \bar{x}) ): 这是空间单元 ( i ) 的属性值(如GDP、房价、发病率)与全局平均值 ( \bar{x} ) 的偏差。它衡量了单元 ( i ) 自身相对于整体水平是“高”还是“低”。这是一个标准化前的中心化处理。
  2. ( \sum w_{ij}(x_j - \bar{x}) ): 这是对单元 ( i ) 的所有邻居 ( j ) 的偏差进行加权求和。权重 ( w_{ij} ) 来自我们预先定义的空间权重矩阵 ( W ),它量化了单元 ( i ) 与每个邻居 ( j ) 的“邻近”程度。这部分代表了单元 ( i ) 的“邻里环境”。如果邻居们的值普遍高于均值,这个和为正;普遍低于均值,则为负。
  3. ( S^2 ): 这是所有单元属性值的方差。它是一个全局的标准化因子,目的是消除数据自身量纲和离散程度的影响,使得不同数据集计算出的 ( I_i ) 具有可比性。

所以,局部莫兰指数 ( I_i ) 的本质,可以理解为“单元自身偏离均值的程度”“其邻居平均偏离均值的程度”的乘积,再经过全局方差标准化。这个乘积关系至关重要:

  • 如果自身高,邻居们也高(两者皆正),乘积为正且较大,指示高-高聚集
  • 如果自身低,邻居们也低(两者皆负),负负得正,乘积也为正且较大,指示低-低聚集
  • 如果自身高,但邻居们都低(一正一负),乘积为负,指示高-低异常
  • 如果自身低,但邻居们都高(一负一正),乘积同样为负,指示低-高异常

注意:这里对“正负”的解读是基于上述乘积逻辑的简化。实际上,( I_i ) 的正负直接由这个乘积项决定,而其数值大小则反映了这种局部空间关联的强度。

2.2 空间关联的四种基本类型

基于局部莫兰指数 ( I_i ) 的符号(正/负)及其统计显著性,我们可以将每个空间单元归类到四种典型的局部空间关联模式中。这通常通过绘制“莫兰散点图”“LISA聚类图”来可视化。

  1. 高-高聚集 (HH): 单元自身属性值高,其周边邻居的属性值也高。这是典型的“热点区”。例如,城市中心商务区及其周边形成的房价高地。
  2. 低-低聚集 (LL): 单元自身属性值低,其周边邻居的属性值也低。这是典型的“冷点区”。例如,经济发展相对滞后的连片乡村地区。
  3. 高-低异常 (HL): 单元自身属性值高,但其周边邻居的属性值低。这表现为一个高值点被低值区域包围,是“鹤立鸡群”式的空间异常值。例如,一个建立在欠发达地区的豪华度假村。
  4. 低-高异常 (LH): 单元自身属性值低,但其周边邻居的属性值高。这表现为一个低值点被高值区域包围,是“洼地”式的空间异常值。例如,城市繁华商圈中的一个待改造的老旧小区。

不显著:除了以上四种,还有大量单元的局部莫兰指数在统计上不显著,意味着其与周边区域没有表现出稳定的、非随机的空间关联模式,可以视为随机分布。

2.3 与全局莫兰指数的关系

一个常见的误解是,将局部莫兰指数简单视为全局指数的分解。它们的关系需要辩证看待:

  • 求和关系:在数学上,所有局部莫兰指数 ( I_i ) 的加权平均值(通常以某种形式)等于全局莫兰指数 ( I )。这体现了全局模式是局部模式的综合体现。
  • 解释的独立性绝不能用全局指数的正负来直接推断局部模式。一个正的全局莫兰指数,可能源于少数几个强烈的HH集群和LL集群,同时混杂着大量的HL、LH异常点。反之,一个接近0的全局指数,也可能内部存在相互抵消的局部聚集模式。因此,局部分析是不可替代的。

3. 计算前的核心准备:空间权重矩阵构建

在计算局部莫兰指数之前,最基础、也最考验研究者空间思维的一步,就是构建合适的空间权重矩阵(Spatial Weight Matrix)( W )。这个 ( n \times n ) 的矩阵定义了所有空间单元之间的“邻居关系”及其强度,它直接决定了“局部”的范围和影响力大小。选择不当,会直接导致分析结果失真。

3.1 常见空间权重定义方法

  1. 邻接权重(Contiguity Weight):

    • Rook邻接:仅共享边(边界)的单元被视为邻居。像国际象棋中的“车”。
    • Queen邻接:共享边或顶点的单元都被视为邻居。像国际象棋中的“后”,更常用,因为它能捕捉到对角接触。
    • 操作:通常将邻居权重设为1,非邻居设为0,然后进行行标准化(使每行权重之和为1)。
    • 适用场景:适用于具有清晰、完整边界的面状数据(如行政区划)。对于不规则形状或存在飞地的区域,需要谨慎处理。
  2. 距离权重(Distance-based Weight):

    • K最近邻(KNN):为每个单元选择距离最近的K个单元作为邻居。确保了每个单元都有相同数量的邻居,避免了边缘单元邻居过少的问题。
    • 距离阈值:设定一个临界距离 ( d ),凡质心距离小于 ( d ) 的单元互为邻居。距离越近,权重越大,常用反距离(如 ( 1/d_{ij} ) 或 ( 1/d_{ij}^2 ) )或距离衰减函数(如高斯核函数)来定义权重。
    • 适用场景:适用于点数据或面状数据但关注的是质心距离。需要根据研究问题的空间尺度合理设定K值或阈值 ( d )。
  3. 经济社会权重:

    • 基于非地理距离,如经济联系强度、交通流量、人口流动量等构建权重矩阵。这更符合某些社会经济现象的空间相互作用逻辑。
    • 挑战:这类数据往往难以获取,且矩阵可能不对称(从i到j的流量不等于从j到i的)。

3.2 权重矩阵的标准化

原始的权重矩阵通常需要标准化,以方便解释和比较。最常用的是行标准化(Row Standardization),即将权重矩阵 ( W ) 的每一行元素除以其行和 ( \sum_j w_{ij} )。这样,标准化后的权重 ( w_{ij}^* ) 满足 ( \sum_j w_{ij}^* = 1 )。

  • 意义:行标准化后,局部莫兰指数公式中的加权求和 ( \sum w_{ij}(x_j - \bar{x}) ) 就变成了邻居偏差的加权平均值。这使得 ( I_i ) 更容易解释:它衡量的是单元自身偏差与邻居平均偏差的关系。
  • 注意:行标准化会改变权重的对称性,并可能影响统计推断。在计算显著性时,软件通常会处理这一点。

3.3 实操心得与避坑指南

  • “邻居”定义是主观的,但必须有据:没有绝对正确的权重矩阵。你的选择必须基于对研究现象的理论认知空间过程的理解。例如,研究空气污染扩散,风向和距离衰减可能是关键;研究零售店竞争,步行时间或车行时间可能比直线距离更合适。在论文中必须详细说明并论证你的权重选择。
  • 敏感性分析必不可少:不要只依赖一种权重矩阵做结论。尝试2-3种不同的定义方式(如Queen邻接、K=4的KNN、30公里阈值反距离),观察LISA聚类图的核心热点/冷点区域是否稳定。如果主要结论在不同权重下都成立,那么你的发现就更稳健。
  • 处理岛屿单元:有些空间单元可能没有邻居(如孤岛)。在邻接权重下,这些单元的行和为零,无法行标准化。通常的解决方法是将其排除在局部莫兰分析之外,或在距离权重中为其指定一个足够大的阈值以确保有邻居。
  • 空间权重矩阵的存储与检查:在Python的libpysal或R的spdep包中,权重对象创建后,务必使用print()或相关摘要函数查看基本信息,如最小/最大邻居数、岛屿单元数量、权重分布等,确保其符合你的预期。

4. 局部莫兰指数的计算与显著性检验

理解了原理和准备好了权重矩阵,我们就可以进入实际计算阶段。这里以Python的esda库和R的spdep库为例,详解计算流程和背后的统计检验逻辑。

4.1 基于Python (PySAL库) 的完整计算流程

假设我们有一个GeoDataFramegdf,包含几何列geometry和我们要分析的属性列value

import libpysal as ps from esda.moran import Moran_Local import geopandas as gpd import numpy as np import matplotlib.pyplot as plt import seaborn as sns # 1. 读取空间数据 gdf = gpd.read_file('your_data.geojson') # 确保数据有正确的投影(用于距离计算)或已是地理坐标系 # gdf = gdf.to_crs('EPSG:xxxx') # 2. 构建空间权重矩阵 (以Queen邻接为例) w = ps.weights.Queen.from_dataframe(gdf) # 检查权重矩阵 print(f"岛屿单元数量: {w.islands}") # 希望为 [] print(f"邻居数示例: {w.cardinalities}") # 查看前几个单元的邻居数 # 3. 计算局部莫兰指数 # 提取属性值数组,确保没有NaN值 y = gdf['value'].values # 计算Moran_Local对象 moran_local = Moran_Local(y, w, transformation='r', permutations=9999) # transformation='r': 行标准化权重 # permutations=9999: 使用9999次随机置换进行显著性检验 # 4. 将结果存回GeoDataFrame gdf['I'] = moran_local.Is # 局部莫兰指数值 gdf['p_sim'] = moran_local.p_sim # 基于随机置换的p值 gdf['z_sim'] = moran_local.z_sim # 基于随机置换的z得分(近似标准正态) # 5. 根据p值和I值判断聚类类型 (α=0.05) # 首先创建一个默认类别‘不显著’ gdf['cluster_type'] = '不显著' # 定义显著性水平 alpha = 0.05 # 判断条件 signif = gdf['p_sim'] < alpha gdf.loc[signif & (moran_local.q==1), 'cluster_type'] = '高-高' # q=1: 自身高,邻居高 gdf.loc[signif & (moran_local.q==2), 'cluster_type'] = '低-低' # q=2: 自身低,邻居低 gdf.loc[signif & (moran_local.q==3), 'cluster_type'] = '低-高' # q=3: 自身低,邻居高 gdf.loc[signif & (moran_local.q==4), 'cluster_type'] = '高-低' # q=4: 自身高,邻居低

关键参数解析

  • permutations: 这是显著性检验的蒙特卡洛模拟次数。默认9999次是一个较稳健的选择。次数越多,p值估计越稳定,但计算量越大。对于探索性分析,999或1999次也可接受。
  • transformation:'r'代表行标准化。这是最常用的设置。
  • moran_local.q: 这个属性直接给出了基于莫兰散点图四个象限的类别(1:HH, 2:LH, 3:LL, 4:HL)。注意,esda中的象限编号可能与某些文献的惯例不同,务必结合I值和实际数据检查。

4.2 基于R (spdep包) 的完整计算流程

library(sf) library(spdep) library(ggplot2) # 1. 读取空间数据 sf_data <- st_read('your_data.geojson') # 2. 构建空间权重矩阵 (以Queen邻接为例) # 首先获取多边形邻接关系 nb <- poly2nb(sf_data, queen=TRUE) # 检查是否有岛屿单元 which(card(nb) == 0) # 创建行标准化权重列表对象 w_list <- nb2listw(nb, style="W") # style="W" 代表行标准化 # 3. 计算局部莫兰指数 local_moran <- localmoran(sf_data$value, w_list, alternative="two.sided", nsim=9999) # alternative="two.sided": 双尾检验 # nsim=9999: 随机置换次数 # 4. 将结果合并到数据中 sf_data$I_i <- local_moran[, "Ii"] # 局部莫兰指数值 sf_data$p_val <- local_moran[, "Pr(z != E(Ii))"] # 近似正态假设下的p值 # 注意:localmoran默认提供的是基于正态近似的p值。若需模拟p值,需使用`localmoran_perm` # 5. 计算期望值和方差,用于判断象限 (HH, HL, LH, LL) y <- sf_data$value y_lag <- lag.listw(w_list, y) # 计算空间滞后项(邻居的加权平均) mean_y <- mean(y) # 中心化 y_centered <- y - mean_y y_lag_centered <- y_lag - mean_y # 6. 根据中心化值和空间滞后项中心化值判断类型 sf_data$cluster_type <- "不显著" signif <- sf_data$p_val < 0.05 sf_data$cluster_type[signif & y_centered > 0 & y_lag_centered > 0] <- "高-高" sf_data$cluster_type[signif & y_centered < 0 & y_lag_centered < 0] <- "低-低" sf_data$cluster_type[signif & y_centered < 0 & y_lag_centered > 0] <- "低-高" sf_data$cluster_type[signif & y_centered > 0 & y_lag_centered < 0] <- "高-低"

4.3 显著性检验:理解p值的来源

为什么需要显著性检验?因为即使数据完全随机分布,由于随机波动,也可能计算出非零的 ( I_i ) 值。检验的目的是判断观察到的 ( I_i ) 是否显著区别于“空间随机”的零假设。

两种主要方法:

  1. 基于正态分布的近似检验:在数据满足一定条件(如正态性)且样本量较大时,可以推导出 ( I_i ) 的近似抽样分布,进而计算z得分和p值。R的localmoran()默认输出此结果。缺点:对数据分布敏感,在样本量小或数据非正态时可能不可靠。
  2. 基于随机置换的蒙特卡洛模拟:这是更稳健、更推荐的方法。其步骤如下:
    • 零假设:属性值在空间上的分布是随机的。
    • 操作:保持空间单元位置和权重矩阵不变,将属性值随机打乱(置换)并重新计算 ( I_i ),重复多次(如9999次)。
    • 构建经验分布:这9999次随机置换计算出的 ( I_i ) 值,构成了在零假设下 ( I_i ) 可能取值的经验分布。
    • 计算模拟p值:将实际观测到的 ( I_i ) 值与这个经验分布比较。如果观测值落在经验分布的极端位置(如最高的2.5%或最低的2.5%),则认为它不太可能由随机过程产生,从而拒绝零假设。p值 = (排名位置) / (置换次数+1)。

实操心得务必使用蒙特卡洛模拟(permutations参数)。正态近似虽然快,但其假设在实际数据中常常不成立。9999次置换在现代计算机上对于中等规模数据(如几百个面单元)通常可以在可接受时间内完成。报告结果时,应注明使用的检验方法和置换次数。

5. 结果可视化:莫兰散点图与LISA聚类图

计算出结果后,我们需要两种核心图形来理解和展示它。

5.1 莫兰散点图

莫兰散点图是理解局部空间自相关最直观的工具。

  • 横轴 (x): 标准化后的原始属性值 ( z )。
  • 纵轴 (y): 标准化后的空间滞后值 ( W_z )(即邻居加权平均值)。
  • 四个象限:分别对应HH(第一象限)、LH(第二象限)、LL(第三象限)、HL(第四象限)。
  • 斜率:所有点的回归线斜率,等于全局莫兰指数
# Python 绘制莫兰散点图 fig, ax = plt.subplots(figsize=(10, 8)) # 绘制散点 scatter = ax.scatter(moran_local.z, moran_local.wz, c=moran_local.q, cmap='Set1', alpha=0.8, edgecolors='k') # 添加象限线 ax.axhline(y=0, color='k', linestyle='--', alpha=0.5) ax.axvline(x=0, color='k', linestyle='--', alpha=0.5) # 添加回归线 b, a = np.polyfit(moran_local.z, moran_local.wz, 1) ax.plot(moran_local.z, a + b * moran_local.z, 'r-', label=f'斜率 (Global I): {b:.3f}') ax.set_xlabel('标准化属性值 (z)') ax.set_ylabel('标准化空间滞后值 (Wz)') ax.set_title('局部莫兰散点图 (LISA)') ax.legend() plt.colorbar(scatter, ax=ax, label='象限 (1:HH,2:LH,3:LL,4:HL)') plt.show()

解读技巧:观察点主要聚集在哪个象限。如果大量点聚集在一三象限(HH, LL),说明存在明显的空间正相关(聚集)。如果大量点分布在二四象限(LH, HL),则说明存在空间负相关(异常)。点的离散程度也反映了局部关联强度的差异。

5.2 LISA聚类图

LISA聚类图是将显著性检验和聚类类型结果直接映射到地理空间上,是最具洞察力的成果图。

# Python 使用geopandas绘制LISA聚类图 fig, ax = plt.subplots(1, 1, figsize=(12, 10)) # 定义颜色映射,为四种显著类型和不显著分别指定颜色 cluster_colors = { '高-高': 'red', '低-低': 'blue', '低-高': 'lightblue', '高-低': 'pink', '不显著': 'lightgrey' } gdf['plot_color'] = gdf['cluster_type'].map(cluster_colors) gdf.plot(color=gdf['plot_color'], ax=ax, edgecolor='black', linewidth=0.5) # 添加图例 from matplotlib.patches import Patch legend_elements = [Patch(facecolor=v, edgecolor='k', label=k) for k, v in cluster_colors.items()] ax.legend(handles=legend_elements, loc='lower left', fontsize=10) ax.set_title('局部空间自相关 (LISA) 聚类图', fontsize=16) ax.set_axis_off() plt.show()

制图与解读要点

  • 颜色选择:遵循制图学惯例,HH(热点)常用红色系,LL(冷点)常用蓝色系,异常点用对比色(如粉色、浅蓝),不显著区域用中性灰色。这有助于读者快速建立认知。
  • 突出显示:可以仅对显著区域着色,不显著区域留白或浅灰,使地图焦点更清晰。
  • 结合背景:一定要将LISA聚类图叠加在基础地理底图或相关背景(如地形、主要道路、河流)上,结合地理背景解释集群的成因。例如,HH集群是否沿交通干线分布?LL集群是否与某种地理限制(如山区)重合?
  • 多图对比:将不同年份或不同变量的LISA图并排展示,可以直观分析时空演变模式。

6. 高级议题与常见陷阱

掌握了基础计算和可视化后,要产出严谨的研究,还必须关注以下高级议题和常见陷阱。

6.1 多重比较问题与修正

当我们对成百上千个空间单元分别进行局部莫兰指数检验时,实际上是在同时进行大量(m次)的假设检验。这会导致多重比较问题:即使所有局部空间关联都是随机的(零假设为真),我们仍可能期望有大约 ( \alpha \times m ) 个单元被错误地判定为显著(第一类错误膨胀)。例如,α=0.05时,检验1000个单元,即使没有真实模式,也可能有50个单元“假显著”。

常用修正方法

  1. Bonferroni校正:将显著性水平调整为 ( \alpha / m )。这是最严格的校正,但可能过于保守,导致很多真实模式被漏掉。
  2. 错误发现率控制:如Benjamini-Hochberg (BH) 方法。它控制的是在所有被拒绝的零假设中,错误拒绝的比例(FDR)。相比Bonferroni,它更灵活,功效更高,在空间分析中应用越来越广泛。
  3. 基于随机置换的联合参考分布:一些高级方法通过构建所有局部统计量的联合分布来评估整体显著性,但实现复杂。

实操建议:在学术研究中,必须报告你是否以及如何处理了多重比较问题。对于探索性分析,可以同时展示未校正和经过FDR校正的结果,并讨论其差异。在软件中,esda库的Moran_Local函数返回的p_sim是未校正的模拟p值,需要自行进行校正。

# Python示例:使用statsmodels进行FDR校正 from statsmodels.stats.multitest import multipletests # moran_local.p_sim 是原始的模拟p值数组 reject_fdr, pvals_fdr, _, _ = multipletests(moran_local.p_sim, alpha=0.05, method='fdr_bh') gdf['p_fdr'] = pvals_fdr gdf['sig_fdr'] = reject_fdr # 然后基于gdf['sig_fdr']和moran_local.q重新定义cluster_type_fdr

6.2 边缘效应与小样本问题

边缘效应是指位于研究区域边界的空间单元,由于其邻居数量可能少于内部单元,导致其局部统计量的计算不稳定、方差增大,进而影响显著性检验的功效。

  • 影响:边界上的单元更容易出现假阴性(真实关联但检验不显著)或假阳性。
  • 缓解方法
    • 使用K最近邻(KNN)权重,确保每个单元都有相同数量的邻居。
    • 在可能的情况下,扩大研究区域范围,将分析区域置于一个更大的缓冲区中(尽管只分析核心区),为边界单元提供更完整的邻居环境。
    • 在解释结果时,对边界区域的显著性格外谨慎,可以结合地理背景判断。

小样本问题是指当研究区域内单元总数很少时(如少于30),基于渐近分布的统计检验可能失效,蒙特卡洛模拟的稳定性也会变差。

  • 对策必须使用蒙特卡洛模拟,并增加置换次数(如99999次)。同时,应明确承认小样本的局限性,对结果的解读保持保守,更多作为描述性探索而非严格的统计推断。

6.3 空间权重矩阵的敏感性再探讨

这是分析中最关键的假设之一,其影响怎么强调都不为过。

  • 不同权重导致不同结果:一个区域用Queen邻接可能是HH集群,用距离阈值可能就不显著了。这并非说明方法无效,而是揭示了空间关联的“尺度依赖性”。
  • 系统性的敏感性分析流程
    1. 定义2-3套合理的权重方案。例如:方案A(Queen邻接)、方案B(K=8的KNN)、方案C(50km距离阈值高斯核)。
    2. 并行计算每种方案下的局部莫兰指数和LISA图。
    3. 比较核心结论:主要的热点区(HH)、冷点区(LL)在空间位置和范围上是否一致?如果几种方案都识别出同一片区域是热点,那么这个热点就非常稳健。如果结果差异很大,你需要回到理论层面,思考哪种“邻近”定义更符合你研究的现象。
    4. 报告时,应展示主要方案的结果,并在方法部分或附录中说明敏感性分析的情况。

6.4 与其他局部空间统计量的关系

局部莫兰指数是LISA家族中最常用的一员,但并非唯一。

  • 局部Getis-Ord Gi/Gi*:这也是识别热点和冷点的常用指标。与局部莫兰指数的关键区别在于:
    • 莫兰I:基于与全局均值的偏差,能识别所有四种模式(HH, HL, LH, LL)。
    • Getis-Ord Gi*:基于局部总和,主要对识别HH和LL集群敏感,对HL和LH异常不敏感。它的零假设是“空间无聚集”,计算时包含了自身值。
    • 如何选择:如果你的研究问题明确关注“热点”和“冷点”,Gi*可能更直观。如果你想全面探测包括异常值在内的所有局部空间模式,局部莫兰指数更合适。在实践中,可以两者都计算,相互印证。

7. 完整案例实操:中国地级市人均GDP的局部空间分析

让我们通过一个假设的完整案例,串联所有步骤。假设我们分析2020年中国各地级市人均GDP的空间分布。

步骤1:数据准备与探索

  • 数据:包含中国地级市行政区划面数据的GeoDataFramegdf_china,以及per_gdp字段。
  • 检查:数据投影(确保为等面积投影如Albers,用于距离计算)、缺失值、数据分布(是否严重偏态?考虑对数变换)。

步骤2:构建空间权重矩阵

  • 考虑到中国地级市面积差异大,Queen邻接可能使新疆、内蒙古等地广人稀的城市邻居过少。我们选择两种方案对比:
    • 方案1:K最近邻,K=8。确保每个城市都有8个最近邻。
    • 方案2:距离阈值权重,阈值设为相邻城市质心距离的80分位数(如300公里),使用高斯核函数衰减。
  • 使用libpysal创建这两个权重矩阵w_knnw_dist

步骤3:计算与检验

  • 分别对w_knnw_dist计算Moran_Local,置换次数9999。
  • 对模拟p值进行FDR校正(BH方法)。
  • 根据校正后的显著性(α=0.05)和象限划分聚类类型。

步骤4:可视化与解读

  • 绘制两张LISA聚类地图(分别对应两种权重方案)。
  • 核心发现可能包括
    • 稳健的HH集群:长三角、珠三角、京津冀城市群在两种权重下都显示为显著的热点区,表明这些区域内部及周边城市形成了高人均GDP的紧密集聚。
    • 稳健的LL集群:西南部分省份的连绵地区在两种权重下均显示为冷点区。
    • 敏感区域:一些中部城市在KNN权重下是HH,在距离权重下不显著。这可能意味着这些城市的高GDP主要依赖于与少数几个特定邻近城市的联系(被KNN捕捉),而非广泛的区域辐射(距离权重要求更广泛的邻近性)。
    • 异常点:可能发现个别资源型城市(如鄂尔多斯)呈现HL模式(自身高,周边低),或个别发达区域内的欠发达县区呈现LH模式。

步骤5:报告与深化

  • 在报告中,首先展示基于更稳健权重(如KNN)的LISA图作为主结果。
  • 在方法部分详细说明权重选择、敏感性分析和多重比较校正过程。
  • 在讨论部分,结合地理、经济政策(如西部大开发、中部崛起)、交通基础设施(高铁网络)等,解释HH和LL集群的形成原因,并探讨HL/LH异常点的特殊性。
  • 可以进一步将局部莫兰指数作为因变量或自变量,构建空间计量经济模型,探究其影响因素或经济后果。

最后的个人体会:局部莫兰指数是一个强大的“空间显微镜”,但它给出的是一幅静态的、描述性的快照。真正的分析功力,体现在权重矩阵构建的理论依据、对多重比较等统计问题的严谨处理、对敏感性分析结果的合理解读,以及最终将统计模式与深层次的地理、经济、社会过程联系起来的能力。记住,软件输出的是“数字”和“颜色”,而研究者需要讲述的是“故事”和“机理”。每次分析前,多花时间思考你的“空间”究竟是如何定义的,这能避免很多后续的麻烦。