截面数据空间计量模型：SAR、SEM、SDM估计与Python/R实战

简介：这份资源面向空间计量经济学初学者与实证研究者，系统整理了截面数据下的主流空间回归估计方法，帮助解决空间依赖建模与模型选择问题。内容覆盖空间滞后模型SLM、空间误差模型SEM、空间杜宾模型SDM及其误差形式，并延伸至自变量空间滞后模型、Kelejian-Prucha模型、一般嵌套空间模型、空间扩展模型与地理加权回归，同时包含拉格朗日乘子检验LM等配套检验，基本构成一套完整的空间计量方法工具箱。压缩包约6.15MB，内含MATLAB脚本与原始数据，并附带jplv7、Elhorst_codes等常用工具包，便于直接运行与二次修改。已有2073人学习下载，适合需要复现经典模型、对比不同估计量或撰写论文实证章节的读者参考。

1. 空间计量模型估计方法：截面数据里那套绕不开的模型体系

手里有一份带地理坐标的截面数据，跑普通 OLS 发现残差莫兰指数显著不为零，这时候空间计量就该上场了。空间滞后、空间误差、杜宾模型，这三个名字几乎出现在每一篇区域经济、房价、环境规制的实证论文里。截面数据没有时间维度，估计量的渐近性质、工具变量的构造、权重矩阵的设定，每一步都比面板数据更敏感。这篇笔记面向已经会用 Stata 或 R 跑回归、但面对空间权重矩阵和拉格朗日乘子检验一头雾水的从业者。我会把从权重矩阵构建到模型选择的完整链路拆开，附上可直接跑的代码和一份模拟截面数据，让你在自己的数据上复现这套流程。空间计量不是玄学，但确实有几个参数设错就全盘翻车的地方，后面会逐个点出来。

2. 截面数据空间计量的理论底座：为什么 OLS 在这里失效

2.1 空间依赖的两种来源与模型对应关系

截面数据里，空间依赖主要来自两个方向。一个是实质性的溢出效应：一个地区的被解释变量直接受到邻近地区被解释变量的影响，比如某城市的房价上涨会带动周边城市房价，这种机制对应空间滞后模型（SAR），形式为 y = ρWy + Xβ + ε。另一个是误差项的空间自相关：遗漏变量本身在空间上聚集，导致残差不再独立，对应空间误差模型（SEM），形式为 y = Xβ + u，u = λWu + ε。杜宾模型（SDM）则是两者的综合，既包含 Wy 也包含 WX，形式为 y = ρWy + Xβ + WXθ + ε。

选哪个模型不是拍脑袋决定的。常见做法是先跑 OLS，然后做拉格朗日乘子检验（LM 检验）及其稳健版本。LM-lag 显著而 LM-error 不显著，倾向 SAR；反过来倾向 SEM；两个都显著，考虑 SDM 或者用稳健 LM 统计量再判断。这个流程在 Anselin 的经典文献里有详细推导，我一般会把它作为模型选择的起点，而不是终点。

2.2 空间权重矩阵的构建：从邻接到距离衰减

权重矩阵 W 是空间计量的核心输入，它决定了“谁影响谁”以及“影响多大”。截面数据常用的 W 有三类：

行标准化是标配操作，让每行权重之和为 1，这样 ρ 和 λ 的解释才是“邻居的加权平均”。不标准化的话，系数会随邻居数量变化，跨区域比较就失去意义。

import numpy as np import pandas as pd from scipy.spatial.distance import cdist # 假设 df 包含 longitude, latitude, y, x1, x2 coords = df[['longitude', 'latitude']].values n = len(df) # 距离矩阵 D = cdist(coords, coords, metric='euclidean') # 距离衰减权重，截断距离取 100 公里（根据坐标单位调整） threshold = 1.0 # 若坐标已投影为公里，此处为 100；若为经纬度需先转换 W = np.where((D > 0) & (D < threshold), 1.0 / D, 0.0) # 行标准化 row_sums = W.sum(axis=1, keepdims=True) W_std = np.divide(W, row_sums, out=np.zeros_like(W), where=row_sums != 0) # 检查是否有孤立点 isolated = np.where(row_sums.flatten() == 0)[0] print(f"孤立点索引: {isolated}")

