121 人关注0 条评论
写回答
查看全部 5 个回答
写回答
叶山Shan Ye
GIS/地质/人文地理/可持续发展
A2A 谢邀，
我和我认识的一些人，刚开始用R做空间分析的时候，也遇到过这个问题。R这种开源的东西，优点是各种包很丰富，缺点是有些包的说明写得很乱，地理加权回归（GWR）的R包其实功能很强大，但大部分说明都不大靠谱。
GWR在R里面可以用好几个不同的包来实现，其中步骤最简单的是spgwr。思路就两步：建立窗口、用窗口扫全局。这其实就是GWR本质上的两步。比如我要在全美国范围内统计某两个（或多个）变量之间的回归关系，我可以做一个全局回归（global regression），但因为这些变量在空间分布上或许会有异质性（heterogeneity），表现在统计结果上就是空间不稳定性（nonstationarity），因此只看全局的统计，可能看不出什么结果来。举个不完全恰当但是很容易领会精神的例子，你比如说，我要分析亚洲范围内，经济发展程度与牛肉销量之间的关系，经济越发达的地方，人们就越吃得起牛肉。可是等我统计到印度的时候，坏了，印度大部分人不吃牛肉，这不是经济状况导致的，这一下就影响了全局统计的参考价值，那怎么办呢？我们可以建立一个窗口（正规说法是带宽窗口，bandwidth window），每次只统计窗口范围内的经济与牛肉销量的关系，然后用这个窗口去扫过全局的范围。等统计到印度的时候，印度内部的各地和印度自己比，吃牛肉的人的比例就不会突然减少，这样就能减少这种空间不稳定性对全局统计的影响。
所以，第一步就是要建立这样一个『窗口』。当然了，首先要安装包，我们要用到的R包有：
library(spgwr)
library(rgdal)
library(sf)
library(spData)
library(sp)
library(lattice)
library(ggplot2)
library(ggthemes)
其中，spgwr是做GWR的包，rgdal是用来读取矢量要素的，sf，sp和spData都是用来处理矢量数据的，别的基本都是画图用。
以下默认你会R和GWR的基本操作。并且，以下只展现方法，不要纠结我的数据和结果，我随便找的数据，这个数据本身没有什么意义，所以做出的统计看起来很『壮观』。
我们先导入数据。这里我用的是美国本土48州各个县（county，也有翻译成郡的）的人口普查数据和农业数据，来源是ESRI Online数据库。为啥用这个数据呢？因为...我电脑里面就存了这么个可以用来做GWR的数据...
我们用rgdal读取数据，然后把它画出来看看
require(rgdal)
usa_agri <- readOGR(dsn = "~/Documents/Spatial", layer = "usa_counties")
plot(usa_agri)
会得到这个东西：
readOGR里面，dsn后面加储存shp的路径（加到文件夹为止），layer后面写shp的文件名（不加.shp）。不喜欢rgdal的同学可以不用，用maptools或者spData等别的处理shp的R包代替。不过如果用maptools，要注意处理一下参考系。
我们看一下这个shp里面的列联表都有什么：
可见，shp里面有3108个县的数据，数据有61种。然后再看data下面有什么：
总之就是各种人口普查的数据，后面截不完图，还有经济、房地产和农业之类的数据。那我们就随便选两个来当变量。我就随便挑了，因变量选AVESIZE12，即2012年各个县农场的平均占地面积。自变量选POP_SQMI，也就是人口密度（每平方英里的人口）。
现在正式建立窗口，调用的是spgwr里面的gwr.sel函数：
bw <- gwr.sel( AVE_SIZE12 ~ POP_SQMI, data = usa_agri, gweight = gwr.Gauss,
verbose = FALSE, method = "cv")
其中~前后分别是因变量和自变量。GWR里因变量只能有1个，但自变量可以选多个，如果需要多个自变量的话，就在代码POP_SQMI之后用+号连接就行。gweight是你的空间加权的函数（随空间距离增大而不断衰减的函数，衰减率由下面要提到的带宽控制），这里用的是比较常用的高斯函数，其余的还有gwr.bisquare等函数可以调用。verbose决定是否汇报制定窗口的过程。method是决定构建带宽窗口模型的方法，这里用的cv指的是cross validation，即交叉验证法，也是最常用的方法，简单说就是把数据分成不同的组，分别用不同的方法来做回归计算，计算完了之后记录下结果，然后打乱重新分组，再回归计算，再看结果，周而复始，最后看哪种计算方法的结果最靠谱，这种方法就是最优解。还有一种很常见的选择最佳拟合模型的方法是AIC optimisation法，把method后面的cv改成aic就可以用。具体AIC optimisation是什么：AIC（赤池信息准则）_百度百科。总之，空间加权函数和带宽窗口构建方法的选择是GWR里面十分重要的步骤。
以上便是固定带宽窗口的示意图。比如我在对佐治亚做GWR，这一轮的regression target是红色的这个县，根据做出来的窗口，圆圈以内的县都要被算为红色县的邻县，其权重根据高斯函数等空间权重函数来赋值，而圆圈以外的县，空间权重都赋为0。
不喜欢固定带宽窗口的同学也可以不用它，而是用符合Tobler地理学第一定律的非固定带宽邻域统计，操作方法是在gwr.sel里面加一个命令adapt = TRUE，这样的情况下，根据你设置的k邻居数，每一轮统计的时候，和本轮对象在k以内相邻的多边形的权重参数会被赋值为0到1之间的一个数，比如下图:
我在对佐治亚做GWR，这一轮的regression target是红色的这个县，那么图上标为1的县就是红色县的1阶邻县，标为2的是2阶（邻县的邻县），标为3的是3阶（邻县的邻县的邻县）。如果用非固定带宽邻域统计，k为3，那么1、2、3都被定义为红色县的邻县，它们的权重从3到1依次增加，会按比例被赋上0和1之间的值，而其它没有标注的县，权重为0。
下一步就是用前一步做出的窗口去扫过全局区域：
gwr_result <- gwr(AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
这一步如果数据量大，可能会要跑一阵，跑完之后我们看看结果里面有什么：
Call:
gwr(formula = AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 205880.3
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu.Max. Global
X.Intercept. 7.3883e+01 2.1081e+02 3.2802e+02 6.6691e+02 8.5705e+03 625.5656
POP_SQMI -8.0085e+01 -4.5983e-01 -1.4704e-01 -7.3703e-02 -2.1859e-03 -0.0426
Number of data points: 3108
Effective number of parameters (residual: 2traceS - traceS'S): 119.6193
Effective degrees of freedom (residual: 2traceS - traceS'S): 2988.381
Sigma (residual: 2traceS - traceS'S): 1048.78
Effective number of parameters (model: traceS): 84.90185
Effective degrees of freedom (model: traceS): 3023.098
Sigma (model: traceS): 1042.741
Sigma (ML): 1028.4
AICc (GWR p. 61, eq 2.33p. 96, eq. 4.21): 52109.55
AIC (GWR p. 96, eq. 4.22): 52017.7
Residual sum of squares: 3287040139
Quasi-global R2: 0.4829366
基本上你做GWR该需要的结果这里都有了。比如窗口大小（Fixed bandwidth）是205880.3，意思是前一步构建的带宽窗口是半径205.88千米的圆。Effective number of parameters显示的是你带宽窗口的大小合不合适。Sigma是残差的标准差，这个值要尽量小。Residual sum of squares（RSS）也是对拟合程度的一个评估值。最重要的是最后那个R2，越靠近1说明统计的拟合度越好。我这里面Sigma很大，R2也不是很大，因为我这里只是呈现方法，用的数据本来就是互不相干、没什么太大意义的，所以不用太纠结。如果你是真正的统计数据要来做GWR，就需要注意这些值了。
然后，我们就可以把每个县的R2画在地图上。首先，前面报告里的这些数据，比如R2，要先自己去生成的GWR结果里面去找，然后自己再算一下每个县的local R2，并把它们赋值到shp里面去：
R语言可以使用spatstat包中的spatstat.geom()函数来实现地理加权回归投影。spatstat.geom()函数接受一个经纬度矢量作为参数，然后根据给定的经纬度坐标，计算出一组投影坐标。这些投影坐标可以用于地理加权回归分析。难，空间回归模型中的回归系数不随空间位置而变化，因此空间回归模型是全局模型。但是由于空间异质性和空间非平稳性，不同空间子区域中自变量和因变量的关系很可能不同，因此就需要处理空间异质性的局部空间回归方法，因此就有了地理加权回归模型（Geographically Weighted Regression, GWR）的提出。地理加权回归同时考虑了空间的相关性与异质性。地理加权回归模型一般形式如下：
第 1 页
防火板首选福建原时态
A级防火板生产厂家，技术领先品质保障，尺寸齐全，厂家直销 。可按需定制。详情欢迎来电咨询
点击立即咨询，了解更多详情
咨询
福建原时态建筑工程 广告
；
其中，与是因变量与自变量在处的观测值，为第个采样点的坐标，是第个采样点上的第个回归参数，点被称为回归点。，当时，。简便记为：
第 2 页
地理加权回归会得出个回归方程，对应每个回归点，都会有一个回归方程。若，则地理加权回归模型退化为普通线性回归模型。
模型回归参数需要通过局部加权最小二乘进行估计。假设在位置的权重为，（），那么位置的参数估计为使达到最小的值。
上式中的为权函数，反映其他观测点的样本对回归点的影响。权函数值越大，影响越大。该值通常由对应两点之间的距离决定。常用的权函数有：
第 3 页
（1）Gauss函数：
（2）bi-square函数：
以上两式中的为点到点的距离。被称为带宽，是需要人工选择的参数。这两个函数中，距离越大，函数值都越小。这说明选择这两个函数时，都假设观测点越远，影响越小。
第 4 页
地理加权回归模型将数据的空间位置嵌入到回归参数中，因此考虑了空间的异质性。同时，由于不同位置的观测点对回归参数的影响大小不同（通常离回归点越近，影响越大），因此该模型也考虑了空间相关性。地理加权回归的使用，也当同时以空间相关性与异质性为前提。如果没有空间相关性，那么该模型就缺乏合理性。如果没有空间异质性，那么该模型就缺乏必要性。
【案例72】
以2011年北京、天津、唐山各县（区）疾病发病率和同期各县（区）的人口密度、人均GDP、年平均风速、光照强度
第 5 页
、相对湿度、年降水量等数据为例，研究社会经济和气象因素对该疾病发病率在各地区影响的不稳定性。
本节所用的数据为2011年某地区某疾病的发病率（C:\Example\Data\7.3WGR\地区汇总.csv），命名为“地区汇总.csv”。此处只展示前10条
表 73 2011年某地区的某疾病发病情况
code
precp
relHum
sunShn
wndspd
PopuDen
PerGdp
rate
110101
58.06842
49.82301
206.0755
2.195716
15318.95
0.002438
347.1859
110102
62.59835
52.05048
202.8901
2.169381
14956.04
0.003675
447.7185
110105
57.65938
50.06015
205.5302
2.189571
8528.913
0.00096
1139.069
110106
55.73981
50.20182
205.5477
2.215026
7698.981
0.000419
1759.943
110107
52.39305
50.77886
205.6759
2.255836
8249.254
0.000541
1561.691
110108
52.13075
51.01494
207.5023
2.255693
8524.685
0.001021
1010.973
110109
44.66608
53.70504
218.3655
2.473628
224.6719
0.00038
758.6155
110111
46.62295
52.79277
208.7774
2.333494
532.4592
0.000457
2082.544
110112
55.23023
51.21378
201.2761
2.186171
1454.807
0.00036
2406.252
110113
50.22722
55.34219
195.5819
2.014273
967.6505
0.001211
1643.1
第 6 页
表 74各变量的含义
变量名称
变量含义
单位
code
地理编码
无
precp
降水量
毫米
relHum
相对湿度
%
sunShn
日照强度
瓦/平方米
wndspd
风速
千米/时
popuDen
人口密度
千人/平方公里
PerGdp
人均国内生产总值
千元
rate
发病率
此外，还需要的数据是包含该地区所有区县的地图文件（.dbf文件和.shp文件
从光盘中获取，C:\Example\Data\Geodata\JJT）。该文件也可以通过arcgis软件从全国各区县地图中选择生成。
（2）采用R语言建立地理加权回归模型
第一步，加载如下程序包，代码如下：
library(spgwr)
第二步，导入所需的数据，代码如下：
hData <- read.csv("C: /Example/Data/7.3WGR/地区汇总.csv ")#导入发病率和影响因素的数据
dbf <- read.dbf("C: /Example/Data/Ge
odata /JJT.dbf")#导入地图的数据（dbf格式）
第三步，将导入的两组数据合并，代码如下：
Data <- merge(hData,dbf, by.x="code" , by.y = "CNTY_CODE" , all.x =T)
第四步，确定带宽，采用gwr.sel函数。所使用的参数如下：
formula：模型公式，用于指出因变量与自变量；
data：自变量与因变量取值的数据集；
coords:代表空间观测值位置的坐标矩阵。
代码如下：
col.bw <- gwr.sel(rate ~ PopuDen + PerGdp+precp+relHum+sunShn+wndspd, data=data, coords=cbind(data$x, data$y)) #利用交叉验证选择最优带宽
第五步，生成地理加权回归模型，采用gwr函数，使用的各参数意义如下：
formula：模型公式，用于指出因变量与自变量；
data：自变量与因变量取值的数据集；
coords:代表空间观测值位置的坐标矩阵；
第 10 页
bandwidth：带宽，由上步gwr.sel生成；
gweight:不指定时，默认使用高斯函数确定权重矩阵；
hatmatrix:如果为TRUE，帽子矩阵作为结果的一部分返回。
代码如下：
col.gauss <- gwr(rate ~ PopuDen + PerGdp+precp+relHum+sunShn+wndspd, data=data, coords=cbind(data$x, data$y), bandwidth=col.bw, hatmatrix=TRUE)
col.gauss
此外，地理加权回归还经常使用bi-square权函数，该方法和Gauss权函数方法相似，这里只将实现代码列出，代码如下：
col.d <- gwr.sel(rate ~ PopuDen + PerGdp+precp+relHum+sunShn+wndspd, data=data, coords=cbind(data$x, data$y), gweight=gwr.bisquare) #确定带宽
col.bisq <- gwr(rate ~ PopuDen + PerGdp+precp+relHum+sunShn+wndspd, data=data, coords=cbind(data$x, data$y), bandwidth=col.d, gweight=gwr.bisquare, hatmatrix=T) #建立地理加权回归模型
第 12 页
col.bisq #结果展示
（3）结果分析
建立的地理加权回归截距和系数统计如表 75所示：
表 75 加权回归系数统计表
变量
最小
四分之一分位数
中位数
四分之三分位数
最大值
全局
截距
-3478.00
-2391.00
896.50
4100.00
11420.00
1223.24
PopuDen
-0.05
-0.04
-0.04
-0.03
0.01
-0.04
PerGdp
-191700.00
-106600.00
-24080.00
-1769.00
44410.00
-37167.70
precp
-0.68
-0.34
-0.26
0.28
2.29
-0.08
relHum
-133.90
-72.10
-21.84
32.26
52.46
-31.49
sunShn
-16.88
0.50
0.84
1.45
2.94
0.66
wndspd
-91.81
357.90
625.00
667.20
874.70
775.70
从计算结果可以看出，由于地理加权回归得出了影响因素在每个地区的影响系数，各个因素对每个地区的影响程度并不相同，若系数的变化范围较大，说明该影响因素总体上影响程度有很大的不稳定性，若系数的变化范围较小，说明该影响因素总体上影响比较稳定。从该案例可以看出，当月人均国内生产总值对疾病的发病率影响最大，并且在大部分地区呈负相关关系，即当月人均国内生产总值越高，疾病发病率越低，此外该变量的系数变化范围十分大，说明该因素在不同地区的影响程度有很大差别，具有不稳定性。其次，风速对该病的发病率影响也较大，但风速与该病的
第 14 页
发病率大部分呈正相关，即风速越大，该病的发病率越大，可以推断该病可以通过空间传染。在6个影响因素中人口密度对疾病的发病率影响最低。
地理加权回归模型的R2为0.47，即该模型能解释疾病发病率总变异的47%，比全局普通线性回归的R2（0.23）大一倍。因此对于存在空间相关性的变量，应该使用地理加权回归进行计算。
空间回归与地理加权回归的比较：
空间回归与地理加权回归都是在经典的回归模型上考虑事物的空间属性，从而发展出来的。
空间回归通常只考虑空间相关性，而地
理加权回归同时考虑了空间相关性与异质性。
从本质上说，两种模型考虑空间相关性的角度也不同。空间回归模型实质是考虑值的空间相关性。而地理加权回归模型考虑的是数量关系或规律的空间相关性。
与只考虑空间相关性的空间回归模型相比，同时考虑空间相关性与异质性的地理加权回归模型显然考虑得更为细致。但是地理加权回归模型的求解更为复杂，且结果更难解读。比如需要逐个解读n个回归方程。这意味着通过地理加权回归较难掌握全局的规律。因此，空间回归模型更善于刻画全局规律。
