Applied Spatial Statistics（八）：Python 空间回归实战——从 Moran‘s I 到空间滞后效应分解

1. 空间回归分析入门从Morans I到模型选择第一次接触空间回归时我被那些专业术语搞得晕头转向。直到用真实数据跑了一遍完整流程才发现原来空间统计分析可以这么直观。咱们就以2020年美国大选数据为例手把手教你用Python的pysal库完成从空间自相关检验到模型选择的完整分析。先说说为什么要做空间回归。传统回归分析有个重要假设样本之间相互独立。但地理数据往往打脸 - 相邻地区的经济指标、人口特征都存在明显相关性。我刚开始做县域经济分析时就犯过这个错用普通OLS回归结果残差图一看好家伙高值区和低值区都扎堆出现这明显违背了独立性假设。空间回归主要解决两类问题空间误差模型(SEM)处理被忽略的空间依赖性空间滞后模型(SLM)研究空间溢出效应安装环境特别简单pip install pysal esda splot geopandas数据准备阶段要注意几个坑空间权重矩阵要提前标准化岛屿地区(没有邻居的单元)需要特殊处理变量最好先做标准化# 数据加载与预处理示例 import geopandas as gpd from libpysal.weights import Queen # 读取选举数据和地理边界 voting gpd.read_file(https://raw.github.com/Ziqi-Li/gis5122/master/data/voting_2020.csv) shp gpd.read_file(https://raw.github.com/Ziqi-Li/gis5122/master/data/us_counties.geojson) # 创建Queen邻接权重矩阵 w Queen.from_dataframe(shp, silence_warningsTrue) w.transform r # 行标准化2. 空间自相关检测Morans I实战Morans I是检测空间自相关的利器取值在-1到1之间接近1相似值聚集(高高或低低)接近-1相异值聚集(高低交错)接近0随机分布计算选举数据的Morans Ifrom esda.moran import Moran from splot.esda import plot_moran # 先跑普通OLS模型 ols spreg.OLS(y, X, name_ynew_pct_dem, name_xvariable_names) # 计算残差的空间自相关 ols_moran Moran(ols.u, w, permutations199) plot_moran(ols_moran, figsize(10,4))上次分析教育投入数据时Morans I值达到0.65p值0.001说明空间自相关极强。这时候再用传统回归就是掩耳盗铃了。空间自相关可能来自测量误差的空间传播遗漏了空间结构的变量真实的空间交互作用3. 模型选择SEM vs SLM当存在空间自相关时我们主要考虑两种模型3.1 空间误差模型(SEM)适用于误差项存在空间依赖的情况。公式为 y Xβ u, u λWu ε# SEM模型拟合 sem spreg.ML_Error(y, X, ww, name_xvariable_names, name_ynew_pct_dem) print(sem.summary) # 检查过滤后的残差 sem_moran Moran(sem.e_filtered, w, permutations199) plot_moran(sem_moran, zstandardTrue, figsize(10,4))3.2 空间滞后模型(SLM)研究空间溢出效应公式为 y ρWy Xβ ε# SLM模型拟合 slm spreg.ML_Lag(y, X, ww, name_ynew_pct_dem, name_xvariable_names) print(slm.summary) # 检查模型残差 slm_moran Moran(slm.u, w, permutations199) plot_moran(slm_moran, zstandardTrue, figsize(10,4))选择依据拉格朗日乘数检验(LM test)稳健LM检验模型拟合指标(AIC/BIC)残差的空间自相关程度在大选数据中误差模型的稳健LM统计量(1186)远高于滞后模型(18)且AIC更低因此SEM更合适。4. 空间效应分解直接效应与间接效应空间滞后模型最有趣的部分是可以量化空间溢出效应。以大学毕业生比例(pct_bach)对民主党得票率的影响为例# 计算总效应矩阵 rho slm.rho[0] # 获取估计的rho值 beta slm.betas[6] # pct_bach的系数 I np.identity(w.n) S np.linalg.inv(I - rho * w.full()[0]) * beta # 提取特定县的影响 leon_idx shp_voting[shp_voting[NAME_x] Leon].index[0] total_effect S[leon_idx, leon_idx] # 总效应 direct_effect beta # 直接影响 indirect_effect total_effect - direct_effect # 间接影响效应分解结果直接影响0.83莱昂县本地效应总效应0.90包含空间溢出间接效应0.07来自邻近县的影响可视化空间溢出效应shp_voting[leon_effect] S[:, leon_idx] ax shp_voting.plot(columnleon_effect, legendTrue, figsize(12,8)) plt.title(空间溢出效应莱昂县大学毕业生比例增加1%的影响)从图中可以清晰看出效应随距离衰减明显相邻县影响约0.1-0.2%远离县影响接近零5. 实战经验与常见问题踩过几次坑后我总结了一些实用经验权重矩阵选择Queen邻接 vs Rook邻接距离阈值法 vs K最近邻经济距离权重(如GDP差异)模型诊断# 空间Hausman检验 from spreg import hausman hausman(sem, slm) # 异方差检验 from spreg import BreuschPagan BreuschPagan(sem)性能优化大数据集使用稀疏矩阵考虑GMM估计替代ML并行化排列检验结果解释陷阱不要混淆空间滞后项和空间误差项注意内生性问题空间效应的非线性检验最后提醒初学者几个常见错误忘记标准化权重矩阵忽略岛屿单元处理误读空间滞后项系数忽视异方差的影响空间回归是个强大的工具但需要结合地理学理论和实际业务知识来解释结果。我在分析房价数据时就犯过单纯依赖统计指标的错后来结合区位理论才发现模型漏掉了关键的空间交互特征。