这段代码先算欧氏距离矩阵，再按距离倒数构造权重，截断距离 threshold 需要根据你的坐标单位来定。如果经纬度直接算距离，threshold 设 1.0 大约对应 111 公里，但更稳妥的做法是先用投影转换把经纬度转成平面坐标。行标准化用np.divide配合where避免除零。最后检查孤立点，如果有地区没有任何邻居，后续估计会报错，需要手动给它指定最近邻或者从样本中剔除。

2.3 截面数据估计方法的选择：ML、GMM 与贝叶斯

截面数据的空间模型估计，主流有三条路。最大似然（ML）是默认选项，Stata 的spregress和 R 的spdep都支持，性质好但计算量大，n 超过几千时雅可比矩阵的行列式计算会拖慢速度。广义矩估计（GMM）用 Wy 的高阶滞后作为工具变量，计算快，但工具变量的有效性依赖权重矩阵设定正确。贝叶斯方法通过 MCMC 抽样，适合小样本和复杂层级结构，但调参和收敛诊断需要经验。

我一般会先用 ML 跑基准结果，再用 GMM 做稳健性检验。如果两者系数符号和显著性差异很大，说明权重矩阵可能设错了，或者模型存在设定偏误。截面数据没有时间维度，无法用固定效应吸收空间异质性，所以对权重矩阵的敏感性比面板数据更高。

3. 用 Python 和 R 跑通 SAR、SEM、SDM 的完整流程

3.1 数据准备与探索性空间数据分析

先加载数据，做探索性空间数据分析（ESDA）。莫兰指数是最基础的全局自相关指标，公式是 I = (n/S0) * (y'Wy / y'y)，其中 S0 是权重矩阵所有元素之和。莫兰指数显著为正，说明高值和高值聚集，负值说明高值和低值相邻。

import libpysal from esda.moran import Moran import geopandas as gpd # 读取 shapefile 或带坐标的 csv gdf = gpd.read_file('county_data.shp') # 构建邻接权重 w = libpysal.weights.Queen.from_dataframe(gdf) w.transform = 'r' # 行标准化 # 计算莫兰指数 y = gdf['y'].values moran = Moran(y, w) print(f"Moran's I: {moran.I:.4f}, p-value: {moran.p_sim:.4f}")

Queen.from_dataframe基于共享边界或顶点构建邻接关系，w.transform = 'r'执行行标准化。moran.p_sim是 999 次随机置换得到的伪 p 值，比正态近似更稳健。如果莫兰指数不显著，空间计量可能不是必需的，但也不绝对——局部自相关可能被全局指标掩盖，可以进一步看局部莫兰指数（LISA）聚类图。

3.2 SAR 模型的 ML 估计与参数解读

空间滞后模型用spregress在 Stata 里一行命令就能跑，Python 这边用spreg库。下面用 ML 估计 SAR：

from spreg import ML_Lag # y 和 X 准备好，w 是行标准化后的权重矩阵 model_sar = ML_Lag(y, X, w=w, name_y='y', name_x=['x1', 'x2']) print(model_sar.summary)

输出里重点关注三个部分。一是 ρ（空间自回归系数），显著为正说明邻居的 y 对本地区 y 有正向溢出。二是 β 系数，解释时要注意 SAR 模型存在反馈效应：X 对 y 的边际影响不是 β，而是 (I - ρW)^(-1)β 的对角线元素，直接拿 β 说事会低估直接效应、忽略间接效应。三是伪 R² 和对数似然值，用于模型比较。

我一般会额外计算直接效应、间接效应和总效应。直接效应是 X 变化对本地区 y 的平均影响，间接效应是通过邻居反馈回来的影响，总效应是两者之和。spreg不直接输出这些，需要手动算：

import numpy as np rho = model_sar.rho I = np.eye(n) S = np.linalg.inv(I - rho * W_std) direct = np.mean(np.diag(S)) * model_sar.betas[1] # 以 x1 为例 total = np.mean(S) * model_sar.betas[1] indirect = total - direct print(f"直接效应: {direct:.4f}, 间接效应: {indirect:.4f}, 总效应: {total:.4f}")

这段代码里S是空间乘数矩阵，np.diag(S)取对角线元素求平均得到直接效应，np.mean(S)是总效应，两者相减是间接效应。注意model_sar.betas[1]对应第一个解释变量的系数，索引要跟你的 X 列顺序对齐。

