用R进行空间关联性分析
「全局溢出」当一个区域的特征变化影响到所有区域的结果时,就会产生全局溢出效应。这甚至适用于区域本身,因为影响可以传递到邻居并返回到自己的区域(反馈)。具体来说,全球溢出效应影响到邻居、邻居到邻居、邻居到邻居等等。
「局部溢出」是指影响只落在附近或近邻的情况,在它们影响邻邻区域之前就消失了。
对应全局与局部溢出,存在全局与局部自相关检验。全局自相关检验指标主要有 moran'I 指数、Geary 指数 C 统计量以及 Getis-Ord global G 统计量;局部自相检验指标主要有局部moran指数、Local Getis-Ord Gi 和 Gi* 统计量。此外,也可以通过图示的方法来显示空间的依赖关系。
下面将对相关指标的计算以及可视化的方法进行演示(所用数据均可通过公共号后台回复「地图可视化R」获取)。
全局空间自相关检验
首先导入数据和矩阵,并进行适当转换
library(readxl)
library("spdep")
# 设置工作路径
setwd('E:\\空间计量专题\\R-空间计量')
# 导入经济变量数据
cdata <- read_xlsx("E:/空间计量专题/R-空间计量/cdata.xlsx")
# 导入自定义矩阵并做适当格式转换
w1 <- read_xlsx("w1.xlsx", col_names = FALSE)
w2 <- as.matrix(w1)
w2[1:5, 1:5]
# 转换格式并标准化
w <- mat2listw(w2, style="W")
导入shp文件,合并数据
library(rgdal)
# 导入shp文件
shpt <- readOGR("广东地级市.shp")
# 合并数据
cdatashpt <- merge(shpt, cdata, by = "city")
moran'I 指数
1.Moran’s I 统计量的计算
>>> moran(cdatashpt$gdp2017, listw=w, n=length(cdatashpt$gdp2017), S0=Szero(w))
$I
0.065546870128173
$K
2.99333822437931
I
Moran’s I 统计量, K
变量的峰度
当数据中出现孤岛或者缺失值时,可通过下面的子选项进行调整:
zero.policy
默认值为空,使用全局选项值;如果为真,则将0分配给没有邻居的区域的滞后值,如果为假,则分配为NA
NAOK
如果 TRUE,则 x 中的任何 NA 或 NaN 或 Inf 值都将传递给外部函数。如果 FALSE ,则存在 NA 或 NaN 或 Inf 值将被视为错误。
2.Monte-Carlo simulation of Moran's I
>>> # Monte-Carlo simulation of Moran's I
>>> set.seed(12345)
>>> moran.mc(cdatashpt$gdp2017, listw = w, nsim = 999, alternative = 'greater')
Monte-Carlo simulation of Moran I
data: cdatashpt$gdp2017
weights: w
number of simulations + 1: 1000
statistic = 0.065547, observed rank = 790, p-value = 0.21
alternative hypothesis: greater
3.moran散点图
>>> # Moran 散点图
>>> moran.plot(cdatashpt$gdp2017, w, zero.policy=NULL, spChk=NULL, labels=TRUE, xlab=NULL, ylab=NULL, quiet=NULL)
labels
给具有高影响度量值的点添加字符标签,如果设置为FALSE,则不会为具有较大影响的点绘制标签
4.Moran’s I 空间自相关检验
>>> # Moran’s I test for spatial autocorrelation
>>> moran.test(cdatashpt$gdp2017, w)
Moran I test under randomisation
data: cdatashpt$gdp2017
weights: w
Moran I statistic standard deviate = 0.76045, p-value = 0.2235
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.06554687 -0.05000000 0.02308712
Geary 检验
>>> geary.test(cdatashpt$gdp2017, listw=w, randomisation=TRUE, alternative="greater")
Geary C test under randomisation
data: cdatashpt$gdp2017
weights: w
Geary C statistic standard deviate = 1.2898, p-value = 0.09856
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
0.79614300 1.00000000 0.02498145
Geary’s C 检验
1.计算 Geary’s C 统计量
>>> geary(cdatashpt$gdp2017, listw=w, n=length(w), n1=length(w)-1, S0=Szero(w))
$C
0.0796142995607693
$K
0.427619746339901
2.Monte-Carlo simulation of Geary's C
>>> geary.mc(cdatashpt$gdp2017, listw=w, nsim=999, alternative="greater")
Monte-Carlo simulation of Geary C
data: cdatashpt$gdp2017
weights: w
number of simulations + 1: 1000
statistic = 0.79614, observed rank = 116, p-value = 0.116
alternative hypothesis: greater
3.模拟分布图
>>> set.seed(12345)
>>> gdpgeary <- geary.mc(cdatashpt$gdp2017, listw=w, nsim=999, alternative="greater")
>>> plot(gdpgeary, type='l', col='orange')
>>> gdpgeary.dens = density(gdpgeary$res)
>>> polygon(gdpgeary.dens, col="gray")
>>> abline(v=gdpgeary$statistic, col='orange',lwd=2)
Getis-Ord global G 检验
>>> globalG.test(cdatashpt$gdp2017, listw=w, alternative="greater")
Getis-Ord global G statistic
data: cdatashpt$gdp2017
weights: w
standard deviate = -1.1248, p-value = 0.8697
alternative hypothesis: greater
sample estimates:
Global G statistic Expectation Variance
4.948903e-02 5.000000e-02 2.063625e-07
空间相关图
>>> w.nb <- w$neighbours
>>> spcorrI = sp.correlogram(w.nb, cdatashpt$gdp2017, order = 2, method = "I", style = "W", randomisation = TRUE)
>>> spcorrI
Spatial correlogram for cdatashpt$gdp2017
method: Moran's I
estimate expectation variance standard deviate Pr(I) two sided
1 (21) 0.065547 -0.050000 0.023087 0.7605 0.4470
2 (21) -0.130212 -0.050000 0.017101 -0.6134 0.5396
>>> plot(spcorrI, main="Spatial correlogram of gdp2017")
局部空间自相关检验
局部moran检验
>>> localmoran(cdatashpt$gdp2017, listw=w, alternative = "greater")
A localmoran: 21 × 5 of type dbl
Ii E.Ii Var.I Z.Ii Pr(z > 0)
1 -0.01370494 -0.05 0.1146316 0.107200142 0.4573151
2 0.65060843 -0.05 0.4279122 1.071021124 0.1420800
3 -0.32033200 -0.05 0.4279122 -0.413256930 0.6602908
4 0.01710619 -0.05 0.4279122 0.102585331 0.4591460
5 0.05648861 -0.05 0.1146316 0.314522027 0.3765623
6 0.48390574 -0.05 0.1929517 1.215458873 0.1120956
7 -0.48582426 -0.05 0.1929517 -0.992172269 0.8394433
8 0.10421133 -0.05 0.1929517 0.351068574 0.3627685
9 -0.26267734 -0.05 0.1146316 -0.628158331 0.7350499
10 -0.44651083 -0.05 0.1929517 -0.902673594 0.8166504
11 -0.38816462 -0.05 0.2712719 -0.649270674 0.7419183
12 0.02013525 -0.05 0.1929517 0.159665854 0.4365722
13 -0.03877614 -0.05 0.1459596 0.029378254 0.4882815
14 0.41014395 -0.05 0.2712719 0.883469050 0.1884914
15 0.02782689 -0.05 0.8978331 0.082135687 0.4672694
16 0.11878306 -0.05 0.2712719 0.324060781 0.3729460
17 0.56506605 -0.05 0.2712719 1.180917005 0.1188178
18 0.47054422 -0.05 0.1929517 1.185040809 0.1180007
19 -0.05142908 -0.05 0.2712719 -0.002743819 0.5010946
20 0.06940159 -0.05 0.1929517 0.271822740 0.3928792
21 0.38968220 -0.05 0.1459596 1.150860065 0.1248949
Local Getis-Ord Gi 和 Gi*统计量
>>> localG(cdatashpt$gdp2017, listw=w)
[1] -0.05909956 1.09552796 -0.05815482 -0.09536069 -1.69368469 -1.94012537
[7] 0.69295805 0.34995778 0.96964875 -0.34538087 -0.15784709 0.51802708
[13] 0.20826225 -0.87609799 -0.19862319 -1.09970094 -1.06316301 -1.42853393
[19] 0.38484915 1.17379925 -1.43159536
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = cdatashpt$gdp2017, listw = w)
attr(,"class")
[1] "localG"
基于残差项的moran检验
>>> ols <- lm(gdp2017 ~ kj2017 + l2017 + ks2017 + pe2017 + inex2017 + new_inc2017 + pri_en2017 + high_stu2017, data = cdatashpt)
>>> lm.morantest(ols, listw = w, alternative = "two.sided")
Global Moran I for regression residuals
data:
model: lm(formula = gdp2017 ~ kj2017 + l2017 + ks2017 + pe2017 +
inex2017 + new_inc2017 + pri_en2017 + high_stu2017, data = cdatashpt)
weights: w
Moran I statistic standard deviate = 0.95884, p-value = 0.3376
alternative hypothesis: two.sided
sample estimates:
Observed Moran I Expectation Variance
0.001873225 -0.123283278 0.017037907
散点图的绘制
>>> moran.plot(ols$residuals, w)
局部moran检验
localmoran(ols$residuals, w)
A localmoran: 21 × 5 of type dbl
Ii E.Ii Var.I Z.Ii Pr(z > 0)
1 -0.0363243420 -0.05 0.1168494 0.040006912 0.48404381
2 0.8370573121 -0.05 0.4404798 1.336560700 0.09068304
3 -1.1222105809 -0.05 0.4404798 -1.615537694 0.94690285
4 0.2373930487 -0.05 0.4404798 0.433025295 0.33249820
5 -0.1380790378 -0.05 0.1168494 -0.257667335 0.60166817
6 -0.1081663716 -0.05 0.1977570 -0.130799494 0.55203304
7 0.1159985869 -0.05 0.1977570 0.373283232 0.35446883
8 0.7082098711 -0.05 0.1977570 1.704996632 0.04409753
9 0.0612709718 -0.05 0.1168494 0.325513260 0.37239632
10 0.4204181464 -0.05 0.1977570 1.057835549 0.14506521
11 0.3774028474 -0.05 0.2786646 0.809648516 0.20907111
12 -0.0868802622 -0.05 0.1977570 -0.082933137 0.53304765
13 0.1244368071 -0.05 0.1492124 0.451580986 0.32578544
14 0.1874982314 -0.05 0.2786646 0.449903626 0.32638997
15 -0.5517955762 -0.05 0.9259254 -0.521481405 0.69898427
16 -0.1211737778 -0.05 0.2786646 -0.134827702 0.55362595
17 -0.0030386117 -0.05 0.2786646 0.088961079 0.46455642
18 -0.4168336708 -0.05 0.1977570 -0.824903760 0.79528688
19 -0.0539797929 -0.05 0.2786646 -0.007539102 0.50300764
20 -0.3916838150 -0.05 0.1977570 -0.768348944 0.77886005
21 -0.0001822581 -0.05 0.1492124 0.128967879 0.44869153
其他相关检验类似,只需将变量替换为残差项即可,不再演示。此外,这里仅演示了各命令的常见用法,需要具体了解的话,可以点击「阅读原文」,这里给出了各命令的详细说明。
更多精彩内容,扫码关注微信公众号!
推荐阅读
-
我用 Python 从美食网站上抓取了 3032 份食谱并对其进行了分析,它闻起来真香!
-
使用 R 语言进行简单的主成分分析 (PCA)
-
数据可视化(VII):用 Pandas 对香港酒店数据进行高级分析,包括相关系数、协方差、数据离散化、数据透视表和其他精美的可视化。
-
原来python还可以这么玩】用python反向抓取网易云评论进行情感分析
-
用 R 语言实现随机前沿分析 SFA、数据包络分析 DEA、*弃置水文学 FDH 和 BOOTSTRAP 方法
-
系统评估--用 R 语言实现数据包络分析的 DEA (VII)
-
纯干货分享 | 研发效能提升——敏捷需求篇-而敏捷需求是提升效能的方式中不可或缺的模块之一。 云智慧的敏捷教练——Iris Xu近期在公司做了一场分享,主题为「敏捷需求挖掘和组织方法,交付更高业务价值的产品」。Iris具有丰富的团队敏捷转型实施经验,完成了企业多个团队从传统模式到敏捷转型的落地和实施,积淀了很多的经验。 这次分享主要包含以下2个部分: 第一部分是用户影响地图 第二部分是事件驱动的业务分析Event driven business analysis(以下简称EDBA) 用户影响地图,是一种从业务目标到产品需求映射的需求挖掘和组织的方法。 在软件开发过程中可能会遇到一些问题,比如大家使用不同的业务语言、技术语言,造成角色间的沟通阻碍,还会导致一些问题,比如需求误解、需求传递错误等;这会直接导致产品的功能需求和要实现的业务目标不是映射关系。 但在交付期间,研发人员必须要将这些需求实现交付,他们实则并不清楚这些功能需求产生的原因是什么、要解决客户的哪些痛点。研发人员往往只是拿到了解决方案,需要把它实现,但没有和业务侧一起去思考解决方案是否正确,能否真正的帮助客户解决问题。而用户影响地图通常是能够连接业务目标和产品功能的一种手段。 我们在每次迭代里加入的假设,也就是功能需求。首先把它先实现,再逐步去验证我们每一个小目标是否已经实现,再看下一个目标要是什么。那影响地图就是在这个过程中帮我们不断地去梳理目标和功能之间的关系。 我们在软件开发中可能存在的一些问题 针对这些问题,我们如何避免?先简单介绍做敏捷转型的常规思路: 先做团队级的敏捷,首先把产品、开发、测试人员,还有一些更后端的人员比如交互运维的同学放在一起,组成一个特训团队做交付。这个团队要包含交付过程中所涉及的所有角色。 接着业务敏捷要打通整个业务环节和研发侧的一个交付。上图中可以看到在敏捷中需求是分层管理的,第一层是业务需求,在这个层级是以用户目标和业务目标作为输入进行规划,同时需要去考虑客户的诉求。业务人员通过获取到的业务需求,进一步的和团队一起将其分解为产品需求。所以业务需求其实是我们真正去发布和运营的单元,它可以被独立发布到我们的生产环境上。我们的产品需求其实就是产品的具体功能,它是我们集成和测试的对象,也就是我们最终去部署到系统上的一个基本单元。产品需求再到了我们的开发团队,映射到迭代计划会上要把它分解为相应的技术任务,包括我们平时所说的比如一些前端的开发、后端的开发、测试都是相应的技术任务。所以业务敏捷要达到的目标是需要去持续顺畅高质量的交付业务价值。 将这几个点串起来,形成金字塔结构。最上层我们会把业务目标放在整个金字塔的塔尖。这个业务目标是通过用户的目标以及北极星指标确立的。确认业务目标后再去梳理相应的业务流程,最后生产。另外产品需求包含了操作流程和业务规则,具需求交付时间、工程时间以及我们的一些质量标准的要求。 谈到用户影响的地图,在敏捷江湖上其实有一个传说,大家都有一个说法叫做敏捷需求的“任督二脉”。用户影响地图其实就是任脉,在黑客马拉松上用过的用户故事地图其实叫督脉。所以说用户影响地图是在用户故事地图之前,先帮我们去梳理出我们要做哪些东西。当我们真正识别出我们要实现的业务活动之后,用户故事地图才去梳理我们整个的业务工作流,以及每个工作流节点下所要包含的具体功能和用户故事。所以说用户影响地图需要解决的问题,我们包括以下这些: 首先是范围蔓延,我们在整张地图上,功能和对应的业务目标是要去有一个映射的。这就避免了一些在我们比如有很多干系人参与的会议上,那大家都有不同想法些立场,会提出很多需求(正确以及错误的需求)。这个时候我们会依据目标去看这些需求是否真的是会影响我们的目标。 这里提到的错误需求,比如是利益相关的人提出的、客户认为产品应该有的、某个产品经理需求分析师认为可以有的....但是这些功能在用户影响地图中匹配不到对应目标的话,就需要降低优先级或弃掉。另外,通常我们去制定解决方案的时候,会考虑较完美的实现,导致解决方案括很多的功能。这个时候关键目标至关重要,会帮助我们梳理筛选、确定优先级。 看一下用户影响到地图概貌 总共分为一个三层的结构: 第一层why,你的业务目标哪个是最重要的,为什么?涉及到的角色有哪些? 第二层how ,怎样产生影响?影响用户角色什么样的行为? (不需要去列出所有的影响,基于业务目标) 第三层what,最关键的是在梳理需求时不需一次把所有细节想全,这通常团队中经常遇到的问题。 我们用这个例子来看一下 这是一个客服中心的影响地图,业务目标是 3个月内不增加客服人数的前提下能支持1.5倍的用户数。此业务目标设定是符合 smart 原则的,specific非常的具体,miserable 是可以衡量的,action reoriented是面向活动的, real list 也是很实际的。 量化的目标会指引我们接下来的行动,梳理一个业务目标,尽量去量化,比如 :我们通过打造一条什么样的流水线,能够提高整个部署的效率,时间是原来的 1/2 。这样才是一个能量化的有意义的目标。 回到这幅图, how 层级识别出来的内容,客服角色:想要对它施加的影响,把客户引导到论坛上,帮助客户更容易的跟踪问题,更快速的去定位问题。初级用户:方论坛上找到问题。高级用户:在论坛上回答问题。通过我们这些用户角色,进行活动,完成在不增加客户客服人数的前提下支持更多的用户数量。 最后一个层级,才是我们日常接触比较多的真正的功能的特性和需求,比如引导到客户到论坛上,其实这个产品就需要有一个常见问题的论坛的链接。这个层次需要我们团队进一步地在交付,在每个迭代之前做进一步的梳理,细化成相应的用户故事。 这个是云智慧团队中,自己做的影响地图的范例,可以看下整个的层级结构。序号表示优先级。 那我们用户影响地图可以总结为:
-
用 R-tmap 制作带指南针和比例尺的空间地图
-
通过对用户访问论坛习惯的调查结果进行分析,发现每个人都喜欢用这种方式浏览论坛
-
用 R 进行中介分析