做空间计量这几年,我踩过最多的坑几乎都集中在“不知道用哪个命令、权重矩阵到底怎么构建”这些最基础的地方。很多人打开Stata,第一反应是去搜“stata下载”“stata安装包”,装完就急着跑回归,结果连莫兰指数都没算,就开始做空间杜宾模型,最后审稿人一问“你的空间依赖性是真实存在的吗”就直接卡住。今天这篇就把莫兰指数和Stata里16个sp系列空间计量命令一次性讲透,从数据准备到结果解读,给你一条能直接照做的路径。
先说结论:在Stata里做空间计量,你不需要学一堆杂七杂八的外部命令,核心就围绕官方sp系列命令展开,再配合几个经典的第三方spat命令做诊断。官方命令负责建权重矩阵、估计模型、分解效应,spat系列负责算莫兰指数和做空间依赖诊断。两者搭配起来,横截面数据也好、面板数据也好,基本都能覆盖日常研究的全部需求。
1. 莫兰指数:动手建模前必须先回答的问题
1.1 莫兰指数到底在检验什么
空间计量和普通计量最大的区别,就是承认“邻居会影响你”。一个地区的房价高,往往不是因为那个地区本身条件好,而是因为隔壁地区房价高、经济活跃,这种空间溢出效应在普通OLS里是被忽略的。莫兰指数(Moran's I)就是用来回答“这种空间依赖性到底存不存在”的统计量。
它的计算逻辑很直白。你先要有一个空间权重矩阵W,W的第i行第j列表示地区j是不是地区i的邻居,或者地区i和地区j之间的距离有多近。莫兰指数的公式是:
[ I = \frac{n}{\sum_{i}\sum_{j}w_{ij}} \cdot \frac{\sum_{i}\sum_{j}w_{ij}(x_i-\bar{x})(x_j-\bar{x})}{\sum_{i}(x_i-\bar{x})^2} ]
看着复杂,其实核心思想就是:把某个变量x的取值,和“邻居的x取值”做相关分析。如果I为正且显著,说明高值地区旁边也是高值,低值旁边也是低值,存在正向空间聚集;如果I为负,说明高值旁边是低值,空间分布存在排斥效应;如果I接近0,说明变量在空间上是随机分布的,做不做空间计量就值得商榷。
用生活场景打个比方:你在一个菜市场里观察摊位的价格,如果贵的摊位总是扎堆在一起,便宜的也总是聚在另一个角落,那这个菜市场的价格就有正向空间自相关。如果贵的和便宜的交替分布,那就是负相关。如果价格高低完全随机,那就没有空间自相关。莫兰指数做的正是这件事,只不过把“摊位”换成了“地区”。
1.2 全局莫兰指数的Stata实现步骤
先提醒一句:计算莫兰指数之前,你必须有空间权重矩阵。很多人忽略了这个前置步骤,直接拿变量跑spatgsa,结果报错“variable not found”或者干脆算出一个无法解释的数值。
在Stata里计算全局莫兰指数,我推荐直接用第三方命令spatgsa。使用方法如下:
ssc install spatwmat, replace ssc install spatgsa, replace先安装命令,然后准备一个包含地区ID、变量值和权重矩阵的数据。权重矩阵可以用spatwmat生成,它支持从shapefile文件读取邻接关系,也支持根据经纬度计算距离权重。
use "yourdata.dta", clear spatwmat using "yourdata_shape.dta", name(W) xcoord(x) ycoord(y) spatgsa var1 var2, weights(W) moran这里的var1、var2是你想检验的变量,weights(W)指定权重矩阵,moran选项输出莫兰指数。运行结果会给出Moran's I的值、期望值、标准差和z统计量。z值大于1.65(单侧5%显著性水平)基本可以认为存在显著的空间自相关。
如果你用的是官方sp矩阵,也可以从spreg回归后的estat moran得到残差的莫兰指数,这个后面第4节会细讲。
1.3 局部莫兰指数与LISA:找出“热点”在哪里
全局莫兰指数只能告诉你“整体上有没有空间聚集”,但它不会告诉你具体是哪个地区在带动这个聚集。这时候就要用局部莫兰指数(Local Moran's I),也就是LISA。
局部莫兰的计算逻辑类似于把全局莫兰拆解到每个地区,对每个地区i计算它和邻居之间的局部相关性。Stata里实现这个的思路比较简单,用spatlsa命令:
ssc install spatlsa, replace spatlsa var1, weights(W) moran这个命令不仅会输出每个地区的局部莫兰指数值,还能生成莫兰散点图。散点图分成四个象限:第一象限是高-高聚集区,也就是所谓的热点;第三象限是低-低聚集区,是冷点;第二象限是低-高区域,第四象限是高-低区域,这两种都是空间异质性比较明显的地方。搞研究时,如果审稿人问你“空间聚集主要发生在哪些区域”,你就可以把LISA图贴出来,配合显著性地图说明。
实操时要注意,spatlsa默认的显著性检验是基于随机置换的,所以跑的Bootstrap次数越多越稳定。如果数据量比较大,建议把permutation次数设成999以上,得到的结果才经得住推敲。
2. 为什么建议直接使用Stata官方sp系列命令
2.1 官方sp命令 vs 第三方spat命令
很多人一搜“Stata空间计量”,搜出来的全是spatwmat、spatreg这类第三方命令。这些命令确实老牌,在Stata 15之前几乎是唯一选择。但Stata 15之后出了一套官方空间计量命令,也就是sp开头的那一批,包括spset、spmatrix、spgenerate、spreg、spxtreg等。这套命令和第三方命令最大的区别是:它把空间权重矩阵、空间滞后变量、回归估计整合到了一个统一的框架里,数据管理更规范,输出也更友好。
我个人的建议是:能用官方命令做估计的,尽量用官方命令;第三方spat系列命令,主要用于做莫兰指数和空间依赖诊断这些官方命令覆盖得不够细的部分。两者不是对立关系,而是互补关系。
2.2 理解官方sp系列的操作主线
官方sp命令的设计逻辑其实就三步:定义空间数据、构建空间权重矩阵、估计空间模型。这三个步骤对应三组命令。
第一步是用spset命令告诉Stata“这份数据是空间数据”,需要指定唯一ID变量、坐标变量(经度纬度)或者链接外部shapefile。如果不先spset,后面的spmatrix和spreg根本没法用。
第二步是用spmatrix命令创建和管理空间权重矩阵。spmatrix create contiguity会基于邻接关系生成0-1权重矩阵,spmatrix create idistance会基于逆距离生成权重矩阵。创建后可以用spmatrix save保存,spmatrix dir查看有哪些矩阵,spmatrix export导出。
第三步是估计模型。横截面数据用spreg,面板数据用spxtreg,工具变量场景用spivreg。回归完之后还能用estat moran检验残差是否还有空间自相关,用estat impact分解直接效应和间接效应。
这套流程的好处是权重矩阵和估计命令高度集成,不容易出错。很多第三方命令之间依赖关系混乱,一个矩阵格式不匹配就会报各种莫名其妙的错误。
2.3 牢记这3个基础命令,其他都顺理成章
开始学习sp系列时,不要想着一次把16个命令全记住,先掌握三个核心命令即可。
第一个是spset,这是空间数据声明命令。第二个是spmatrix create,这是权重矩阵构建命令。第三个是spreg,这是模型估计命令。把这三个命令串起来,一个最简空间自回归模型就出来了:
use "data.dta", clear spset id, modify spmatrix create contiguity W spreg x1 x2 x3, dv(y) ml weights(W)运行结果会报告空间自回归系数rho,如果rho显著为正,说明存在正向空间溢出效应。学到这里,空间计量的核心操作就已经完成了。剩下的命令,更多是在这个基础上做扩展,比如处理面板数据、加入工具变量、做更复杂的权重矩阵等。
3. 16个sp系列命令逐一拆解
3.1 数据准备类:spset、spbalance、spgenerate
spset是进入空间计量世界的第一道门。这个命令的作用是声明数据的空间属性,也就是告诉Stata“哪个变量是唯一ID,哪个变量是经度,哪个变量是纬度”。
spset id, coords(xcoord ycoord) modify如果你有shapefile,可以写成spset id, modify replace,Stata会自动匹配坐标信息。spset执行后可以用spset summarize查看空间数据概况,比如有多少个唯一多边形、坐标变量的范围。
spbalance是用来处理空间面板数据平衡性的。面板数据经常存在某些年份缺值的情况,缺值会导致权重矩阵和模型估计出现问题。spbalance的作用就是检查并转换数据,让每个时空单元在时间维度上保持平衡。它有几个选项:fill表示填充缺失的观测,base()表示指定基准年份,generate()可以把是否平衡的信息存成新变量。
spgenerate是一个被低估的神器。它的功能是生成空间滞后变量。所谓空间滞后,就是邻居变量的加权平均。比如你想创建一个“邻居人均GDP”的变量,可以先创建权重矩阵W,然后:
spgenerate lag_gdp = W * gdp这里的spgenerate会把W矩阵和gdp变量做矩阵乘法,生成一个lag_gdp变量。这个变量在做空间滞后模型的稳健性检验、画莫兰散点图、甚至做空间异质性分析时都很有用。
3.2 权重矩阵类:spmatrix、spdistance
spmatrix是整个官方sp系列的地基。它负责创建、导入、导出、管理空间权重矩阵。常用的创建方式有:
* 邻接矩阵(根据共同边界或顶点相连) spmatrix create contiguity W1 * 逆距离矩阵(距离越近权重越大,可设定阈值) spmatrix create idistance W2, cutoff(100) power(1)邻接矩阵适合处理行政区域数据,比如省、市、县这类有明显边界的数据;逆距离矩阵适合处理像企业、银行网点这类没有行政边界但有经纬度坐标的数据。spmatrix还支持spmatrix import从外部文件导入自定义矩阵,比如从GeoDa导出的gal文件或者从ArcGIS导出的权重矩阵。这个功能非常实用,因为很多审稿人要求你用不同权重矩阵做稳健性检验,你就需要反复创建和切换矩阵。
spdistance是计算两两地区之间距离的命令。它会生成一个距离矩阵,你可以把它当作构建自定义权重矩阵的基础。通常情况下,如果你用spmatrix create idistance,就不需要单独用spdistance,但当你要做距离衰减阈值分析、或者要检验距离的某种非线性影响时,spdistance能给你更多灵活性。它生成的距离矩阵可以用spmatrix save保存,也可以导出成数据文件供其他软件使用。
3.3 横截面估计类:spreg、spregcs、spivreg
spreg是官方横截面空间回归命令,支持三种主要模型:空间自回归模型(SAR)、空间误差模型(SEM)、空间杜宾模型(SDM)。它的语法很直观:
* 空间自回归模型 spreg y x1 x2, dv(y) ml weights(W) * 空间误差模型 spreg y x1 x2, error(1) dv(y) ml weights(W) * 空间杜宾模型 spreg y x1 x2, dv(y) ml weights(W) durbin(x1 x2)选择哪个模型,主要看研究假设。如果你认为被解释变量存在空间溢出效应,比如一个地区的房价会影响邻近地区房价,用SAR;如果你认为误差项存在空间相关,可能是遗漏了某个空间相关的变量,用SEM;如果你认为解释变量的空间滞后也影响被解释变量,用SDM。常用的模型选择策略是先用spatdiag做LM检验,再结合经济学理论决定。
spregcs是spreg的扩展,它把几种常见空间模型纳入一个更灵活的框架,可以通过选项组合空间滞后和空间误差项,还能输出直接效应、间接效应和总效应的分解。这个命令特别适合做政策评估类研究,因为你不仅要关心某个变量的回归系数,还要关心它对本地区的影响(直接效应)和对邻近地区的影响(间接效应)。
spivreg是空间工具变量回归命令。当你的解释变量存在内生性时,普通spreg会得到有偏估计。spivreg允许你指定工具变量,同时处理空间滞后项和内生解释变量。语法类似ivregress,但多了一个weights(W)选项来指定空间权重矩阵。这个命令用得相对少,但一旦遇到内生性问题,它就是救场的存在。
3.4 面板估计类:spxtreg、spxtregs
spxtreg是面板数据的空间回归命令。语法和xtreg非常相似,区别在于多了权重矩阵和空间效应选项。它支持固定效应和随机效应模型:
* 固定效应空间自回归模型 spxtreg y x1 x2, fe dv(y) weights(W) * 随机效应空间自回归模型 spxtreg y x1 x2, re dv(y) weights(W)面板数据的核心优势是能控制个体异质性。固定效应模型主要消除不随时间变化的地区特征影响,随机效应模型则假设个体效应与解释变量不相关。选择fe还是re可以用传统的Hausman检验,但要注意,空间滞后项的存在可能影响检验统计量的渐进性质,建议配合理论判断。
spxtregs是spxtreg的扩展版本,它支持更复杂的空间效应设定,包括空间固定效应、时间固定效应以及空间误差项。如果你的数据时间跨度比较长,面板结构复杂,spxtregs会是更好的选择。比如同时控制地区和年份固定效应的空间杜宾模型,spxtregs就可以通过选项组合实现。
3.5 经典诊断类:spatwmat、spatgsa、spatlsa、spatdiag、spatcorr、spatreg
这六个命令来自第三方,但它们在空间计量诊断环节价值极高。spatwmat用于生成空间权重矩阵,支持从shapefile读取邻接关系,也支持基于坐标计算距离矩阵。虽然官方spmatrix功能更全,但spatwmat生成的矩阵格式是spat系列命令通用的,所以在做莫兰检验时我常常用spatwmat生成权重矩阵再喂给spatgsa。spatgsa是全局莫兰和Geary指数的计算命令,输出结果简洁清晰,是做空间自相关初步检验的首选。spatlsa用于局部莫兰指数计算和莫兰散点图绘制,能识别出高-高聚集、低-低聚集、高-低离群等不同类型。spatdiag是OLS回归后的空间依赖诊断命令,对OLS残差做LM检验和稳健LM检验,帮你判断该用SAR还是SEM。spatcorr可以计算空间相关图,在不同距离范围内考察空间自相关随距离变化的趋势。spatreg是第三方空间回归估计命令,虽然估计方法相对传统,但对于某些特殊模型设定仍有自己的优势,适合作为稳健性检验的补充。
这个夜间诊断组合拳的思路是:先用spatgsa确认存在空间依赖,再用spatlsa观察聚集位置,然后跑OLS,用spatdiag判断模型形式,最后再用spreg或spatreg正式估计。
4. 完整实操案例:从莫兰指数到空间回归模型
4.1 数据准备和权重矩阵构建
为了讲清楚操作路径,我用一个模拟的横截面案例来演示。假设有100个地区,每个地区有GDP、教育支出等变量,以及经纬度坐标。研究目的是检验GDP是否存在空间聚集,以及某解释变量对GDP的影响是否具有空间外溢。
第一步,先把数据读入Stata,并用spset声明空间属性。由于没有shapefile,我使用经纬度坐标来声明:
use "spatial_data.dta", clear spset id, coords(xcoord ycoord) modify这一步至关重要。如果不做spset,后续一切sp开头命令都会报错。spset执行后,用spdescribe查看数据状态,确认ID无重复、坐标变量有效。
第二步,创建空间权重矩阵。我选择基于距离的逆距离矩阵,因为100个地区分布比较分散,单纯用邻接矩阵会导致很多地区没有邻居,权重矩阵过于稀疏:
spmatrix create idistance W, cutoff(50) power(1)cutoff(50)表示只有距离在50单位以内的地区才算邻居,power(1)表示权重取距离的倒数,即距离越近权重越大。
4.2 全局莫兰指数检验与结果解读
直接用官方spmatrix创建的矩阵计算莫兰指数,一个快捷方法是先用spgenerate生成空间滞后变量,再用corr命令计算Moran's I。更标准的做法是使用spatgsa。两者各有利弊,spatgsa更正式,能给出期望值和标准化统计量。
* 方式一:用spatgsa spatwmat using "spatial_data.dta", name(W) xcoord(xcoord) ycoord(ycoord) spatgsa gdp, weights(W) moran结果可以看到Moran's I约为0.35,z值约为6.8,p值小于0.001,说明GDP在空间上存在显著的正向聚集效应。这意味着GDP高的地区倾向于和GDP高的地区相邻,空间计量建模的合理性得到了初步验证。
如果你想使用官方spmatrix创建的权重矩阵来完成检验,可以通过:
spgenerate Wgdp = W * gdp corr gdp Wgdp这个相关系数虽然不等于莫兰指数,但也能反映空间滞后的相关程度。要获得比较严格的检验结果,还是建议使用spatgsa。
4.3 从OLS残差诊断到模型选择
在莫兰指数显著之后,下一步不是直接套模型,而是先做一个普通OLS回归,然后对残差做空间依赖诊断。这一步的作用是判断应该采用SAR、SEM还是SDM。
reg gdp edu invest predict e, resid spatdiag, weights(W)spatdiag会输出多种检验统计量,包括Moran's I、LM-error、LM-lag以及对应的稳健版本。如果LM-lag显著而LM-error不显著,优先考虑空间滞后模型SAR;如果LM-error显著而LM-lag不显著,优先考虑空间误差模型SEM;如果两者都显著,再看稳健版本,稳健lm-lag显著则选SAR,稳健lm-error显著则选SEM。如果两个稳健统计量都显著,就要考虑SDM或者更复杂的SAC模型。
本例中,LM-lag和稳健LM-lag均显著,LM-error不显著,因此选择空间自回归模型SAR是合理的。
4.4 估计空间回归模型并解读直接与间接效应
选定SAR模型后,使用官方spreg命令估计:
spreg gdp edu invest, dv(gdp) ml weights(W) estat impact, nose结果中的rho为空间自回归系数,显著为正,说明GDP存在正向的空间溢出效应。edu和invest的回归系数是直接的边际效应,但空间模型里解释变量对被解释变量的总影响不仅限于本地区,还有对邻居地区的间接影响。estat impact会输出直接效应、间接效应和总效应。
假设edu的直接效应是0.42,间接效应是0.18,总效应是0.60。你可以在论文里这样描述:教育支出每提高1%,不仅会显著促进本地区GDP增长0.42%,还会通过空间溢出效应促进邻近地区GDP增长0.18%。这就是空间计量区别于普通回归的核心价值。
如果你的数据能支撑面板模型,可以把流程换成spset配合spxtreg,其他诊断思路完全一致。面板数据还能控制地区固定效应,缓解遗漏变量偏误。
5. 常见问题与排查技巧
5.1 shapefile导入报错:spset无法识别坐标
很多人拿到一份shapefile,直接在Stata里use,然后spset,结果报错“variable not found”。原因很可能是Stata没有直接读取shapefile,你需要先把shapefile转换为Stata格式。办法是用shp2dta命令:
shp2dta using "province.shp", data("province_data.dta") coords("province_coords.dta") gencentroids(centroids)这个命令会把shapefile的属性数据和坐标数据分开导出。然后再用spset去连接。转换时要注意shapefile的坐标系,如果是经纬度,spset里不用额外设置;如果投影坐标,建议先处理坐标单位,保证距离计算符合实际。
5.2 权重矩阵行标准化警告
在做空间计量时,Stata经常会提示“row-standardization recommended”,提醒你把权重矩阵标准化。这是一个关键问题。行标准化的目的是让每行权重之和为1,这样空间滞后变量的含义就是“邻居变量的加权平均”,解释起来更直观。
如果你用的是spmatrix创建的矩阵,可以在创建后使用spmatrix normalize W, row进行行标准化。如果你用的是spatwmat,可以给spatgsa或spatdiag指定标准化选项,比如standardize。忽视这个步骤可能导致莫兰指数计算出来的值偏大或偏小,显著性检验也会失真。
5.3 面板数据spxtreg报错:面板ID变量不匹配
spxtreg对数据结构很挑剔。它的要求是:经spset声明的空间ID必须和面板ID完全对应,且不能在时间维度上有缺失。常见的报错是“panel data must be strongly balanced”或者“spatial panel data are not balanced”。解决办法是先用xtset设置面板结构,然后用xtdescribe查看是否有缺失年份,再用spbalance调整。
还有一些情况,你的地区ID和形状数据里的ID变量名字不同,导致spset时无法自动匹配。解决办法是在spset前先看一眼shapefile数据的ID变量,比如用describe查看,再用rename改成一致的名字,或者用spset的id选项显式指定。
5.4 局部莫兰指数的结果无法保存
spatlsa输出的是每个地区的局部莫兰值、p值和散点图,但它不像普通回归那样自动把所有结果存入e()或r()。你如果要导出LISA表格,需要手动把结果保存下来。一个技巧是使用spatlsa的graph选项生成莫兰散点图,然后用Stata的graph save保存图片。如果要导出数值,可以在spatlsa运行的临时数据区里手动复制,或者用log记录来截取结果。
另一个更高效的方式是直接用spgenerate计算空间滞后变量,然后用常见的统计命令画散点图:
spgenerate Wgdp = W * gdp scatter Wgdp gdp, mlabel(id)这个散点图本质上就是莫兰散点图,横轴是本地区值,纵轴是邻居平均值。斜率为正且越陡,空间正向相关越强。这个方法的好处是可以用到官方命令生成的权重矩阵,不需要来回切换工具。
5.5 一个容易忽视的细节:莫兰指数对权重矩阵高度敏感
最后分享一个容易被忽视的细节:莫兰指数的结果对权重矩阵的选择非常敏感。同一份数据,用邻接矩阵算出来可能显著,用逆距离矩阵算出来可能不显著。这不是软件问题,而是空间计量的固有特性。所以在论文里,一定要报告你所用的权重矩阵类型和构建方式,并且建议做多种权重矩阵的稳健性检验,让结论不那么依赖单一阵设定。我在实际研究中,通常会同时使用邻接矩阵和不同的距离阈值逆距离矩阵,如果莫兰指数和空间回归系数在多种设定下都保持一致的方向和显著性,结果才算是真正站得住脚。如果结果在不同矩阵下差异很大,就需要仔细分析是什么原因造成的,是矩阵太稀疏,还是距离阈值选得不合理,而不是直接挑一个“好看”的结果发出去。这个习惯帮我避开了不少审稿质疑。