利用shp2mat包生成空间权重矩阵，然后用spatial.cross输出所有回归结果。

# 需要预先(安装)加载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")

# 设定回归公式

> 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)

# 将OLS、SAR、SEM、SAC和SDM五个模型结果输出为屏幕文本表格

> screenreg(list(ols, sar, sem, sac, sdm),

                custom.model.names=c("OLS","SAR","SEM","SAC",   "SDM"))

# 所有模型一次性估计,并将全部结果写出到工作文件夹的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)

# 空间杜宾模型效应分解,无显著性判断

> 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

版权申明

本号所有图文资料，除特别说明外，其版权归浙江工业大学王庆喜教授领衔的“思达工作室”所有。转发、引用请注明原始出处。

网络链接

1、科学网博客：http://blog.sciencenet.cn/u/workshopofsdar

2、经管之家（原人大经济论坛）账号：R语言区域经济

公众号功能

1、致力于打造中国R语言空间数据分析第壹号；

2、探讨R语言、ArcGIS、matlab、GeoDa、Stata、SPSS等软件在空间数据分析中的应用；

3、讨论产业集聚、空间溢出、区域创新网络、城市群发展等热门研究主题。