企业官网建设流程全解析

你是不是也遇到过这种情况：标准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β + ε_i

y_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的标准误也是错的。空间数据的信息量比独立数据低，误差项的正相关会低估标准误，导致假阳性。

所以我现在的习惯是：只要数据带地理坐标或者区域编码，先跑一个Moran's I检验，看看因变量或模型残差是否存在显著的空间自相关。如果Moran's 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-Logit，稳健性检验用LPM+spregress结果做补充，时间允许再用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 * subsidy

spgenerate的语法是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_sub

F统计量如果小于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 impact

estat impact会输出每个变量的直接效应、间接效应和总效应。直接效应反映本地解释变量对本地结局的影响，间接效应代表本地解释变量变化通过空间传导对邻居结局的影响。

这是一个典型的输出结构：

LPM结果如果和2SRI-Logit的符号、显著性基本一致，那这组实证结果就比较扎实了。如果不一致，优先检查权重矩阵和工具变量，不要急着下结论。

5. 论文汇报：系数、边际效应和直接/间接效应

5.1 Logit系数不能直接解释概率变化

回归表格里的系数是Logit尺度上的对数几率比，审稿人不会满足于"系数为正且显著"。你需要给出经济意义上的量化解释。

最简单的做法是margins：

margins, dydx(R_D size subsidy) post

这个命令计算平均边际效应，含义是平均而言，R_D每增加一单位，专利概率变化多少个百分点。但如果模型里包含Wy，这个"平均边际效应"严格来说只是直接效应的近似，因为空间滞后项的存在意味着任何解释变量的变化都会通过空间乘子产生连锁反应。

5.2 规范汇报：LPM的直接/间接/总效应

在包含空间项的非线性模型里，直接效应、间接效应、总效应的完整计算需要模拟整个空间乘子结构，Stata的标准margins做不到这一点。所以我的处理办法是分层次汇报：

• 主回归用2SRI-Logit的系数和平均边际效应，回答"是否存在空间效应、方向如何"；
• 效应量用LPM+spregress的estat impact结果回答，因为只有线性框架能直接给出直接/间接/总效应；
• 两种方法交叉印证，比只报一种可信度高得多。

论文里我一般放一张主表（2SRI-Logit结果）和一张附表（LPM效应分解），两边的核心结论必须一致。

5.3 审稿人常问的三个空间问题

整理完结果后，我通常会预演审稿人会怎么追问：

问题一："权重矩阵W怎么选的，换一种W结论还成立吗？"

回答：正文用逆距离矩阵，附表放邻接矩阵的结果，两个矩阵下核心变量的系数和显著性保持一致。这就是权重矩阵敏感性分析。

问题二："Wy会不会内生？"

回答：采用控制函数法处理内生性，第一阶段工具变量是自变量的空间滞后WX，DWH检验的v_hat系数显著，说明控制函数法适用。同时报告第一阶段F值大于10，排除弱工具变量。

问题三："为什么用Logit不用Probit？"

回答：两者在本研究样本下结论一致。主回归用Logit是出于系数可解释性，稳健性检验里用Probit重跑一遍。如果有R包MCMC结果，也可以作为补充证据。

5.4 完整的空间Logit汇报组合清单

我现在的论文模板固定包含这五项内容：

1. 基准Logit结果（不包含空间项，作为参照）；
2. 2SRI-Logit主回归结果（含Wy和v_hat）；
3. LPM+spregress的效应分解表（直接、间接、总效应）；
4. 权重矩阵敏感性分析（不同W下的ρ和核心变量系数）；
5. 工具变量诊断（第一阶段F值、DWH检验）。

这五项全部到位，审稿人在空间维度上基本挑不出大毛病。

6. 实测避坑：收敛失败、无邻居个案与其他隐藏问题

6.1 "完美预测"与收敛失败

空间Logit最常见的报错是日志输出里连续弹出"perfect predictions"或"not concave"警告。

出现原因通常是：因变量分布极端不平衡，比如只有3%的企业申请了绿色专利。此时Wy的加入让某个变量的线性组合能够完美区分0和1，极大似然估计无法在有限参数空间内收敛。

我的处理步骤是：

• 先看因变量分布，如果0或1的比例低于5%，考虑改用probit，Probit的尾部更平滑，收敛难度低于Logit；
• 减少自变量个数，或者剔除组内方差几乎为零的变量；
• 改进权重矩阵，使用连续距离权重替换离散邻接权重，平滑程度更高；
• 如果以上都不行，就把主回归改成LPM+spregress，2SRI-Logit降级为符号稳健性检验。

6.2 无邻居个案处理

邻接矩阵的一个常见问题：有些区域周边没有邻居，比如海岛、边远县域。行标准化之后这些行的权重之和为0，spgenerate生成的空间滞后变量对无邻居观测就是缺失值。Stata在估计时默认删除缺失值，样本量会悄悄减少。

排查方法是在提取滞后变量后检查缺失值：

count if missing(Wy)

如果缺失个数很少，直接删掉并说明；如果很多，建议改用逆距离矩阵，因为只要空间有其他点，距离权重就不会为零。

6.3 逆距离矩阵的衰减参数不能拍脑袋

spmatrix create idistance W, dband(50)里的50公里阈值，以及默认的1/d距离衰减函数，都不应该凭空设定。

我的经验是：先用坐标数据计算样本点的空间分布，画出不同距离阈值下Moran's I值的变化趋势，找到空间自相关最强的距离范围，再据此设定dband。如果自相关在100公里内持续存在，dband就设100；如果到了50公里后Moran's 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的输出看明白，再开始跑回归，这会帮你避开一半以上的坑。