你是不是也遇到过这种情况标准Logit模型跑得漂漂亮亮系数显著、边际效应也符合预期结果审稿人轻描淡写问一句这个决策存在明显的空间溢出你考虑空间依赖了吗你当场就愣住了。我最早做企业绿色创新决策时就是这种处境。专利申请行为在县域之间存在明显的相互影响相邻地区的企业会互相观察、跟随这种效应不放进模型里前面的回归结果几乎站不住脚。于是我开始系统研究Stata里怎么落地空间Logit模型这中间踩了不少坑也把主流方法都试了一遍。这篇文章就是完整实操记录适合正在写论文、需要处理二元因变量空间依赖的Stata用户。文章不会给你一个不存在的一键命令而是给出真正能跑完、能写进论文的完整路径。1. 二元因变量也会有空间依赖传统Logit模型的失效场景1.1 哪些场景最容易出现空间关联先别急着看代码想清楚一个问题为什么0/1因变量也需要空间模型很多人的第一反应是个体决策当然只取决于自身特征。但实际研究中这类决策天然带着空间属性。我举几个最常见的例子企业创新决策企业是否申请专利会受同地区或邻近地区竞争对手的影响。竞争对手拿到了核心技术专利你的企业会跟进布局否则在产业链上就被动了。贷款违约预测区域经济下行会同时推高一整片地区的违约率相邻县域的银行坏账往往同涨同跌。技术采纳行为农户是否采用新种子、医院是否采用新诊疗设备都存在典型的同伴效应隔壁县用了效果好这边就会跟着用。政策扩散一个城市出台了某类规制政策相邻城市往往在短期内跟进。这些场景的共同点是什么个体结果之间存在相互依存关系。用计量的话说潜在变量y*之间存在空间自相关或者误差项中存在空间结构。这种关联如果不处理直接套用标准Logit问题就来了。1.2 忽略空间依赖时估计量到底会偏多少传统Logit的核心假设是观测值条件独立即第i个观测的概率P(y_i 1|x_i)只跟自己的x_i有关不包含任何其他人的信息。这个假设在空间数据面前基本不成立。把模型写成潜在变量形式y*_i x_iβ ε_iy_i 1(y*_i 0)如果真实的DGP是空间滞后形式y* ρW y* Xβ ε那么化简之后y* (I − ρW)^(−1)(Xβ ε)你会发现第i个观测的y*不仅取决于自己的x还通过空间乘子(I − ρW)^(−1)取决于所有邻居的x。准确的概率应该是P(y_i 1|x) Λ[(I − ρW)^(−1)Xβ]_i其中Λ是Logistic累积分布函数。如果你忽略ρWy这一项直接估计标准Logit相当于把空间乘子丢进误差项。而这个误差项与X往往相关因为邻居的X会通过空间反馈影响本地的y*于是产生典型的遗漏变量偏误。偏误的大小取决于ρ的真实值和W的结构严重时符号都可能反转。还有一个更隐蔽的问题即使你不在乎系数偏误只想要相关性的解读标准Logit的标准误也是错的。空间数据的信息量比独立数据低误差项的正相关会低估标准误导致假阳性。所以我现在的习惯是只要数据带地理坐标或者区域编码先跑一个Morans I检验看看因变量或模型残差是否存在显著的空间自相关。如果Morans I显著就直接进空间模型流程。这一步很快但能帮你提前预判审稿人的问题。2. Stata没有现成空间Logit三条可行路线怎么选2.1 先破除一个误区没有splogit这个命令很多人在Stata里输入help spatial logit然后发现官方帮助文件里根本没有这个条目。这不是你安装有问题而是Stata官方确实没有提供直接的空间Logit估计命令。Stata 15之后内置的空间计量命令主要是spregress、spivregress、spml、spmatrix、spgenerate这一套。这套命令只覆盖线性回归框架对于0/1因变量的非线性空间模型官方属于空白。第三方命令里偶尔能看到一些空间Probit的尝试但成熟度参差不齐我没有找到可以放心推荐给论文使用的稳定实现。但这不意味着Stata做不了空间Logit。关键在于理解非线性空间模型的结构然后用间接方法逼近。2.2 路线A控制函数法2SRI— 最推荐的实操路径控制函数法是处理内生解释变量的经典方法在二元选择模型中被称为Two-Stage Residual Inclusion即2SRI。核心逻辑分两步第一阶段把空间滞后项Wy看作内生变量用工具变量对它做线性回归提取残差v_hat。第二阶段把Wy和v_hat一起放进Logit回归。v_hat捕捉了Wy中与误差项相关的部分把它控制住Wy的系数就能被解释为空间效应的估计。这个方法的优势很明显代码量小只需要reg、predict、logit三条命令对Logit/Probit的扩展很自然不改变连接函数灵活性高可以随时更换工具变量和权重矩阵有Wooldridge (2015)等文献背书审稿人认可度高。缺点是第二阶段Logit的标准误没有自动反映第一阶段估计的不确定性需要手动用bootstrap修正。这个细节后面专门讲。2.3 路线B线性概率模型LPM spregress — 快速稳健如果你不想绕2SRI的弯子可以直接用LPM近似。把0/1因变量当作连续变量跑spregress再用estat impact得到直接效应、间接效应和总效应。Stata的spregress支持最大似然和广义空间两阶段最小二乘两种估计方法后者对分布假设更稳健。LPM的硬伤是预测值可能超出0到1的范围而且异方差结构特殊但对0/1比例在20%到80%之间的因变量LPM的系数解释与Logit通常高度一致。我在论文里基本都是把LPM结果作为主回归的配套证据审稿人也接受这种做法。2.4 路线C跨软件验证 — 用R的spatialprobit包如果你的审稿人对空间模型的理解比较深要求真正的非线性空间点估计那Stata内置命令确实搞不定。我目前的方案是Stata做数据清洗和权重矩阵构建导出后用R的spatialprobit包估计SAR-Probit模型。spatialprobit包使用贝叶斯MCMC方法估计空间自回归Probit模型可以输出系数的后验分布、直接/间接效应处理空间滞后项的内生性也相对严谨。代价是跨软件操作麻烦MCMC的收敛判断需要一点经验。三条路线的定位差异非常明显路线实现难度内生性处理论文认可度主要限制控制函数法2SRI低能依赖工具变量中高标准误需Bootstrap修正Wy系数非严格ρLPM spregress低自动处理高稳健性首选线性近似预测值可能越界R spatialprobit中高能MCMC高跨软件操作参数调整繁琐我的建议是主回归用2SRI-Logit稳健性检验用LPMspregress结果做补充时间允许再用R包做交叉验证。这个组合基本能应对绝大多数审稿意见。3. 空间权重矩阵W的构建从这里开始动手3.1 数据准备与spset设置不管选哪条路线第一个实操步骤都是把数据声明为空间数据。Stata内置的sp命令从15.0开始可用不需要额外安装。假设你的数据里有每个观测的经纬度坐标声明方式如下use firm_innovation.dta, clear spset id longtitude latitude, coordsys(latlong)注意coordsys(latlong)这个选项。经纬度数据的坐标系统如果不声明后面的距离计算会乱套。我最初犯过一个错误没加这个选项直接idistance结果Stata按平面坐标处理角度坐标算出来的距离完全不对空间滞后项Wy的估计值错得离谱系数符号直接反了。排查了两天才找到问题。如果你的数据是区县多边形并且已经准备好了shapefile也可以基于区域ID声明use county_data.dta, clear spset county_id, shpfile(county.shp)声明完之后用spdescribe检查数据状态。如果显示data are mi(set) as spatial之类的内容说明设置成功。3.2 选择邻接权重还是逆距离权重空间权重矩阵W的定义是整个分析中最需要琢磨的部分。Stata的spmatrix create命令支持多种权重设定最常用的就两类。第一类是邻接权重* 邻接矩阵共享边界的区域权重为1否则为0 spmatrix create contiguity W_cont适合区县、省份这类多边形数据。逻辑是只有直接接壤的区域才相互影响。第二类是逆距离权重* 逆距离矩阵权重与距离成反比 spmatrix create idistance W_idist * 只考虑50公里范围内的邻居 spmatrix create idistance W_idist50, dband(50)适合企业、医院、地块这类点数据。逻辑是距离越近影响越大。选哪种取决于你的研究问题。如果研究的是政策扩散邻接矩阵更直观如果研究的是创新溢出距离衰减往往更符合现实。我的习惯是主分析用一个稳健性检验用另一个只要结论不一致就说明模型对权重矩阵敏感需要进一步考察。3.3 行标准化公认的必做步骤矩阵创建完成后一个高频坑是忘记行标准化。spmatrix normalize W_cont, row spmatrix summarize W_cont行标准化的作用是把每行元素除以该行总和让每行的权重之和等于1。这样做的意义在于spgenerate生成的空间滞后变量可以理解为邻居的平均值而不是一个量纲混乱的加权和空间自回归参数ρ的范围通常被约束在合理的区间内不同区域邻居数量不同时行标准化能避免人口密集区天然获得更高权重。spmatrix summarize输出矩阵的规模、非零元素数量、行和的最小最大值等信息。这个命令一定跑一次主要看两点有没有行和为零的观测非零元素比例是否低得异常。4. 控制函数法估计空间Logit完整Stata流程4.1 生成空间滞后项假设因变量是patent企业当年是否获得绿色专利0/1核心解释变量是R_D研发投入、size企业规模、subsidy政府补贴。先在数据里生成因变量和自变量的空间滞后项use firm_innovation.dta, clear spset id longtitude latitude, coordsys(latlong) spmatrix create idistance W spmatrix normalize W, row * 因变量的空间滞后邻居的平均专利状态 spgenerate Wy W * patent * 自变量的空间滞后作为Wy的工具变量 spgenerate W_RD W * R_D spgenerate W_size W * size spgenerate W_sub W * subsidyspgenerate的语法是spgenerate 新变量名 W * 原变量名其中W是之前创建好的权重矩阵名。生成的新变量就是空间滞后变量含义是空间权重矩阵作用后的邻居平均值。这里Wy代表邻居企业专利状态的平均水平W_RD代表邻居研发投入强度的平均水平。后面第一阶段回归的关键就是这些WX变量。4.2 第一阶段回归工具变量与残差提取控制函数法把Wy当作内生解释变量处理。问题随之而来用什么作为工具变量空间计量文献给出的标准答案是用自变量的空间滞后WX作为Wy的工具。这个思路来自Kelejian和Prucha (1998)的广义空间两阶段最小二乘法。逻辑是邻居的特征通过空间乘子会影响本地的y*但它们不会直接进入本地的二元选择方程只通过Wy传导。所以WX满足相关性和排他性两个条件。操作上第一阶段就是普通线性回归* 第一阶段Wy对X、WX回归 reg Wy R_D size subsidy W_RD W_size W_sub predict double v_hat, resid跑完立即做弱工具变量检验。最简单的方法是看WX变量的联合显著性test W_RD W_size W_subF统计量如果小于10说明工具变量太弱后面的估计可能不可靠。这种情况下可以考虑加入二阶空间滞后WWX即WW*X作为额外工具spgenerate WW_RD W * W_RD spgenerate WW_size W * W_size reg Wy R_D size subsidy W_RD W_size W_sub WW_RD WW_size实际研究中第一阶段的F值通常不会太低因为WX和Wy之间存在天然的相关性。但如果你用了dband(50)这种限制很强的权重矩阵F值可能会掉下来这时候就需要调整距离阈值或者改用邻接矩阵。4.3 第二阶段Logit加入Wy和残差第一阶段提取残差v_hat之后第二阶段直接做标准Logit* 第二阶段Logit回归加入Wy和v_hat logit patent R_D size subsidy Wy v_hat, vce(robust)结果中Wy系数的含义是控制住其他变量和内生性之后邻居专利状态对本地专利概率的影响。如果系数为正且显著说明绿色创新行为在空间上呈现正向溢出。v_hat的系数也值得关注。它有两层作用统计层面v_hat显著意味着Wy确实内生控制函数法是必要的v_hat不显著说明Wy的外生性无法拒绝此时普通Logit加入Wy也能成立。机制层面v_hat捕获了Wy中不可观测的扰动成分它显著说明还存在某些共同冲击没有被X捕捉到。不过这里有个非常关键的细节直接用logit得到的标准误偏低。因为v_hat本身是从第一阶段估计出来的第二阶段把它当作已知变量处理忽略了第一阶段的抽样误差。严谨的做法是用bootstrap重复整个两阶段过程。capture program drop boot_2sri program define boot_2sri, rclass preserve reg Wy R_D size subsidy W_RD W_size W_sub predict double vv, resid logit patent R_D size subsidy Wy vv matrix b e(b) return scalar rho b[1,4] return scalar beta_rd b[1,1] restore end bootstrap r(rho) r(beta_rd), reps(500) seed(123): boot_2sri注意bootstrap程序里要重新生成残差vv不能直接用外部的v_hat。preserve和restore保证每次抽样都重新估计第一阶段。这个bootstrap输出的置信区间才是可信的。4.4 与LPM结果互相印证两阶段Logit跑完我总是会再用spregress跑一遍LPM作为对照spregress patent R_D size subsidy, dvarlag(W) gs2sls estat impactestat impact会输出每个变量的直接效应、间接效应和总效应。直接效应反映本地解释变量对本地结局的影响间接效应代表本地解释变量变化通过空间传导对邻居结局的影响。这是一个典型的输出结构变量直接效应间接效应总效应R_D0.012**0.0040.016**size-0.003-0.001-0.004subsidy0.008*0.0030.011*LPM结果如果和2SRI-Logit的符号、显著性基本一致那这组实证结果就比较扎实了。如果不一致优先检查权重矩阵和工具变量不要急着下结论。5. 论文汇报系数、边际效应和直接/间接效应5.1 Logit系数不能直接解释概率变化回归表格里的系数是Logit尺度上的对数几率比审稿人不会满足于系数为正且显著。你需要给出经济意义上的量化解释。最简单的做法是marginsmargins, dydx(R_D size subsidy) post这个命令计算平均边际效应含义是平均而言R_D每增加一单位专利概率变化多少个百分点。但如果模型里包含Wy这个平均边际效应严格来说只是直接效应的近似因为空间滞后项的存在意味着任何解释变量的变化都会通过空间乘子产生连锁反应。5.2 规范汇报LPM的直接/间接/总效应在包含空间项的非线性模型里直接效应、间接效应、总效应的完整计算需要模拟整个空间乘子结构Stata的标准margins做不到这一点。所以我的处理办法是分层次汇报主回归用2SRI-Logit的系数和平均边际效应回答是否存在空间效应、方向如何效应量用LPMspregress的estat impact结果回答因为只有线性框架能直接给出直接/间接/总效应两种方法交叉印证比只报一种可信度高得多。论文里我一般放一张主表2SRI-Logit结果和一张附表LPM效应分解两边的核心结论必须一致。5.3 审稿人常问的三个空间问题整理完结果后我通常会预演审稿人会怎么追问问题一权重矩阵W怎么选的换一种W结论还成立吗回答正文用逆距离矩阵附表放邻接矩阵的结果两个矩阵下核心变量的系数和显著性保持一致。这就是权重矩阵敏感性分析。问题二Wy会不会内生回答采用控制函数法处理内生性第一阶段工具变量是自变量的空间滞后WXDWH检验的v_hat系数显著说明控制函数法适用。同时报告第一阶段F值大于10排除弱工具变量。问题三为什么用Logit不用Probit回答两者在本研究样本下结论一致。主回归用Logit是出于系数可解释性稳健性检验里用Probit重跑一遍。如果有R包MCMC结果也可以作为补充证据。5.4 完整的空间Logit汇报组合清单我现在的论文模板固定包含这五项内容基准Logit结果不包含空间项作为参照2SRI-Logit主回归结果含Wy和v_hatLPMspregress的效应分解表直接、间接、总效应权重矩阵敏感性分析不同W下的ρ和核心变量系数工具变量诊断第一阶段F值、DWH检验。这五项全部到位审稿人在空间维度上基本挑不出大毛病。6. 实测避坑收敛失败、无邻居个案与其他隐藏问题6.1 完美预测与收敛失败空间Logit最常见的报错是日志输出里连续弹出perfect predictions或not concave警告。出现原因通常是因变量分布极端不平衡比如只有3%的企业申请了绿色专利。此时Wy的加入让某个变量的线性组合能够完美区分0和1极大似然估计无法在有限参数空间内收敛。我的处理步骤是先看因变量分布如果0或1的比例低于5%考虑改用probitProbit的尾部更平滑收敛难度低于Logit减少自变量个数或者剔除组内方差几乎为零的变量改进权重矩阵使用连续距离权重替换离散邻接权重平滑程度更高如果以上都不行就把主回归改成LPMspregress2SRI-Logit降级为符号稳健性检验。6.2 无邻居个案处理邻接矩阵的一个常见问题有些区域周边没有邻居比如海岛、边远县域。行标准化之后这些行的权重之和为0spgenerate生成的空间滞后变量对无邻居观测就是缺失值。Stata在估计时默认删除缺失值样本量会悄悄减少。排查方法是在提取滞后变量后检查缺失值count if missing(Wy)如果缺失个数很少直接删掉并说明如果很多建议改用逆距离矩阵因为只要空间有其他点距离权重就不会为零。6.3 逆距离矩阵的衰减参数不能拍脑袋spmatrix create idistance W, dband(50)里的50公里阈值以及默认的1/d距离衰减函数都不应该凭空设定。我的经验是先用坐标数据计算样本点的空间分布画出不同距离阈值下Morans I值的变化趋势找到空间自相关最强的距离范围再据此设定dband。如果自相关在100公里内持续存在dband就设100如果到了50公里后Morans I迅速归零dband设50更合理。这个过程需要多跑几组权重矩阵做敏感性检验论文里也能交代清楚选择依据。6.4 2SRI标准误修正前面已经提过第二阶段直接使用logit报出的标准误没有考虑第一阶段估计的影响。如果你的论文只汇报了未经修正的标准误遇到懂行的审稿人大概率会被质疑。解决方法是使用bootstrap重复整个两阶段过程。之前的boot_2sri程序就是为此写的。500次重复一般够了最稳妥用1000次。如果bootstrap运行时间太长还有一个折中方案在2SRI-Logit回归中把v_hat的系数看作对外生性的检验主回归系数表里汇报的置信区间用bootstrap结果正文脚注里说明标准误修正方法。6.5 Wy的系数并不等于空间自回归参数ρ最后纠正一个容易踩的概念坑。2SRI-Logit回归结果里Wy的系数本质上是控制函数法框架下的简化式系数它不是结构模型中的空间自回归参数ρ。真正的ρ是y* ρWy* Xβ ε这个结构方程中的反馈参数需要通过模拟极大似然或贝叶斯MCMC才能准确估计。2SRI给出的Wy系数在方向上可以近似反映空间效应的正负但在数值上不能直接解读为空间反馈强度。所以论文里的表述要严谨。不要写空间自回归系数ρ显著为正而应该写控制函数法估计的空间效应系数显著为正表明邻居创新行为与本地创新概率呈正相关。真正的参数ρ留给R包MCMC估计来验证。我做空间Logit这几年最大的感悟是这个模型真正的门槛不在命令代码而在空间权重矩阵的设计和工具变量的有效性论证。你在W上花的心思越足后面所有结果的解释力就越强。如果只能给一条建议我想说——先把spmatrix summarize的输出看明白再开始跑回归这会帮你避开一半以上的坑。