Stata:一文遍览Stata官方空间计量命令:sp系列命令

发布时间:2020-10-09 阅读 52

Stata 连享会   主页 || 视频 || 推文

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

课程详情 https://gitee.com/arlionn/Course   |   lianxh.cn

课程主页 https://gitee.com/arlionn/Course

作者:潘星宇 (清华大学)


目录


Stata15 版本新添加了空间计量的官方命令, 本文将结合一些操作实例带领大家系统地了解stata15的官方空间计量命令—— sp 家族命令。

1. 命令概览

Stata15 提供的空间计量命令主要基于空间自回归模型展开,主要有以下几类:

A. 空间数据准备

  • spshape2dtashapefile 文件转换为 Stata 格式
  • spset 将现有数据声明为空间数据
  • spbalance 使面板数据具有很强的平衡性
  • spcompress 压缩状态格式的 shapefile 文件

B. 空间权重矩阵的构建和空间滞后项的生成

  • spmatrix 创建、操作和导入/导出权重矩阵
  • spgenerate 生成空间滞后变量(wx)

C. 空间计量的回归命令

  • spregress 截面 SAR 模型
  • spivregress 工具变量 SAR 模型
  • spxtregress 面板 SAR 模型

D. 空间计量的检验命令

  • estat moran 基于moran指数的检验
  • spregress postestimation 截面 SAR 模型的检验命令
  • spivregress postestimation 工具变量 SAR 模型的检验命令
  • spxtregress postestimation 面板 SAR 模型的检验命令

2. 准备工作:空间数据准备与空间权重矩阵的构建

让我们从一个简单的例子开始,研究二氧化碳排放以及与贫困的关系。有关 二氧化碳排放 的官方统计数据可在线获取 ; 我们仅为此保留了 2014 年的数据(并删除了给出区域总数的行)以生成此文件 LA-CO2-2014.dta。 我们还将使用多维贫困指数(IMD),这是一种结合了英国人口普查的信息的贫困数据。英国的 shapefile 文件可以从这里下载(未包含北爱尔兰)。 Shapefile 通常有一个 .zip 压缩文件,其中包含几个文件组成部分。对于 Shapefile 文件的介绍可以参考 空间权重矩阵的构建

下面我们开始来进行数据的导入和声明工作。

/* 数据的导入和声明 */

*- 解压缩 Shapefile 文件
unzipfile "Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain.zip" // 解压缩文件内容

*- 数据的导入
spshape2dta Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain

*- 数据的声明
use Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain.dta , clear
encode lad16cd, generate(lacode)
spset lacode, modify replace

在上述过程中,Stata 通过 spshape2dta 命令读取 shapefile 数据并构造了两个 .dta 文件。 其中一个更大文件名为 ……_shp.dta,其中包含了边界和位置的详细信息。 另一个文件名 为 …….dta,包含了 shapefile 数据本身带有的行政区划编码 "lad16cd" 。我们打开不含 _shp 后缀的文件,将 "lad16cd" 编码转换为数值型lacode的数值变量,然后使用 spset 命令将其声明为空间数据。这一步非常重要,因为声明为空间数据后,我们就可以进一步构建空间权重矩阵乃至进行空间回归的操作了。

然后我们将二氧化碳排放数据和多维贫困数据按照代码进行空间匹配

/* 数据匹配 */
merge 1:1 lad16cd using "LA-CO2-2014.dta" // 匹配二氧化碳排放数据
keep if _merge==3
drop _merge
save "CO2-merged.dta", replace  //这一步我们就完成了给二氧化碳排放数据的编码

import excel "File_10_ID2015_Local_Authority_District_Summaries.xlsx", ///
       sheet("IMD") firstrow clear 
keep LocalAuthorityDistrictcode2 IMDAveragescore
rename LocalAuthorityDistrictcode2 lad16cd
merge 1:1 lad16cd using "CO2-merged.dta" //匹配多维贫困数据
keep if _merge==3
drop _merge

接下来我们来进行空间权重矩阵的构建。spmatrix 命令是一个 Stata 官方提供的比较好用的构建空间权重矩阵的命令。关于空间权重矩阵构建的其他知识可以参考空间权重矩阵的构建

我们有两个权重矩阵构建的选择:第一个是假设与接壤的辖区可以相互影响。这有时称为adjacency matrix 或者 contiguity matrix (邻接矩阵)。第二种选择是使相关性与反距离成比例,它使用shapefile变量_CX和_CY来计算局部权限中点之间的距离。我们将尝试两者并进行比较。

