数据和代码下载请点击“阅读原文” 提取码: t2c9
【导读】 |
利用spdep包中的poly2nb等函数构建空间邻接关系&＃xff0c;然后用moran.test函数计算莫兰指数&＃xff0c;可实现批量计算&＃xff0c;这是GeoDa等软件所不具备的优势。 |
这个学期在上研究生的区域经济研究方法课&＃xff0c;正好利用这个机会讲一下空间计量经济分析&＃xff1a;一是空间自相关检验&＃xff0c;二是空间横截面数据分析&＃xff0c;三是空间面板数据分析。先从空间自相关分析开始。
计算莫兰指数需要准备&＃xff1a;一列观测数据(属性数据)&＃xff0c;以及表达观测单元空间相邻或距离远近关系的权重数据(矩阵)&＃xff0c;关键是后者。
空间权重数据(矩阵)需要基于shape文件生成。这里要注意的是&＃xff0c;确保属性数据和权重数据对应的空间单元顺序一致。在R语言里&＃xff0c;其做法为&＃xff1a;首先用rgdal包中的readOGR函数读入shape文件为spatial对象&＃xff0c;其次用poly2nb函数基于spatial对象得到各单元间的邻接关系&＃xff0c;即nb对象&＃xff0c;最后用nb2listw函数将nb对象转成listw对象&＃xff0c;就可以输入相应函数进行空间自相关和空间回归分析了。在构建空间权重时&＃xff0c;还有其他函数&＃xff0c;如基于距离范围找邻居的dnearneigh函数、找k个最近邻居的knearneigh函数、计算邻居直线或大圆距离的nbdists函数等&＃xff0c;可参看spdep包说明&＃xff0c;这里采用最常用的邻接关系来获得空间权重。
这里以中国大陆31个省级单元2000-2010年的规模以上工业产出数据为例&＃xff0c;演示如何计算莫兰指数。先计算2000年规模以上工业总产值的空间自相关莫兰指数。
首先读入属性数据&＃xff1a;
代码1&＃xff1a;准备属性数据 |
# install needed packages if not > install.packages(c("spdep", "rgdal")) > library(spdep) > library(rgdal) # set your own working directory if not > setwd("E:\\Spatial_Econometrics_for_Graduates\\R") > mydata 中国省份规模以上工业生产数据_2000_2010.csv", stringsAsFactors &＃61; F) # view data > head(mydata) 省份年份工业总产值固定资产净值全部从业人员 1 北京 2000 2565.38 1517.20 113.13 2 天津 2000 2606.38 1561.39 120.19 3 河北 2000 3426.05 2285.01 269.75 4 山西 2000 1216.86 1473.69 183.56 5 内蒙古 2000 748.97 946.87 85.34 6 辽宁 2000 4249.46 3449.08 295.18 |
接着准备权重数据。这里需要注意三个细节&＃xff1a;一是原始的shape文件包含了港澳台&＃xff0c;需要先去掉&＃xff1b;二是shape文件中的省份排列顺序与属性数据不一致&＃xff0c;需要调整为与属性数据一致&＃xff1b;三是海南没有邻居&＃xff0c;这里指定广东为其邻居&＃xff0c;同时需要将二者的空间邻接关系调整为对称。
代码2&＃xff1a;准备空间权重数据 |
# generate spatial weights object > shp stringsAsFactors &＃61; F, encoding &＃61; "UTF-8") # remove Hongkong, Macau and Taiwan > shp 香港", "澳门", "台湾")) # make observation order in spatial data same as attribute data > row.names(shp) > shp 年份 &＃61;&＃61; 2000, "省份"],] # get nb object from polygon object, need a little while > mynb # look into mynb, find Hainan has no neighbour > mynb 1 region with no links: 海南 # make Guangdong and Hainan to be neighbours each other > number.hainan 海南") > number.guangdong 广东") > mynb[[number.hainan]] # make neighbour relation symmetric > mynb # see if nb object is really symmetric > is.symmetric.nb(mynb) [1] TRUE # transform nb object to nblistw object, row standardisation as default > mylistw |
万事大吉&＃xff01;最后是计算莫兰指数。需要注意的是&＃xff0c;莫兰指数的检验是基于标准化后的莫兰指数统计量的分布来进行的。在计算这一分布的方差时&＃xff0c;有两种做法&＃xff1a;一是基于随机排列模拟计算&＃xff0c;二是假设数据为正态分布根据理论公式推算。moran.test函数默认通过随机模拟来进行显著性检验。如其参数设置randomisation &＃61; FALSE&＃xff0c;则为正态分布假设。由于空间自相关一般为正相关&＃xff0c;因此默认的检验方向为“greater”。
代码3&＃xff1a;计算莫兰指数 |
# calculate Moran&＃39;s index under assumption of normality > x 年份 &＃61;&＃61; 2000, "工业总产值"] > moran.test(x, mylistw, randomisation &＃61; FALSE, alternative &＃61; "two.sided") Moran I statistic standard deviate &＃61; 1.6897, p-value &＃61; 0.09108 alternative hypothesis: two.sided sample estimates: Moran I statistic Expectation Variance 0.16690302 -0.03333333 0.01404318 |
为了解莫兰指数的检验机制&＃xff0c;我们下面来随机模拟一下。为减少随机误差&＃xff0c;模拟次数为9999次&＃xff0c;加上实际数据给出的一次&＃xff0c;共10000次。
代码4&＃xff1a;莫兰指数模拟 |
# Moran&＃39;s index random simulation > set.seed(12345) > morperm > morperm$res[10000] [1] 0.166903 # plot simulated Moran&＃39;s index > morp > md > plot(md, main&＃61;"Moran&＃39;s I Permutation Test", xlab &＃61; "Reference Distribution", xlim &＃61; c(-0.4, 0.55), ylim &＃61; c(0, 4), lwd &＃61; 2, col &＃61; 2) > hist(morp, freq &＃61; F, add &＃61; T) > abline(v &＃61; morperm$statistic, lwd &＃61; 2, col &＃61; 4) |
将属性数据标准化后&＃xff0c;利用lag.listw函数计算其空间滞后值&＃xff0c;可画出莫兰散点图。同时还可发现&＃xff0c;莫兰指数实际上等于空间滞后值对其本身值(需要事先标准化)的回归系数。
代码5&＃xff1a;莫兰散点图 |
# Moran&＃39;s Index scatter plot > production > moran.plot(production, mylistw) # Moran&＃39;s index equals standardized regression coefficient > wx > lm(wx ~ x) (Intercept) x 0.1027 0.1669 |
最后&＃xff0c;让我们批量计算2000-2010年全部年份的莫兰指数&＃xff0c;绘图表示&＃xff1a;
代码6&＃xff1a;批量计算2000-2010年的莫兰指数 |
# batch calculate Moran&＃39;s indexes for all years > moran.all 工业总产值, mydata$年份, moran.test, mylistw) > statistic > p.value > moran.all > plot(2000:2010, moran.all$statistic, type &＃61; "o", xlab &＃61; "year", ylab &＃61; "Moran&＃39;s I") > text(2000:2010, moran.all$statistic, round(moran.all$statistic, 3), cex &＃61; 0.8) |
数据下载请关注微信公众号&＃xff1a;“思达区域经济研究方法”&＃xff0c;SDAR-workshop 交流请加微信群&＃xff1a;“思达区域经济研究方法”&＃xff0c;需邀请加入&＃xff0c;请注明“单位名称&＃43;姓名” 扫码或长按&＃xff0c;关注该微信号
|
版权申明 本号所有图文资料&＃xff0c;除特别说明外&＃xff0c;其版权归浙江工业大学王庆喜领衔的“思达工作室”所有。转发、引用请注明原始出处。 |
网络链接 1、网易博客&＃xff1a;http://wqx1976.blog.163.com/ 2、人大经济论坛账号&＃xff1a;R语言区域经济 3、知乎账号&＃xff1a;sdar, https://www.zhihu.com/people/sdar-77 |
公众号功能 1、致力于打造中国R语言空间数据分析第壹号&＃xff1b; 2、探讨R语言、ArcGIS、matlab、GeoDa、Stata、SPSS等软件在空间数据分析中的应用&＃xff1b; 3、讨论产业集聚、空间溢出、区域创新网络、城市群发展等热门研究主题。 |
![[大整数乘法] java代码实现](https://img1.php1.cn/3cd4a/24c6f/9f3/0133bb25da242824.jpeg)
京公网安备 11010802041100号