简介这份资源面向空间计量经济学初学者与实证研究者系统整理了截面数据下的主流空间回归估计方法帮助解决空间依赖建模与模型选择问题。内容覆盖空间滞后模型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β uu λWu ε。杜宾模型SDM则是两者的综合既包含 Wy 也包含 WX形式为 y ρWy Xβ WXθ ε。选哪个模型不是拍脑袋决定的。常见做法是先跑 OLS然后做拉格朗日乘子检验LM 检验及其稳健版本。LM-lag 显著而 LM-error 不显著倾向 SAR反过来倾向 SEM两个都显著考虑 SDM 或者用稳健 LM 统计量再判断。这个流程在 Anselin 的经典文献里有详细推导我一般会把它作为模型选择的起点而不是终点。2.2 空间权重矩阵的构建从邻接到距离衰减权重矩阵 W 是空间计量的核心输入它决定了“谁影响谁”以及“影响多大”。截面数据常用的 W 有三类类型构造方式适用场景注意事项邻接矩阵共享边界为 1否则 0行政区划数据岛屿需手动指定邻居距离阈值距离小于 d 为 1城市群、县域d 的选择影响结果距离衰减1/d 或 1/d²连续空间过程需设截断距离避免权重过大行标准化是标配操作让每行权重之和为 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, metriceuclidean) # 距离衰减权重截断距离取 100 公里根据坐标单位调整 threshold 1.0 # 若坐标已投影为公里此处为 100若为经纬度需先转换 W np.where((D 0) (D threshold), 1.0 / D, 0.0) # 行标准化 row_sums W.sum(axis1, keepdimsTrue) W_std np.divide(W, row_sums, outnp.zeros_like(W), whererow_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) * (yWy / yy)其中 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(fMorans 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 估计 SARfrom spreg import ML_Lag # y 和 X 准备好w 是行标准化后的权重矩阵 model_sar ML_Lag(y, X, ww, name_yy, 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, ww, name_yy, name_x[x1, x2]) print(fSEM 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, ww, name_yy, name_x[x1, x2, Wx1, Wx2]) print(fSDM 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, df2) # 两个 WX 项约束 print(fLR 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, df2) print(fLR test SDM vs SEM: {lr_sem:.2f}, p{p_sem:.4f})自由度取约束个数SDM 比 SAR 多两个 WX 系数所以 df2。如果 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_boot500): effects [] for _ in range(n_boot): idx choice(len(y), len(y), replaceTrue) # 注意重抽样后需重新构建 W 的子矩阵此处简化示意 model ML_Lag(y[idx], X[idx], wW[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(axis0), effects.std(axis0)这段代码是简化示意实际重抽样时权重矩阵的子矩阵需要重新行标准化否则行和不为 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”的时候有底气回答。希望帮到你。本文还有配套的精品资源点击获取
截面数据空间计量模型：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β uu λWu ε。杜宾模型SDM则是两者的综合既包含 Wy 也包含 WX形式为 y ρWy Xβ WXθ ε。选哪个模型不是拍脑袋决定的。常见做法是先跑 OLS然后做拉格朗日乘子检验LM 检验及其稳健版本。LM-lag 显著而 LM-error 不显著倾向 SAR反过来倾向 SEM两个都显著考虑 SDM 或者用稳健 LM 统计量再判断。这个流程在 Anselin 的经典文献里有详细推导我一般会把它作为模型选择的起点而不是终点。2.2 空间权重矩阵的构建从邻接到距离衰减权重矩阵 W 是空间计量的核心输入它决定了“谁影响谁”以及“影响多大”。截面数据常用的 W 有三类类型构造方式适用场景注意事项邻接矩阵共享边界为 1否则 0行政区划数据岛屿需手动指定邻居距离阈值距离小于 d 为 1城市群、县域d 的选择影响结果距离衰减1/d 或 1/d²连续空间过程需设截断距离避免权重过大行标准化是标配操作让每行权重之和为 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, metriceuclidean) # 距离衰减权重截断距离取 100 公里根据坐标单位调整 threshold 1.0 # 若坐标已投影为公里此处为 100若为经纬度需先转换 W np.where((D 0) (D threshold), 1.0 / D, 0.0) # 行标准化 row_sums W.sum(axis1, keepdimsTrue) W_std np.divide(W, row_sums, outnp.zeros_like(W), whererow_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) * (yWy / yy)其中 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(fMorans 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 估计 SARfrom spreg import ML_Lag # y 和 X 准备好w 是行标准化后的权重矩阵 model_sar ML_Lag(y, X, ww, name_yy, 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, ww, name_yy, name_x[x1, x2]) print(fSEM 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, ww, name_yy, name_x[x1, x2, Wx1, Wx2]) print(fSDM 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, df2) # 两个 WX 项约束 print(fLR 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, df2) print(fLR test SDM vs SEM: {lr_sem:.2f}, p{p_sem:.4f})自由度取约束个数SDM 比 SAR 多两个 WX 系数所以 df2。如果 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_boot500): effects [] for _ in range(n_boot): idx choice(len(y), len(y), replaceTrue) # 注意重抽样后需重新构建 W 的子矩阵此处简化示意 model ML_Lag(y[idx], X[idx], wW[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(axis0), effects.std(axis0)这段代码是简化示意实际重抽样时权重矩阵的子矩阵需要重新行标准化否则行和不为 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”的时候有底气回答。希望帮到你。本文还有配套的精品资源点击获取