/* 构建空间权重矩阵 */
spmatrix create contiguity W, rook
spmatrix create idistance W2

然后我们就可以分别得到一个边相邻的邻接矩阵(命名为 W )和反距离矩阵(命名为 W2)。

3. 截面空间回归:回归与图示

准备工作完成后,我们就可以进行空间回归工作了。需要注意的是,stata15 官方命令只支持 SAR, 即 空间滞后模型 因此需要更多模型拓展的小伙伴还是需要其他外部命令的协助。

我们首先用最简单的 OLS 模型进行基准回归用以对比:

/* 基准回归 */
generate logdom = log(domestictotal)
generate logpop = log(population)
generate pc = logdom-logpop // 生成人均二氧化碳排放量的对数指标

regress pc IMDAveragescore  // 基准回归
predict regres, residuals      

多维贫困指数与人均二氧化碳排放量的系数为 -0.0082(95%置信区间为 -0.0066 到 -0.0097,p <0.001):贫困指数越高的地区的排放量越低。

在这里我们可以用残差的可视化图来反映空间模型的适用性问题:

/* 残差的可视化 */
scatter regres domest, msymbol(Oh) name(regres, replace)
grmap regres, clnumber(9) fcolor(BuYlRd) // 在地图上绘制残差
name(regresmap, replace)

图中淡黄色表示残差为零(完全拟合),深蓝色表示较大的负残差(模型高估数据)和较暗红色表示大的正残差(模型低估了数据)。从图中可以看出残差在空间上是聚类的。有一些红色和蓝色的连片的大区域(即高值-低值聚集的区域)。伦敦周边的蓝色地区表明,尽管一些当局的贫困程度相当高,而另一些贫困程度较低,但人均二氧化碳排放总体上低于拟合值,而在英格兰北部边缘则相反,人均二氧化碳排放总体要高于拟合值。

我们再采用 spregress 命令进行空间滞后模型的回归。若附加 gs2sls 选项,则执行广义两阶段最小二乘法,语法格式如下:

spregress depvar [indepvars] [if] [in], gs2sls [gs2sls_options] 

各个选项的含义如下:

  • dvarlag(spmatname) 选项为必要选项,表示要调用的空间权重矩阵
  • ivarlag(spmatname: varlist) 表示采用的工具变量
  • errorlag(spmatname) 表示误差项也具有空间相关性(可以理解为SEM模型的选项)

我们亦可附加 ml 选项以便执行 MLE 估计,具体语法如下:

spregress depvar [indepvars] [if] [in], ml [ml_options] 

具体估计结果如下:

*-空间 SAR 模型回归(相邻矩阵)
spregress pc IMDAveragescore, gs2sls dvarlag(W)

Spatial autoregressive model                   Number of obs     =        326
GS2SLS estimates                               Wald chi2(2)      =     154.82
                                               Prob > chi2       =     0.0000
                                               Pseudo R2         =     0.2928
								
-----------------------------------------------------------------------------
             pc |     Coef.   Std. Err.     z    P>|z|   [95% Conf. Interval]
----------------+------------------------------------------------------------
IMDAveragescore | -.0062574   .0008112   -7.71   0.000  -.0078472   -.0046675
          _cons |  .5901449   .0273667   21.56   0.000   .5365072    .6437826
----------------+------------------------------------------------------------
W							
----------------+------------------------------------------------------------
             pc |  .1716607   .0343417    5.00   0.000   .1043522    .2389692
-----------------------------------------------------------------------------
Wald test of spatial terms:          chi2(1) = 24.99      Prob > chi2 = 0.0000

在相邻空间权重矩阵下,多维贫困指数与人均二氧化碳排放量的相关系数为 -0.0063 ( 95% 置信区间 -0.0047至 -0.0078,z = -7.7,伪R2 = 0.29);同样,贫困指数越高的地区的排放量越低。其空间自相关系数为 0.1716607,说明人均二氧化碳排放与相邻地区的排放量呈现出正相关关系。

Stata
*-空间 SAR 模型回归(反距离矩阵)
spregress pc IMDAveragescore, gs2sls dvarlag(W)

Spatial autoregressive model                    Number of obs     =        326
GS2SLS estimates                                Wald chi2(2)      =     189.68
                                                Prob > chi2       =     0.0000
                                                Pseudo R2         =     0.3868
								
-------------------------------------------------------------------------------
             pc |      Coef.   Std. Err.      z    P>|z|   [95% Conf. Interval]