3.3 SEM 与 SDM 的估计及 LR 检验对比

空间误差模型用ML_Error，杜宾模型用ML_Lag加上 WX 项。SDM 的估计可以手动构造 WX 矩阵后放进ML_Lag，也可以用spreg的ML_Lag配合扩展矩阵。

from spreg import ML_Error # SEM model_sem = ML_Error(y, X, w=w, name_y='y', name_x=['x1', 'x2']) print(f"SEM lambda: {model_sem.lam:.4f}, Log-L: {model_sem.logll:.2f}") # SDM: 构造 [X, WX] WX = W_std @ X X_sdm = np.hstack([X, WX]) model_sdm = ML_Lag(y, X_sdm, w=w, name_y='y', name_x=['x1', 'x2', 'Wx1', 'Wx2']) print(f"SDM rho: {model_sdm.rho:.4f}, Log-L: {model_sdm.logll:.2f}")

模型比较用似然比检验（LR）。SDM 嵌套 SAR（原假设 θ=0）和 SEM（原假设 θ + ρβ = 0），所以可以用 LR 统计量判断 SDM 是否可以简化为 SAR 或 SEM。LR = 2*(logL_SDM - logL_restricted)，自由度是约束个数。如果 LR 不显著，选简洁的 SAR 或 SEM；显著则保留 SDM。

from scipy.stats import chi2 lr_sar = 2 * (model_sdm.logll - model_sar.logll) p_sar = 1 - chi2.cdf(lr_sar, df=2) # 两个 WX 项约束 print(f"LR test SDM vs SAR: {lr_sar:.2f}, p={p_sar:.4f}") lr_sem = 2 * (model_sdm.logll - model_sem.logll) p_sem = 1 - chi2.cdf(lr_sem, df=2) print(f"LR test SDM vs SEM: {lr_sem:.2f}, p={p_sem:.4f}")

自由度取约束个数，SDM 比 SAR 多两个 WX 系数，所以 df=2。如果 p 值小于 0.05，拒绝原假设，SDM 更合适。这套检验流程在 LeSage 和 Pace 的教材里有完整推导，我一般会把它作为模型选择的最终依据，而不是只看 LM 检验。

4. 避坑与排查：截面空间计量里那些让人翻车的细节

4.1 权重矩阵行标准化后出现孤立点

现象：估计时提示“matrix is singular”或者 ρ 的估计值接近 1，标准误爆炸。原因：某些地区没有邻居，行标准化后整行全为零，W 矩阵出现零行，导致 (I - ρW) 不可逆。解决：在构建 W 时检查row_sums，对孤立点手动指定最近邻，或者用 k 近邻权重替代距离阈值权重。k 近邻保证每个地区至少有 k 个邻居，不会出现零行。

4.2 莫兰指数显著但回归残差仍有空间自相关

现象：跑了 SAR 或 SEM，残差的莫兰指数依然显著。原因：模型设定不完整，可能遗漏了 WX 项（应该用 SDM），或者权重矩阵与实际溢出机制不匹配。解决：先跑 SDM 看 WX 系数是否显著，如果显著则保留 SDM；如果不显著但残差仍有自相关，尝试换权重矩阵，比如从邻接换成距离衰减，或者用经济距离矩阵（如 GDP 倒数）做稳健性检验。

4.3 直接效应和间接效应的标准误无法直接获取

现象：spreg只输出 β 和 ρ 的协方差矩阵，直接效应和间接效应是 β 和 ρ 的非线性函数，标准误需要 Delta 方法或 Bootstrap。原因：很多人直接拿 β 的显著性说事，忽略了反馈效应。解决：用 Bootstrap 重抽样计算效应分布。每次重抽样后重新估计模型，计算直接效应和间接效应，重复 500 次取标准差。计算量大但结果可靠。

from numpy.random import choice def bootstrap_effects(y, X, W, n_boot=500): effects = [] for _ in range(n_boot): idx = choice(len(y), len(y), replace=True) # 注意：重抽样后需重新构建 W 的子矩阵，此处简化示意 model = ML_Lag(y[idx], X[idx], w=W[np.ix_(idx, idx)]) rho = model.rho S = np.linalg.inv(np.eye(len(idx)) - rho * W[np.ix_(idx, idx)]) direct = np.mean(np.diag(S)) * model.betas[1] total = np.mean(S) * model.betas[1] effects.append([direct, total - direct]) effects = np.array(effects) return effects.mean(axis=0), effects.std(axis=0)

