博文
一招搞定空间横截面数据回归模型
||
【导读】 |
利用shp2mat包生成空间权重矩阵，然后用spatial.cross输出所有回归结果。 |
一、使用空间计量经济模型的动机 |
一种方法从来都是为了解决问题，而不是故弄玄虚。所以，在准备用空间回归分析之前，我们得先问自己，为什么要用它？这有两大考虑：一是模型更可靠，二是为了识别空间效应。在一个空间样本集中，样本点之间是相互影响的，这种影响表现在数据上就是Y的空间自相关性，一般用莫兰指数来衡量。空间自相关来源有三：或是Y之间相互影响，或是毗邻的X影响本身的Y，或是模型中忽略的因素存在空间关联性。
根据这三种关联机制，建模的思路也很直接：如果是邻居的Yj影响自身的Yi（反过来Yi也会影响Yj），那就把邻居的Yj值平均后视为新的自变量LY，加到X中去再回归。好比浙江的GDP受到本身投入水平的影响，但也与周边的GDP产出水平有关，因此需要将其毗邻省份，如上海、江苏、安徽、江西、福建的GDP平均后作为新的自变量。每个省份都如此处理，就得到了一列新的变量。
如果是邻居的X影响本身的Y，类似以上做法，把邻居X平均后得到的变量LX加进原有的X再做回归。一般来说，有多少个X，就有多少个LX。
如果模型中应该考虑LY或（和）LX而你没有考虑，统计上看，就等于产生了遗漏变量偏差，因此模型估计是不可靠的。
第三种情况，模型中忽略的因素间存在空间关联性，这种效应将被误差项吸收，造成误差项相关。如果忽略的因素外生性很强，如环境变量或是外生冲击等，其不会造成有偏性或是一致性等问题，因此在大样本下问题不严重，但其会影响估计效率。处理方法是将误差项设定为空间自回归，等于将其分解，一部分为空间自相关部分，则剩下的那部分就是白噪声了。举个例子，假如在城市某区域政府突然要建一个垃圾站，则将会整体拉低那篇区域的房价，使其变动表现出空间自相关性。假如在房价影响模型中你没有考虑这个影响（当然你可以设置虚拟变量或是加入与垃圾站的距离来予以考虑），则其影响就归入误差项中，统计上表现出来就是残差值的莫兰指数显著。
除了模型估计上的考虑外，有时候还想看一下邻居到底对自身有没有影响？如果有，是正向还是负向（通常为正向），影响效应到底有多大？这就是空间效应的识别，这是普通回归模型所做不到的。
基于上述三种空间效应的组合，可以得到7种不同形式的空间回归模型。这些模型中，用得比较多的是空间滞后模型（SLM）、空间误差模型（SEM）和空间杜宾模型（SDM）。
二、直接效应和间接效应 |
如何利用空间回归识别空间效应？这个比较复杂，它不能简答地理解为空间回归系数，如ρ、λ等，需要将总效应分解为直接效应和间接效应，后者即为空间溢出效应。
从样本总体来看，我们看到的是所有观测单元某一个X的平均变动，所造成的Y的平均变动。考虑一个包含3个单元、只有一个X的简单情况。在普通回归模型设定里，单元间是彼此独立的，如果X1（即第一个单元的X，下同）、X2、X3都变动1个单位，Y1、Y2、Y3分别变动1、2、3个单位，则Y的平均变动为( 1+2+3)] / 3=2个单位。但在空间回归模型设定里，我们观测到的是Y可能变动2.5个单位。因为单元间存在空间互动，所以X1变动1个单位，除了造成自己的Y1变动1个单位外（由于存在反馈效应，实际上会大于1），还使得旁边的Y2和Y3分别变动0.1和0.2个单位；X2变动一个单位，除了造成Y2变动2个单位外，还造成旁边的Y1和Y3分别变动0.2和0.3个单位；X3变动一个单位，除了造成Y3变动3个单位外，还造成旁边的Y1和Y2分别变动0.3和0.4个单位。这样，我们在总体上看到的是：X平均变动1个单位，Y平均变动[(1+0.1+0.2)+(2+0.2+0.3)+(3+0.3+0.4)] / 3=2+0.5=2.5个单位。2.5个单位变动的总效应中，2个单位是直接效应，0.5个单位是间接效应。
当然，现实情况远比以上例子复杂。如X1影响Y1、Y2和Y3，后两者又反过来影响Y1。每个人既影响自己，也影响邻居，邻居反过来影响自己，同时又影响邻居的邻居。只要没有孤岛，这种影响就如波浪般绵延不绝。等到了一个均衡状态，空间溢出效应就相对稳定了。
至于分解的计算公式，以及如何判别间接效应的显著性，这里不予详述，反正软件会算出来，我们懂其原理即可，具体请参看Elhorst（2014，p20-25）。
三、空间横截面数据计量经济分析 |
从数据角度看，要做空间回归，需要的不过是在普通数据集上加上空间权重矩阵，然后用软件分析即可。这里关键是空间权重矩阵的生成，工具有GeoDa、Stata、R等，各有千秋，都基于存储地理坐标数据的shape文件。这里以中国31个大陆省份为例，用笔者编写的shp2mat命令生成31*31的空间权重矩阵。这个命令的使用请参看本公众号已经发表的文章《shp2mat：基于shape文件生成空间权重矩阵》。一言以蔽之，无论是生成0-1形式矩阵，还是距离倒数矩阵（包括距离阈值），用这个命令一行代码即可得到n*n矩阵，且在有“孤岛”的情况下自动将其最近单元设为邻居。
经验上讲，如果想顺利运行空间回归，首先观测数据性质要良好，即不要有缺失值、无穷值等；其次是空间权重矩阵最好为稀疏（即有大量0值），每个单元确保至少有1个邻居（即不能出现某行全为0的情况），行标准化之前为对称最好，但这不是必须的。如果你不知道何种方式矩阵为最佳，那么就选最原始的0-1邻接矩阵，并保持对称。
运行以下命令，将数据和准备用来生成权重矩阵的shape文件读入内存。这里数据为中国大陆31个省份的规模以上工业产值和固定资产、劳动力数等投入产出数据，shape文件存储31个省份的地理边界坐标。
代码1：读入属性和空间数据，并生成空间权重矩阵 |
# 需要预先(安装)加载readXL、rgdal、spdep、spatialreg和texreg等包，略 # 读入属性数据并将变量取对数 > mydata <- read_excel("省规工业数据2000.xlsx") > mydata <- transform(mydata, lnproduct=log(工业总产值), lnassets =log(固定资产), lnlabors =log(从业人员)) # 读入shape文件数据 > myshp <- readOGR("China_province31_shp\\province31.shp", encoding="utf-8") # 将shape文件单元调整为与属性数据表一致 > myshp <- myshp[match(mydata$省份, myshp$NAME),] # 生成空间权重矩阵(默认为0-1空间邻接形式),然后予以行标准化 > source("shp2mat.R") > myswm <- shp2mat(myshp) Warning message: #21号单元海南没有毗邻省份，给它配上距离最近的邻居，并给出警告 Unit 21 has no neighbor, find its nearest unit 20 as its neighbor. > myswm <- swm.norm(myswm) # 权重矩阵输出,可供Matlab、Stata等软件使用，输出的矩阵包含行名和列名 > write.csv(myswm, "p31_swm.csv") |
数据准备好后，就可以用spatialreg包中的各种命令进行空间回归分析了。
代码2：横截面数据空间回归分析 |
# 设定回归公式 > fm <- lnproduct ~ lnassets + lnlabors # 将n*n权重矩阵转成listw对象，需明确设定行标准化(W)以适应lmSLX函数,否则其会将截距也空间滞后 > mylistw <- mat2listw(myswm, style="W") # OLS回归 > ols <- lm(fm, mydata) # 空间滞后模型SAR > sar <- lagsarlm(fm, mydata, mylistw) > summary(sar) Estimate Std. Error z value Pr(>|z|) (Intercept) -2.452503 0.633602 -3.8707 0.0001085 lnassets 1.185742 0.154514 7.6740 1.665e-14 lnlabors 0.048116 0.143074 0.3363 0.7366427 Rho: 0.15619, LR test value: 8.3213, p-value: 0.0039183 # 空间误差模型SEM > sem <- errorsarlm(fm, mydata, mylistw) # 空间复合模型SAC > sac <- sacsarlm(fm, mydata, mylistw) # 自变量空间滞后模型SLX,仅空间滞后部分变量 > slx <- lmSLX(fm, mydata, mylistw, Durbin=~lnlabors) # 空间杜宾模型SDM,空间滞后所有变量 > sdm <- lagsarlm(fm, mydata, mylistw, Durbin=T) # 空间杜宾误差模型SDEM > sdem <- errorsarlm(fm, mydata, mylistw, Durbin=T) # 通用嵌套模型GNS > gns <- sacsarlm(fm, mydata, mylistw, Durbin=T) |
以上多个回归结果可以用textreg包中的相应函数表格化输出：
代码3：表格化输出多个回归结果 |
# 将OLS、SAR、SEM、SAC和SDM五个模型结果输出为屏幕文本表格 > screenreg(list(ols, sar, sem, sac, sdm), custom.model.names=c("OLS","SAR","SEM","SAC", "SDM")) |
用笔者自编的函数spatial.cross可一次性运行所有8个回归，并将结果输出为word表格。
代码4：表格化输出多个回归结果 |
# 所有模型一次性估计,并将全部结果写出到工作文件夹的spatialreg.doc表格 > source("spatial_cross.R") > results <- spatial.cross(fm, mydata, mylistw, Durbin=T) # 函数lmSLX输出的SLX类对象不能被texreg包函数识别,将其改为普通的lm对象 > class(results[["SLX"]]) <- "lm" # 将结果表格输出到工作文件夹中的spatialreg.doc文件 > htmlreg(results, file="spatialreg.doc", doctype=TRUE, digits=3) |
四、总效应分解 |
总效应分解比较复杂。以空间杜宾模型为例，运行以下代码，即可将总效应分解为直接效应和间接效应，并通过模拟抽样分布的方法进行统计显著性检验。
代码5：总效应分解 |
# 空间杜宾模型效应分解,无显著性判断 > ev <- eigenw(mylistw) > sdm2 <- lagsarlm(fm, mydata, mylistw, Durbin=T, control=list(pre_eig=ev)) > impacts(sdm2, listw=mylistw) Direct Indirect Total lnassets 1.1369237 0.3835558 1.52047956 lnlabors 0.1091415 -0.1714666 -0.06232517 # 空间杜宾模型效应分解,进行显著性判断,需要用coda包模拟抽样分布 > library(coda) > tr.matc <- trW(as(mylistw, "CsparseMatrix")) > slm.imp <- impacts(sdm2, tr=tr.matc, R=1000) > summary(slm.imp, zstats=TRUE, short=TRUE) ======================================================== Simulation results (asymptotic variance matrix): ======================================================== Simulated standard errors Direct Indirect Total lnassets 0.1659784 0.8736238 0.9580027 lnlabors 0.1480907 0.7537447 0.8076951 Simulated z-values: Direct Indirect Total lnassets 6.8833015 0.5064105 1.6543708 lnlabors 0.6994679 -0.2821460 -0.1350527 Simulated p-values: Direct Indirect Total lnassets 5.8482e-12 0.61257 0.098052 lnlabors 0.48426 0.77783 0.892570 |
参考文献：Elhorst J P. Spatial econometrics: from cross-sectional data to spatial panels[M]. Berlin: Springer, 2014.
数据和代码下载请关注微信公众号：“思达区域经济研究方法”，SDAR-workshop 扫码或长按，关注该微信号 |
版权申明 本号所有图文资料，除特别说明外，其版权归浙江工业大学王庆喜教授领衔的“思达工作室”所有。转发、引用请注明原始出处。 |
网络链接 1、科学网博客：http://blog.sciencenet.cn/u/workshopofsdar 2、经管之家（原人大经济论坛）账号：R语言区域经济 |
公众号功能 1、致力于打造中国R语言空间数据分析第壹号； 2、探讨R语言、ArcGIS、matlab、GeoDa、Stata、SPSS等软件在空间数据分析中的应用； 3、讨论产业集聚、空间溢出、区域创新网络、城市群发展等热门研究主题。 |
https://blog.sciencenet.cn/blog-3376208-1222645.html
上一篇：shp2mat：基于shape文件生成空间权重矩阵