----------------+--------------------------------------------------------------
IMDAveragescore |  -.0084929   .0007202   -11.79   0.000  -.0099046   -.0070813
          _cons |   .8576374   .0248013    34.58   0.000   .8090278     .906247
----------------+--------------------------------------------------------------
W2              							
----------------+--------------------------------------------------------------
             pc |  -.3431917   .0439964    -7.80   0.000   -.429423   -.2569604
-------------------------------------------------------------------------------
Wald test of spatial terms:          chi2(1) = 60.85      Prob > chi2 = 0.0000

在反距离空间权重矩阵下,多维贫困指数与人均二氧化碳排放量的相关系数为 -0.0085 ( 95% 置信区间 -0.0071 至 - 0.0099,z = -11.7,伪 -R2 = 0.37 )。

*-空间SAR模型的残差可视化

// 相邻矩阵的残差可视化
spregress pc IMDAveragescore, gs2sls dvarlag(W)
predict spres, residuals
scatter spres domest, msymbol(Oh) name(spres, replace)
grmap spres, clnumber(9) fcolor(BuYlRd) name(spresmap, replace) title("SAR (adjacency) residuals")
twoway (scatter spres regres, ms(Oh) ///
		ytitle("SAR (adjacency) residuals") ///
		xtitle("Linear regression residuals")) ///
	   (function y=x, range(regres)), ///
	    legend(off) name(resvres, replace)
	   
// 反距离矩阵的残差可视化
spregress pc IMDAveragescore, gs2sls dvarlag(W2)
predict spres2, residuals
scatter spres2 domest, msymbol(Oh) name(spres2, replace)
grmap spres2, clnumber(9) fcolor(BuYlRd) ///
      name(spresmap2, replace) ///
	  title("SAR (inverse-distance) residuals")
twoway (scatter spres2 regres, ms(Oh) ///
	    ytitle("SAR (inverse-distance) residuals") ///
		xtitle("Linear regression residuals"))     ///
	   (function y=x, range(regres)), ///
	    legend(off) name(resvres2, replace)
相邻空间权重矩阵的残差分布
相邻空间权重矩阵的残差分布
反距离空间权重矩阵的残差分布
反距离空间权重矩阵的残差分布

哪种类型的相关矩阵是更合意的?伪 R2 表明反距离矩阵更适合数据,但答案实际上取决于我们对数据上下文的理解以及我们试图回答的问题。英国的地方行政区的规模大不相同,因此两个行政区可能会被另一个行政区分开,从而不连续,但他们在地理上非常接近,这一点特别是在城市地区非常明显。而面积非常大,人口却比较稀少的农村地区可能有相反的差异:地域上相连但中心之间的距离很远。从图中我们可以看到反距离矩阵更好地适用于那些邻近矩阵模型的红色和低估的大型农村地区,以及伦敦及周边地区。

4. 截面空间回归:检验命令

这一节主要介绍两个命令: estat impact以及 lrtest ,前者用于对空间效应的分解,后者为似然比检验,检验模型的选择。

效应分解:

  • 平均直接影响: 衡量每个观测单位自变量对于因变量的平均总影响
  • 平均间接影响: 衡量某个观测单位对所有其他观测单位的影响
*-效应分解
spregress pc IMDAveragescore, gs2sls dvarlag(W)
estat impact 

------------------------------------------------------------------------------
                |            Delta-Met
                |      dy/dx   Std. Err.     z    P>|z|   [95% Conf. Interval]
----------------+-------------------------------------------------------------
direct        
IMDAveragescore |  -.0062816   .0008096   -7.76   0.000  -.0078685   -.0046947
----------------+-------------------------------------------------------------
indirect      
IMDAveragescore |  -.0009863   .0002044   -4.83   0.000  -.0013868   -.0005857
----------------+-------------------------------------------------------------
total         
IMDAveragescore |  -.0072679    .000856   -8.49   0.000  -.0089457   -.0055901
------------------------------------------------------------------------------

可以发现,在相邻空间权重矩阵下。直接效应为 -0.0062816,即贫困指数每上升1个单位,直接人均二氧化碳排放下降0.63%;间接效应为 -0.0009863,表示本地人均二氧化碳排放量的变化引起的周边地区的二氧化碳排放量变化。

模型的选择:我们在模型中加入

. spregress pc IMDAveragescore, ml errorlag(W) dvarlag(W)
. estimates store normalreg
. spregress pc IMDAveragescore, ml dvarlag(W)
. lrtest normalreg

Likelihood-ratio test                                 LR chi2(1)  =     59.46
(Assumption: normalreg nested in .)                   Prob > chi2 =    0.0000