这段代码是简化示意，实际重抽样时权重矩阵的子矩阵需要重新行标准化，否则行和不为 1。Bootstrap 的另一个坑是重抽样后可能出现新的孤立点，需要在循环里加判断。

4.4 截面数据样本量过小导致 ML 估计不收敛

现象：n 小于 50 时，ML 估计的 ρ 标准误很大，或者优化算法不收敛。原因：截面数据的渐近性质依赖大样本，小样本下雅可比矩阵的行列式计算不稳定。解决：改用贝叶斯估计，给 ρ 设一个合理的先验（如均匀分布 [-1, 1]），通过 MCMC 抽样得到后验分布。R 的spBayes或 Python 的pymc都可以做。或者用 GMM，工具变量用 W²y、W³y，计算快且对小样本更稳健。

4.5 权重矩阵的截断距离选择影响结论

现象：距离阈值从 50 公里换到 100 公里，ρ 的显著性从 0.01 变成 0.10。原因：截断距离决定了邻居数量，距离太近邻居太少，估计不稳定；距离太远邻居太多，溢出效应被稀释。解决：做敏感性分析，画 ρ 随截断距离变化的曲线，选择 ρ 估计最稳定的一段。常见做法是取距离分布的 25%、50%、75% 分位数分别跑，如果结论一致则稳健，不一致则需在论文里说明局限性。

5. 进阶技巧：用贝叶斯方法处理小样本和模型不确定性

截面数据做空间计量，最头疼的是样本量不够大。全国 31 个省份、2800 多个县，听起来不少，但一旦按区域拆分或者加入交互项，自由度就紧张了。贝叶斯方法在这种情况下比 ML 更有优势，因为它不依赖大样本渐近理论，先验信息可以正则化估计。

我一般用 R 的spdep配合MCMCpack做贝叶斯 SAR。核心是给 ρ 设均匀先验，β 和 σ² 用默认的无信息先验，然后跑 10000 次 MCMC，丢弃前 2000 次作为 burn-in。收敛诊断看 Gelman-Rubin 统计量，跑三条链，R-hat 小于 1.1 才算收敛。

library(spdep) library(MCMCpack) # 假设 listw 是行标准化后的权重列表 # y 和 X 已准备好 n <- length(y) I <- diag(n) W <- listw2mat(listw) # 对数似然函数 loglik <- function(rho, y, X, W) { A <- I - rho * W e <- A %*% y - X %*% solve(t(X) %*% X) %*% t(X) %*% A %*% y sigma2 <- sum(e^2) / n log(det(A)) - (n/2) * log(sigma2) } # MCMC 抽样 post <- MCMCmetrop1R(loglik, theta.init = 0.3, y = y, X = X, W = W, mcmc = 10000, burnin = 2000, thin = 1) summary(post)

这段 R 代码用MCMCmetrop1R对 ρ 做 Metropolis 抽样，loglik函数里log(det(A))是雅可比项，sigma2是误差方差。theta.init给 ρ 一个初始值，一般设 0.3 左右。跑完后summary(post)给出后验均值和置信区间。如果 ρ 的后验分布集中在 0 附近，说明空间自相关弱，SAR 可能不比 OLS 好多少。

贝叶斯方法的另一个好处是可以做模型平均。把 SAR、SEM、SDM 分别跑一遍，用贝叶斯因子或者 WAIC 比较，按后验概率加权平均预测值。这样避免了“选一个模型然后假装它是对的”这种常见问题。截面数据没有时间维度，模型不确定性比面板数据更大，模型平均是一个值得投入的方向。

最后说一个我自己的习惯：每次跑完空间计量，我都会把残差的莫兰指数再算一遍，画一张残差的空间分布图。如果残差还有聚集，说明模型没榨干空间信息，要么加 WX 项，要么换权重矩阵。这个习惯帮我省了很多后悔药，也让我在审稿人问“为什么不用 SDM”的时候有底气回答。希望帮到你。

本文还有配套的精品资源，点击获取