从数据角度看，要做空间回归，需要的不过是在普通数据集上加上空间权重矩阵，然后用软件分析即可。这里关键是空间权重矩阵的生成，工具有GeoDa、Stata、R等，各有千秋，都基于存储地理坐标数据的shape文件。这里以中国31个大陆省份为例，用笔者编写的shp2mat命令生成3131的空间权重矩阵。一言以蔽之，无论是生成0-1形式矩阵，还是距离倒数矩阵（包括距离阈值），用这个命令一行代码即可得到nn矩阵，且在有“孤岛”的情况下自动将其最近单元设为邻居。
经验上讲，如果想顺利运行空间回归，首先观测数据性质要良好，即不要有缺失值、无穷值等；其次是空间权重矩阵最好为稀疏（即有大量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)