这里显示加入空间滞后误差项的模型在小于0.01%的水平下显著。因此我们可以认为残差项也存在着非常显著的空间依赖性。因此加入空间滞后误差项是更为合意的。

5. 结合 spgen 命令模型拓展

在更多的情境下,我们也要考虑自变量X也存在着空间依赖的可能性。但 stata 官方并未提供 sdm 的相关命令,这时我们可以采用 spgen 命令生成空间权重矩阵 W 与自变量X的空间滞后项,再采用 spregress 命令进行回归,从而可以得到考虑了自变量空间依赖性的模型结果。需要注意的是:这一方法缺乏对一致性和渐近分布的数理推导,因此在理论上不够严谨,因此使用的时候需要慎重选择。

spshape2dta Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain, replace
spmatrix create contiguity W, rook  \\ 读取 shapefile文件并生成空间权重矩阵
spgenerate Wtransporttotal = W*transporttotal

Variable	Obs	Mean	Std. Dev.	Min	Max
pc	326	0.545492	0.129761	0.024488	0.965792
IMDAverage~e	326	19.46362	8.004457	5.009	41.997
WIMDAverag~e	326	14.35601	7.139981	0	34.3326

spregress  pc IMDAveragescore WIMDAveragescore,gs2sls dvarlag(W) 

Spatial autoregressive model                    Number of obs     =        326
GS2SLS estimates                                Wald chi2(3)      =     141.17
                                                Prob > chi2      =   	0.0000
                                                Pseudo R2       =     	0.3017
------------------------------------------------------------------------------
              pc |     Coef.   Std. Err.     z    P>|z|   [95% Conf. Interval]
-----------------+------------------------------------------------------------
pc               | 			
 IMDAveragescore | -.0080134    .000819   -9.78   0.000  -.0096186   -.0064083
WIMDAveragescore |  .0028504   .0014383    1.98   0.048   .0000314    .0056695
           _cons |   .639564   .0229944   27.81   0.000   .5944958    .6846322
-----------------+------------------------------------------------------------
W                | 			
              pc |  .0488928   .0526389    0.93   0.353  -.0542776    .1520632
------------------------------------------------------------------------------
Wald test of spatial terms:          chi2(1) = 0.86       Prob > chi2 = 0.3530

在反距离空间权重矩阵下,多维贫困指数与人均二氧化碳排放量的相关系数为 -0.0085 (95% 置信区间 -0.0071 至 - 0.0099,z = -11.7,伪 -R2 = 0.37 )。而多维贫困指数也呈现出一定的空间依赖性。但这里我们的空间Wald检验没有通过,因此加入WX的模型在这里并不适用。

6. 面板模型的推广

掌握了上述操作后,空间面板模型的操作就非常简单了,这里只需要同时进行空间数据的声明 spset 以及面板数据的声明 xtset。具体的操作命令如下:

*-空间面板操作命令
copy http://www.stata-press.com/data/r15/homicide_1960_1990.dta .
copy http://www.stata-press.com/data/r15/homicide_1960_1990_shp.dta .
use homicide_1960_1990
xtset _ID year  //声明面板数据
spset           //声明空间数据
spmatrix create contiguity W if year == 1990     //以1990年数据为基准构建空间权重矩阵
spxtregress hrate ln_population ln_pdensity gini i.year, re dvarlag(W) //空间滞后模型回归

具体的参数解释也与空间截面的SAR模型类似,同样地,也可以通过 spgenerate 命令加入自变量的空间滞后项进行空间回归。

参考文献

  1. Spatial Analysis in Stata 15

  2. 空间权重矩阵的构建

  3. Stata15 空间计量手册

  • help spregress
  • help spgenerate
  • help spxtregress
  • help lrtest

相关课程

连享会-直播课 上线了!
http://lianxh.duanshu.com

免费公开课:


课程一览

支持回看,所有课程可以随时购买观看。

专题 嘉宾 直播/回看视频
最新专题 DSGE, 因果推断, 空间计量等
Stata数据清洗 游万海 直播, 2 小时,已上线
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]

Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。


关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。直播间 有很多视频课程,可以随时观看。
  • 连享会-主页知乎专栏,300+ 推文,实证分析不再抓狂。
  • 公众号推文分类: 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类,主流方法介绍一目了然:DID, RDD, IV, GMM, FE, Probit 等。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会主页  lianxh.cn
连享会主页 lianxh.cn

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会学习群-常见问题解答汇总:
https://gitee.com/arlionn/WD