做空间计量的朋友应该都有过这种经历拿到一份研究区域的shp文件打开Stata一脸懵不知道地理数据该怎么喂给模型跑去搜“莫兰指数怎么算”看到一堆老教程让你装spatgsa、写矩阵、手动标准化结果折腾两小时还在报错。其实Stata从15版本开始就把一整套官方空间计量命令集成了进来也就是标题里说的sp系列。这套命令最让我觉得省心的地方在于它把地理文件读取、空间权重矩阵构建、空间回归估计、残差空间自相关检验这整条链路全部统一到了同一个语法框架里不需要你再混搭五六个第三方程序也不用在Excel和Stata之间来回导权重矩阵了。这篇文章我不打算念手册就把这套sp系列命令里最核心的16条挨个讲清楚每条是干嘛的、什么时候用、语法长什么样、有哪些容易翻车的细节。最后再用一个完整例子把从shp文件到莫兰指数检验再到效应分解的流程串一遍。无论你是刚接触空间计量、被xsmle和spatgsa搞到头大的新手还是已经在用第三方命令但想迁到官方体系的老手这篇都适合当一份速查手册来用。1. 把shapefile变成Stata能懂的格式spshape2dta与spset怎么配合空间计量和普通计量最大的区别就是你的数据里多了一份“地理关系”。Stata官方sp系列解决这个问题的思路非常清晰先让你把GIS领域的shapefile转成Stata自己的格式然后在数据文件上打一个“空间标记”之后所有命令都基于这个标记来识别单元之间的邻居关系。这就好比你去房产中介办事得先出示房产证编号让人家确认你名下那套房的位置后面签合同、贷款、缴税都靠这个编号关联不需要每次重新解释一遍你的房子在哪。1.1 spshape2dta地理数据入库的第一步拿到一份city.shp之后第一步是用spshape2dta命令把它转换成Stata数据集。这个命令的写法非常简单spshape2dta city.shp, replace saving(city)运行之后它会生成两个文件city.dta和city_shp.dta。这两个文件的关系需要特别说明一下很多第一次用的人会在这里懵。city.dta是你的分析数据里面包含了shapefile属性表里的所有变量比如城市名称、GDP、人口、年份这些而city_shp.dta是空间几何信息的“附属档案”里面存着每个多边形的顶点坐标、边界线这些地理信息。两个文件靠一个ID变量关联在一起后续所有命令都会自动识别这个关联关系。转换完之后我建议立刻打开数据看一眼use city.dta, clear describe list in 1/5重点确认三件事第一ID变量是否存在且唯一这在后面构建权重矩阵时是命门第二属性表里的变量名是否跟你的研究变量对得上比如原始shp里如果字段是中文名Stata导入后很可能变成乱码或拼音需要提前重命名第三数据里有没有经纬度坐标变量有的shp带经纬度字段有的不带坐标信息在做距离权重矩阵时会用到。1.2 spset给数据打上“空间身份标记”spshape2dta只是完成了格式转换真正让Stata知道你这份数据是“空间数据”的是spset命令。用法也很直白use city.dta, clear spsetspset不带任何参数时作用是检查当前数据是否已经具备空间数据资格如果缺少必要信息它会明确告诉你缺什么。比如它可能会提示你需要设置坐标变量、需要指定ID变量、或者需要指定坐标系。第一次转换完shapefile之后直接跑spset通常它会自动识别一切然后输出一段摘要信息告诉你数据里有多少个空间单元、ID变量叫什么、坐标系是什么类型。这里有一个非常实用的细节如果你的shp文件没有自带的经纬度坐标字段或者你用的是手工整理的数据而不是标准shapefile可以用spset的modify选项手动指定坐标变量spset, modify xc(longitude) yc(latitude) spset, modify coordsys(latlong)coordsys(latlong)的意思是告诉Stata这两个变量是经纬度坐标不是平面投影坐标。这个设置至关重要因为它影响后续计算距离和识别邻接关系时使用的地球曲率算法。如果你把经纬度当成平面坐标处理算出来的邻接关系和距离会有不小的偏差尤其是研究区域跨越大范围的时候。1.3 spbalance面板数据的“体检报告”如果你的数据结构是面板每个空间单元对应多个年份那在跑面板空间模型之前先做一次“平衡性体检”非常有必要。spbalance命令就是干这个的spbalance输出结果会告诉你每个空间单元的观测次数是否一致、哪些单元缺了哪些年份。空间面板模型里非平衡面板虽然也能估计但估计结果的解释会被复杂化。我第一次跑spxtregress时报了一个难懂的错排查了半天才发现是某个城市缺了一年的数据导致空间权重矩阵和样本观测对不上。后来我养成了一个习惯任何面板数据跑回归之前必先xtset加spbalance双体检。2. 空间权重矩阵是灵魂用spmatrix构建“谁和谁挨着”的规则如果说sp系列命令里只能学一个那我一定选spmatrix。空间权重矩阵这东西说白了就是一张表格每个格子填写两个空间单元之间的“关联强度”。为什么说它是灵魂因为后面所有sp命令的回归结果都直接受这张表的影响——同一个模型、同一份数据换一种权重矩阵系数估计可能截然不同。这就像你定义“朋友关系”时按“住得近”、按“经常联系”、按“有什么共同爱好”划分得到的社交网络结构完全是两套后续分析自然不同。2.1 三种最常用的矩阵构建方式spmatrix create是构建权重矩阵的主命令最常用的有三种后缀* 邻接矩阵相邻为1否则为0 spmatrix create contiguity Wcon, replace * 距离倒数矩阵距离越近权重越大 spmatrix create idistance Wdist, replace * K最近邻矩阵每个单元取距离最近的K个邻居 spmatrix create knn Wknn, replace三种矩阵分别对应不同的实际场景。contiguity适合行政区划这种“接壤即关联”的情况比如研究省份间的经济溢出两省是否相邻是最直接的关系刻画idistance适合研究城市或企业这类点状对象比如研究城市间的知识溢出距离衰减是更合理的假设knn适合区域存在“孤岛”或边界不规则的情况强制给每个单元都分配固定数量的邻居能避免某些单元一个邻居都没有。构建完成之后一定要做两件事看摘要、画矩阵图。spmatrix summarize Wconspmatrix summarize会报告矩阵的维度、非零元素个数、密度等信息。我特别关注“是否有全零行”——如果某个空间单元没有任何邻居它的那一行权重就全是0这在后面回归时会导致该单元被静默排除。比如研究海南和大陆省份的邻接关系时海南如果不把广东设成“虚拟邻居”它就是一个标准的孤岛单元。这种问题在summarize输出里通常会以行求和的方式暴露出来一眼就能看到。2.2 自建矩阵的导入导出别让经济距离矩阵难倒你很多研究不使用纯地理的邻接或距离而是用经济距离、人口规模差异、贸易往来强度来构建矩阵这类自定义矩阵需要你自己准备好。spmatrix支持从外部文件导入先看怎么导出spmatrix export Wcon using myW.dta, replace导出后你会得到一个两两配对的长表数据包含三个变量单元i的ID、单元j的ID、对应权重值。自建矩阵时你在Excel里也按这个格式整理第一列源单元、第二列目标单元、第三列权重值保存成CSV后用import delimited读入Stata再执行spmatrix import Wcustom using myW.dta, replace这里有几个特别容易出问题的点我必须强调。第一ID变量的取值必须是1到N的连续整数不能出现跳号或重复否则导入时矩阵行列对不上报错信息还特别难懂第二矩阵必须是对称的无向关系或者至少方向逻辑自洽第三导入前确认数据行数和空间单元数一致。2.3 行标准化的真实含义在构建完矩阵后官方教程通常建议做行标准化但很多人不理解为什么只是机械地照做。行标准化的操作很简单spmatrix normalize Wcon, normalize(row) replace它的数学含义是把每一行的权重除以该行之和使得每行加起来等于1。这样一来空间滞后变量的含义就变得非常直观一个单元的空间滞后值等于所有邻居某个属性的加权平均权重是行标准化后的系数。这在模型解释时价值巨大因为空间自回归系数p就可以直接解读为“邻居加权平均每变化1个单位对自身因变量的影响”。如果不行标准化权重的绝对值量级会影响系数的解释甚至影响收敛。我的经验是除非你有特别坚实的理论理由否则宁可做行标准化也别指望回归程序替你自动处理。3. 16个sp命令全景图从截面SAR到面板SDM该用谁这一节是全文的核心速查表。很多人学sp命令时最大的困惑是命令太多记不住不知道什么时候该用哪个。我把常用的带sp前缀的官方命令按功能分成四组整理成一张总表然后逐个展开。分组命令功能定位典型使用场景数据准备spshape2dta把shapefile转为Stata格式刚拿到地理数据的初始化数据准备spset设定/检查空间数据属性每次分析前的身份确认数据准备spbalance检查面板空间数据是否平衡面板模型估计前的体检权重矩阵spmatrix创建、管理、导入导出权重矩阵构建W矩阵的主要命令权重矩阵spweight兼容早期格式的权重构建处理旧程序遗留数据权重矩阵spatom生成空间滞后变量手工构造解释变量时权重矩阵sppack对多变量批量生成空间组合需要一次处理很多解释变量时距离计算spdistance计算单元间距离并保存自定义距离权重的基础距离计算spdist2计算当前数据与另一数据集的距离匹配最近邻、构建外部对照截面模型spregress截面空间回归(SAR/SEM/SAC)横截面数据的核心估计命令截面模型spivregress面板/截面空间回归工具变量解释变量内生时需要截面模型sppois空间泊松回归因变量是计数数据时面板模型spxtregress面板空间回归面板数据最常用命令面板模型spxtivregress面板空间回归工具变量面板数据且内生系统模型spsur空间似无关回归多个方程误差相关时可视化spgraph空间数据作图辅助展示地理分布常配spmap使用3.1 漂亮的回归主命令spregress、spxtregress与它的兄弟们真正的重头戏是spregress。它针对截面数据估计三类空间模型空间滞后模型SAR、空间误差模型SEM、以及同时包含两者的空间杜宾/自回归组合模型。基本语法是spregress y x1 x2, dvarlag(W) mldvarlag(W)表示把因变量的空间滞后项放进模型等价于在方程右侧加了一个W*y这就是SAR模型。如果你认为空间依赖不在因变量而在误差项用errorlag(W)选项这就是SEM模型。如果你两个都想要可以同时写spregress y x1 x2, dvarlag(W) errorlag(W) ml这等于估计了一个同时含空间滞后和空间误差的广义空间模型也就是常说的SAC。估计方法建议用ml因为ml得到的结果可以使用后续的estat impact和estat moran这些后验工具。当数据结构变成面板时对应的命令是spxtregress。它的语法和spregress几乎一一对应只是需要先xtset声明面板结构xtset city_id year spxtregress y x1 x2, fe dvarlag(W)这里fe是固定效应re是随机效应。特别需要注意的是面板模型里空间权重矩阵不随时间变化也就是说假定空间单元之间的邻接关系在研究期间保持不变——这对大多数区域研究是合理的但在人口迁移频繁、行政区划调整的年代就需要谨慎了。spivregress和spxtivregress则是在上述基础上加入工具变量处理内生解释变量语法基本沿用了Stata的ivregress风格用endog()选项指定内生变量、ivarlag(W: x)指定空间滞后的解释变量。sppois处理计数因变量适用于专利授权数、犯罪事件数、疾病发病人数这类数据。spsur则是进阶玩家的选择它允许同时估计多个方程并且让方程之间的误差相关适合研究住房价格和租金、或者多个污染物排放这种存在系统关联的方程组。3.2 被低估的辅助命令spatom、sppack与距离计算spatom和sppack的名字看着陌生但实际使用率很高。spatom的用途是把某个变量的空间滞后直接生成一个新变量而不是在回归命令里现算。spatom x1, wmatrix(W) gname(lag_x1)生成的lag_x1就等于W乘以x1这个向量的结果即每个单元邻居的加权平均。这在做探索性数据分析时特别好用比如你想画一张“本地区GDP增速及其邻区加权增速”的散点图spatom一步就能完成。sppack则允许你对一组变量批量生成空间滞后避免了写十几个spatom的重复劳动。spdistance和spdist2解决的是“单元之间的距离”问题。前者计算当前空间数据集内部两两单元之间的距离并可以保存为矩阵后者计算当前数据集中每个观测到另一个数据集中每个观测的距离常用于匹配最近邻的对照组。如果你要自己构建“高速公路通行时间距离矩阵”或者“通勤圈矩阵”spdistance生成的距离文件就是你二次计算的底子。3.3 spgraph可视化让邻居关系看得见spgraph是官方提供的空间图形命令它可以把空间对象和变量值画在地图上让你直观地看到变量的空间分布模式。不过说实话在实际使用中我看到更多研究者还是习惯装第三方命令spmap来出图因为spmap拥有更丰富的配色、图例和多图层控制。spgraph完全可以和spmap并存——用官方命令管理数据和建模用spmap做最终的报告图这是目前比较顺手的工作流。4. 莫兰指数在这套框架里的真实位置回归前探索与回归后诊断现在聊回标题里的莫兰指数。很多初学者以为空间计量的标准流程是“先算莫兰指数再跑空间回归”这个理解不能算错但要分清楚这里的莫兰指数到底是对原始变量算的探索性指数还是对回归残差算的诊断性统计量。两者的作用和命令完全不同。4.1 探索性莫兰指数官方命令为什么不直接给你先明确一个很多人搜遍全网都想知道答案的问题Stata的sp系列官方命令里确实没有一条直接叫moran或morani的命令。为什么因为官方把“探索性空间数据分析”和“空间回归模型的残差诊断”作了严格区分。前者属于描述性分析通常用GeoDa、ArcGIS、QGIS、PySAL等专业空间分析软件来完成能顺带画出莫兰散点图和LISA聚类图后者才是sp系列框架内应该用estat moran做的模型诊断。如果你的数据必须全程在Stata里处理想对原始变量y计算全局莫兰指数我的建议是退而求其次用spmatrix export把权重矩阵导出来配合Mata或者直接用spatom生成邻居加权均值再用普通的相关分析去近似。所谓全局莫兰指数本质上就是在衡量“变量值”和“邻居加权值”这两个向量之间的相关程度。顺着这个思路你可以把莫兰指数理解为空间版的相关系数取值在-1到1之间正值说明高值和高值扎堆、低值和低值扎堆负值说明空间上呈现高低交替的格局。4.2 estat moran才是回归后诊断的标准动作绕开探索性的坑回到建模后的诊断。当你在spregress或spxtregress估计完之后标准做法是spregress y x1 x2, dvarlag(W) ml estat moran, errorestat moran, error检验的是模型残差里是否还残留空间自相关。如果检验结果不显著说明模型已经比较充分地吸收掉了数据中的空间依赖如果显著说明你的空间结构没有被完全捕捉可能需要换模型设定比如在SEM基础上再加上解释变量的空间滞后SDM或者调整权重矩阵。这里有一个特别容易误解的点我见过很多同学把estat moran, error的结果直接解读为“因变量的全局莫兰指数”这是不对的。原始变量的莫兰指数高恰恰说明应该用空间模型而残差的莫兰指数是用来检验空间模型是否够格的。一个显著为正的因变量莫兰指数加上一个不显著的残差莫兰指数才是理想的结果组合。4.3 有没有更全面的检验做一次LM诊断更稳妥除了莫兰指数spregress还自带了一套基于拉格朗日乘子的空间依赖诊断工具。标准模型估计完之后可以跑spregress y x1 x2 estat lmspregress不加任何空间项先估计一个普通OLS然后用estat lm输出LM-Lag、LM-Error、稳健版的Robust LM-Lag和Robust LM-Error。这组统计量能帮你判断空间依赖主要来自因变量的滞后还是误差项或者两者都有。我的习惯是拿到数据后先用spregress y x1 x2跑一个不带空间项的基准模型看estat lm的输出再决定模型设定往SAR、SEM还是SAC方向走。这个流程比直接上来就写dvarlag(W)要稳健得多也更容易在论文答辩里讲清楚自己的建模逻辑。5. 一条链跑通完整实例从shp文件清洗到效应分解前面把每条命令都拆开讲了这一节把它们串成一条流水线。我以一个虚拟的“城市空气质量与产业结构”为例走一遍完整的空间计量流程。假设手头有air_city.shp属性表里有城市名、PM2.5年均值、第二产业占比、人口密度、年份数据结构是2020和2021两年。5.1 阶段一数据初始化与权重矩阵* 转换shp文件 spshape2dta air_city.shp, replace saving(air) * 打开分析数据设置空间属性 use air.dta, clear spset * 生成面板结构标识 encode 城市名, gen(city_id) xtset city_id 年份 * 体检空间单元是否平衡 spbalance * 构建邻接权重矩阵并做行标准化 spmatrix create contiguity W, replace spmatrix normalize W, normalize(row) replace spmatrix summarize W这一段跑完数据已经具备了空间身份权重矩阵也就位了。spmatrix summarize的摘要里我确认没有全零行之后才会进入下一步。5.2 阶段二基础检验摸清空间依赖模式* 普通OLS基准模型 LM诊断 spregress pm25 第二产业占比 人口密度 estat lm根据estat lm的结果决定模型方向。假设Robust LM-Lag比Robust LM-Error更显著那优先考虑SAR模型如果两者都显著可能得用SAC。这个决定必须记录在案论文方法部分写“根据LM诊断结果选择空间滞后模型”就是有据可依的。5.3 阶段三正式回归、莫兰指数诊断、效应分解* 估计SAR模型 spregress pm25 第二产业占比 人口密度, dvarlag(W) ml * 残差空间自相关检验 estat moran, error * 直接效应、间接效应溢出效应、总效应 estat impactestat impact是整条链路里最容易被忽略却最有价值的一步。因为在空间滞后模型里回归系数本身不能直接解读为边际效应。原因在于某个城市的解释变量变化会直接影响该城市的因变量同时通过空间滞后项传导到邻居城市邻居城市的变化又会通过同样的机制反馈回来。这个循环迭代最终收敛的结果被estat impact分解成三类效应。直接效应度量的是对本市的平均影响含反馈回路间接效应度量的是对其他城市的平均溢出影响总效应是两者之和。实际写论文时“第二产业占比的间接效应显著为正”这类结论比笼统地说“系数显著”要高级得多也更能体现空间计量的价值。5.4 阶段四面板空间模型收尾因为数据有两年最终模型应该用面板版本xtset city_id 年份 spxtregress pm25 第二产业占比 人口密度, fe dvarlag(W) estat impact注意面板模型里如果你想做SDM空间杜宾模型也就是把解释变量的空间滞后也加进来语法是spxtregress pm25 第二产业占比 人口密度, fe dvarlag(W) ivarlag(W: 第二产业占比 人口密度)ivarlag允许你指定哪些解释变量要生成空间滞后。SDM有一个额外的好处即便真实的数据生成过程是SEMSDM也能给出无偏的系数估计所以很多应用研究直接上SDM作为保守选择。6. 实跑中容易翻车的地方ID顺序、权重对齐与面板平衡最后一个章节重点说我在实际使用sp系列命令时踩过的坑以及帮别人排查时反复遇到的几个问题。这些问题单独看都不大但每一个都可能导致回归无法运行或者结果莫名其妙。6.1 坑位一spset之后ID变量不唯一spmatrix create构建权重矩阵时矩阵的行列顺序严格对应spset里的ID变量排序。如果ID变量不是每个空间单元唯一后面几乎所有命令都会报错或者静默出错。我遇到过一个案例转换完shp后数据里_ID看着是1到100但一数发现有个别重复原因是原始shp的属性表里存在重复要素ID。排查方法很简单提前用duplicates检查。6.2 坑位二权重矩阵的行和样本对不上spset之后如果对数据进行了drop、keep或者merge空间单元集合发生变化权重矩阵却还停留在旧版本运行spregress时就会报类似“observations not found in matrix”的错误。解决办法只有一个任何影响样本范围的操作之后必须重新spmatrix create。我在项目管理里养成了固定习惯把“数据清洗”和“矩阵构建”分成两个严格隔离的模块矩阵构建永远在数据最终定稿后执行。6.3 坑位三estat impact在GSEM估计后不可用spregress支持两种估计方法ml和gsem。gsem在某些复杂模型比如因变量是非正态分布、或者带缺失值处理时会更快但代价是estat impact和estat moran这类后验命令在部分版本上可能不支持。所以如果模型设定允许我建议默认用ml只在遇到收敛问题或复杂面板结构时才考虑gsem会更好。6.4 坑位四面板模型里的空间单元缺失spxtregress要求每个时间点上出现的空间单元集合必须一致否则权重矩阵和样本不匹配。spbalance这个时候就派上用场了。如果spbalance发现确实有城市缺年份可以用tsfill补齐面板或者用fillin把缺失的单位-年份对补出来再考虑缺失值处理。千万不能带着非平衡面板直接跑spxtregress跑出来结果也许有但解释空间效应时会非常难受。6.5 坑位五经纬度坐标与投影坐标系混淆如果你的shp文件原本用的是投影坐标比如北京54、西安80、WGS84 UTM但spset时没有正确指定spmatrix create idistance算出来的距离就不是以公里/米为单位而是以度为单位。这会导致距离衰减参数的含义完全变掉。判断方法是在spset之后看一眼输出的坐标系信息如果不确定直接在spshape2dta之前用GIS软件把shp重新投影成经纬度好过在Stata里反复试。实跑空间计量真正拉开差距的不是会不会敲命令而是知不知道命令背后的数据逻辑。sp系列最大的优点就是把空间分析的门槛降了下来让研究者能把更多精力放在模型设定和结果解释上。我自己从第三方命令迁移过来之后最大的感受是调试时间少了至少一半因为所有环节的报错信息清晰多了。如果你还在用老一套的spatgsa手工流程真心建议花一个下午把官方sp系列跑通你会发现空间计量没有想象中那么麻烦。