地理加权归回模型 (GWR) 参数估计

作者：陈凤 (西安交通大学)

Stata连享会计量专题 || 精品课程 || 简书推文 || 公众号合集

连享会计量方法专题……

文章目录

- 连享会计量方法专题……
1. 地理加权回归模型简介
2. 地理加权回归模型的参数估计方法
- 连享会计量方法专题……
3. 常用的核函数
- 3.1. Gussian kernel function
- 3.2. Bi-square kernel function
- 3.3. K-nearest neighbor kernel function
4. 窗宽h的选择准则
- 4.1. 交叉确认方法(Cross-validation (CV) criterion)
- - 连享会计量方法专题……
- 4.2. 广义交叉确认方法(Generalized cross-validation (GCV) criterion)
- 4.3. AICc信息准则(Corrected Akaike information criterion (AICc))
- - 连享会计量方法专题……
4. 在 R 软件运行地理加权回归模型
5. 参考文献
- 关于我们

1. 地理加权回归模型简介

空间数据在地理学、经济学、环境学、生态学以及气象学等众多领域中广泛存在。根据 Tobler 提出的「地理学第一定律」：任何事物之间都是空间相关的，距离越近的事物之间的空间相关性越大。因此，不同于传统的截面数据，空间数据的空间相关性会导致回归关系的空间非平稳性 (空间异质性)。为了探索空间数据的空间非平稳性， Brunsdon 等 (1996) 首次提出了 地理加权回归模型，设定如下：

Yi=β0(ui,vi)+∑j=1pβj(ui,vi)Xij+εi(1)Y_i=\beta_0(u_i,v_i)+\sum_{j=1}^p\beta_j(u_i,v_i)X_{ij}+\varepsilon_i \tag{1}Yi=β0(ui,vi)+j=1∑pβj(ui,vi)Xij+εi(1)

其中，βj(u,v)(j=0,1,⋯,p)\beta_j(u,v)\ (j=0,1,\cdots,p)βj(u,v) (j=0,1,⋯,p) 为 「空间地理位置函数」。

以某城市的房屋价格 YYY 和房屋面积 XXX 为例，如果不考虑房屋的地理位置信息，可以建立一个简单的线性回归模型：

Yi=βXi+εi(2)Y_i=\beta X_i+\varepsilon_i \tag{2}Yi=βXi+εi(2)

其中，β\betaβ 为房屋的单位面积均价。实际中，处于不同位置的房屋价格可能会相差甚远，但是模型 (2)(2)(2) 却不能反映出这种异质性。因此，为了能够描述不同位置房屋价格的差异性，我们可以建立如下模型：

Yi=β(ui,vi)Xi+εi(3)Y_i=\beta(u_i,v_i)X_i+\varepsilon_i \tag{3}Yi=β(ui,vi)Xi+εi(3)

其中，β(u,v)\beta(u,v)β(u,v) 是地理位置的函数。相比于模型 (2)(2)(2)，模型 (3)(3)(3) 可以反映房屋价格随地理位置的变化而变化的规律。

上述例子说明有必要对空间数据建立地理加权回归模型来探索空间数据的非平稳性。

2. 地理加权回归模型的参数估计方法

根据 Tobler 地理学第一定律，距离越近的事物之间的相关性越大。故对于一个给定的地理位置(u0,v0)(u_0,v_0)(u0,v0)，可以采用局部加权最小二乘来估计 βj(u0,v0)(j=0,1,⋯,p)\beta_j(u_0,v_0)\ (j=0,1,\cdots,p)βj(u0,v0) (j=0,1,⋯,p)，即

min⁡∑i=1n[yi−∑j=1pβj(u0,v0)xij]2wi(u0,v0)(4)\min \sum_{i=1}^n [y_i-\sum_{j=1}^p\beta_j(u_0,v_0)x_{ij}]^2w_i(u_0,v_0) \tag{4} mini=1∑n[yi−j=1∑pβj(u0,v0)xij]2wi(u0,v0)(4)

其中，{wi(u0,v0)}i=1n\{w_i(u_0,v_0)\}_{i=1}^n{wi(u0,v0)}i=1n 是在地理位置 (u0,v0)(u_0,v_0)(u0,v0) 处的空间权重。令 β(u0,v0)=(β0(u0,v0),β1(u0,v0),⋯,βp(u0,v0))T,\boldsymbol \beta(u_0,v_0)=(\beta_0(u_0,v_0),\beta_1(u_0,v_0),\cdots,\beta_p(u_0,v_0))^{\rm T},β(u0,v0)=(β0(u0,v0),β1(u0,v0),⋯,βp(u0,v0))T, 则 β(u0,v0)\boldsymbol \beta(u_0,v_0)β(u0,v0) 在 (u0,v0)(u_0,v_0)(u0,v0) 处的局部最小二乘估计值为

β^(u0,v0)=(XTW(u0,v0)X)−1XTW(u0,v0)Y(5)\hat{\boldsymbol \beta}(u_0,v_0)=\left(\boldsymbol X^{\rm T}\boldsymbol W(u_0,v_0)\boldsymbol X\right)^{-1}\boldsymbol X^{\rm T}\boldsymbol W(u_0,v_0)\boldsymbol Y \tag{5}β^(u0,v0)=(XTW(u0,v0)X)−1XTW(u0,v0)Y(5)

其中，

X=(X0,X1,⋯,Xp),Xj=(x1j,x2j,⋯,xnj)T;\boldsymbol X=(\boldsymbol X_0,\boldsymbol X_1,\cdots,\boldsymbol X_p), \boldsymbol X_j=(x_{1j},x_{2j},\cdots,x_{nj})^{\rm T};X=(X0,X1,⋯,Xp),Xj=(x1j,x2j,⋯,xnj)T;

Y=(Y1,Y2,⋯,Yn)T;\boldsymbol Y=(Y_1,Y_2,\cdots,Y_n)^{\rm T}; Y=(Y1,Y2,⋯,Yn)T;

W(u0,v0)=Diag(w1(u0,v0),w2(u0,v0),⋯,wn(u0,v0)).\boldsymbol W(u_0,v_0)={\rm Diag}\left(w_1(u_0,v_0),w_2(u_0,v_0),\cdots,w_n(u_0,v_0)\right).W(u0,v0)=Diag(w1(u0,v0),w2(u0,v0),⋯,wn(u0,v0)).

令 (u0,v0)=(ui,vi),i=1,2,⋯,n(u_0,v_0)=(u_i,v_i), i=1,2,\cdots,n(u0,v0)=(ui,vi),i=1,2,⋯,n，则可以由公式 (5)(5)(5) 得到回归函数 β(u,v)\boldsymbol \beta(u,v)β(u,v)在所有观测位置处的局部估计值。

注1： βj(u,v)(j=0,1,⋯,p)\beta_j(u,v)\ (j=0,1,\cdots,p)βj(u,v) (j=0,1,⋯,p) 可以在任意位置处被估计。因此，GWR 模型也可以作为空间数据的插值工具。
注2： 在 (u0,v0)(u_0,v_0)(u0,v0) 处，β(u0,v0)\boldsymbol \beta(u_0,v_0)β(u0,v0) 的 GWR 估计值和如下线性模型的最小二乘估计是等价的：

wi(u0,v0)Yi=∑j=1pwi(u0,v0)xijβj(u0,v0)+εˉi,i=1,2,⋯,n.(6)\sqrt{w_i(u_0,v_0)}Y_i=\sum_{j=1}^p\sqrt{w_i(u_0,v_0)}x_{ij}\beta_j(u_0,v_0)+\bar{\varepsilon}_i, i=1,2,\cdots,n. \tag{6}wi(u0,v0)Yi=j=1∑pwi(u0,v0)xijβj(u0,v0)+εˉi,i=1,2,⋯,n.(6)

连享会计量方法专题……

3. 常用的核函数

在核光滑方法中，常用的核函数如下：

3.1. Gussian kernel function

wi(uj,vj)=Kh(dij)=12πexp⁡(−12(dijh)2),w_i(u_j,v_j)=K_h(d_{ij})=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{d_{ij}}{h}\right)^2\right),wi(uj,vj)=Kh(dij)=2π1exp(−21(hdij)2),

其中，hhh 为窗宽，dijd_{ij}dij 为点 (ui,vi)(u_i,v_i)(ui,vi) 和 (uj,vj)(u_j,v_j)(uj,vj) 之间的距离。

3.2. Bi-square kernel function

wi(uj,vj)=Kh(dij)={[1−(dijh)2]2,∣dij∣<h;0,∣dij∣>h.w_i(u_j,v_j)=K_h(d_{ij})=\left\{ \begin{array}{rcl} \left[1-\left(\frac{d_{ij}}{h}\right)^2\right]^2, & \left|d_{ij}\right|<h; \\ 0,& \left|d_{ij}\right|>h. \end{array}\right. wi(uj,vj)=Kh(dij)=⎩⎨⎧[1−(hdij)2]2,0,∣dij∣<h;∣dij∣>h.

给一个hhh，(ui,vi)(u_i,v_i)(ui,vi)处自变量的观测值对(u0,v0)(u_0,v_0)(u0,v0)处因变量的权重wi(u0,v0)w_i(u_0,v_0)wi(u0,v0)如下所示。

3.3. K-nearest neighbor kernel function

给定KKK，dikd_{ik}dik为(ui,vi)(u_i,v_i)(ui,vi)到第KKK个邻近点的距离，则

wi(uj,vj)=Kh(dij)={[1−(dijdik)2]2,∣dij∣<dik;0,∣dij∣>dik.w_i(u_j,v_j)=K_h(d_{ij})=\left\{ \begin{array}{rcl} \left[1-\left(\frac{d_{ij}}{d_{ik}}\right)^2\right]^2, & \left|d_{ij}\right|<d_{ik}; \\ 0,& \left|d_{ij}\right|>d_{ik}. \end{array}\right. wi(uj,vj)=Kh(dij)=⎩⎨⎧[1−(dikdij)2]2,0,∣dij∣<dik;∣dij∣>dik.

对于任意的观测点来说，KKK近邻核函数总是保持有KKK个观测点的空间权重不为零，如下所示。

4. 窗宽h的选择准则

在地理加权回归模型中，常用的最优窗宽选取准则有交叉确认方法、广义交叉确认方法以及AICc信息准则。这三种准则的定义分别如下所示。

4.1. 交叉确认方法(Cross-validation (CV) criterion)

交叉确认方法的具体过程如下：给定一个 hhh, 去掉第 iii 组观测值 (Yi,Xi)(Y_i,X_i)(Yi,Xi)，用剩下的 (n−1)(n-1)(n−1) 组数据在给定的 hhh下进行地理加权回归参数估计，然后得到在 XiX_iXi 处的拟合值Y^(−i)(h)\hat{Y}_{(-i)}(h)Y^(−i)(h)。令

CV(h)=1h∑i=1n(Yi−Y^(−i)(h))2,{\rm CV}(h)=\frac{1}{h}\sum_{i=1}^n\left(Y_i-\hat{Y}_{(-i)}(h)\right)^2,CV(h)=h1i=1∑n(Yi−Y^(−i)(h))2,

则最优窗宽 h0h_0h0 的选取如下：

h0=arg⁡min⁡h>0CV(h)h_0=\arg\min\limits_{h>0} {\rm CV}(h)h0=argh>0minCV(h)

连享会计量方法专题……

4.2. 广义交叉确认方法(Generalized cross-validation (GCV) criterion)

设

Y^(h)=(Y^1(h),Y^2(h),⋯,Y^n(h),)T=L(h)Y,\hat{\boldsymbol Y}(h)=\left(\hat{Y}_1(h),\hat{Y}_2(h),\cdots,\hat{Y}_n(h),\right)^{\rm T}=\boldsymbol L(h)\boldsymbol Y,Y^(h)=(Y^1(h),Y^2(h),⋯,Y^n(h),)T=L(h)Y,

其中，L(h)\boldsymbol L(h)L(h) 为 “帽子” 矩阵，令

GCV(h)=n(n−tr(L(h)))2∑i=1n(Yi−Y^i(h))2,{\rm GCV}(h)=\frac{n}{\left(n-tr(\boldsymbol L(h))\right)^2}\sum_{i=1}^n\left(Y_i-\hat{Y}_{i}(h)\right)^2,GCV(h)=(n−tr(L(h)))2ni=1∑n(Yi−Y^i(h))2,

则最优窗宽 h0h_0h0 的选择标准为：

h0=arg⁡min⁡h>0GCV(h)h_0=\arg\min\limits_{h>0} {\rm GCV}(h)h0=argh>0minGCV(h)

4.3. AICc信息准则(Corrected Akaike information criterion (AICc))

令 Y^(h)=L(h)Y\hat{\boldsymbol Y}(h)=\boldsymbol L(h)\boldsymbol YY^(h)=L(h)Y， ε^=YT(In−L(h))T(In−L(h))Y\hat{\boldsymbol\varepsilon}=\boldsymbol Y^{\rm T}(\boldsymbol I_n-\boldsymbol L(h))^{\rm T}(\boldsymbol I_n-\boldsymbol L(h))\boldsymbol Yε^=YT(In−L(h))T(In−L(h))Y，则有

AICc(h)=log⁡(1nε^Tε^)+n+tr(L(h))n−2−tr(L(h)).{\rm AICc}(h)=\log\left(\frac{1}{n}\hat{\boldsymbol\varepsilon}^{\rm T}\hat{\boldsymbol\varepsilon}\right)+\frac{n+tr(\boldsymbol L(h))}{n-2-tr(\boldsymbol L(h))}.AICc(h)=log(n1ε^Tε^)+n−2−tr(L(h))n+tr(L(h)).

最优窗宽 h0h_0h0 的选取如下：

h0=arg⁡min⁡h>0AICc(h)h_0=\arg\min\limits_{h>0}{\rm AICc}(h)h0=argh>0minAICc(h)

AICc\rm AICcAICc 准则选择最优窗宽如下所示：

注：模拟实验以及经验表明，CV\rm CVCV 和 GCV\rm GCVGCV 准则一般会趋于确定一个稍微偏小的窗宽h0h_0h0，而较小的窗宽会使得回归函数的估计值偏差减小，但是方差会增大。因此，会出现过拟合现象。但是对于 AICc\rm AICcAICc 在很多情况下可以较好的克服过拟合现象，即趋于确定一个更合理的窗宽。

连享会计量方法专题……

4. 在 R 软件运行地理加权回归模型

在R软件中，可以调用 Rpacakge-GWmodel (Lu, B. B, et al. 2014) 来实现地理加权回归模型的参数过程。以都柏林 2014 年的选举数据为例，下面介绍 GWR 在 R 中的实现过程：

# 1. 加载 Rpacakge：
library("GWmodel")# 2. 加载数据
data(DubVoter)# 3. 选择最优窗宽
Dub<cbind(Dub.voter\$DiffAdd,Dub.voter\$LARent,Dub.voter\$SC1,Dub.voter\$Unempl,Dub.voter\$LowEduc,Dub.voter\$Age18_24,Dub.voter\$Age25_44,Dub.voter\$Age45_64,Dub.voter\$GenEl2004)
DubCoord<-cbind(Dub.voter\$X,Dub.voter\$Y) %地理位置
DIS<-gw.dist(dp.locat=DubCoord)% 计算距离矩阵
bw1<bw.gwr(GenEl2004~DiffAdd+LARent+SC1+Unempl+LowEduc+Age18_24+Age25_44+Age45_64, approach="AICc",adaptive=TRUE, data=Dub.voter, kernel = "bisquare",dMat=DIS) %选择bi-square函数作为核函数，使用AICc准则选择最优窗宽# 4. 拟合GWR模型
gwr.res1<gwr.basic(GenEl2004~DiffAdd+LARent+SC1+Unempl+LowEduc+Age18_24+Age25_44+Age45_64, data=Dub.voter, bw=bw1,adaptive=TRUE,kernel = "bisquare", dMat=DIS)# 5. 将局部估计值画在对应的地图上面
library("RColorBrewer")
Mcolor<-1;mypalette <- colorRampPalette(brewer.pal(9, "Greys"))(200)
map.na=list("SpatialPolygonsRescale", layout.north.arrow(),offset=c(329000, 261500), scale=4000,col=1)
map.scale.1=list("SpatialPolygonsRescale",layout.scale.bar(),offset=c(326500, 217000), scale=5000,col=1,fill=c("transparent","black"))
map.scale.2=list("sp.text",c(326500, 217900),"0",cex=0.9,col=1)
map.scale.3=list("sp.text",c(331500, 217900),"5km",cex=0.9,col=1)
map.layout<-list(map.na,map.scale.1,map.scale.2,map.scale.3)
mypalette.9<-brewer.pal(9,"Greys")
spplot(gwr.res1$SDF, "LowEduc",key.space="right",col.regions=mypalette.6, at=c(-8,-6,-4,-2,0,2,4),sp.layout=map.layout)

在都柏林 2014 年选择数据中，使用 AICc 准则确定的最优窗宽为 115，其中变量 LowEduc 的回归系数的局部估计值如下：

5. 参考文献

Brunsdon, C. E, Fotheringham, A. S. and Charlton, M. E., 1999. Some notes on parametric significance test for geographically weighted regression. Journal of Regional Science, 39 (3): 497–524. [PDF]
梅长林, 王宁. 近代回归分析方法 [M]: 北京：科学出版社, 2012.
Lu, B. B., Harris, P., Charlton, M. and Brunsdon, C., 2014. The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models. Geo-spatial Information Science, 17 (2): 85–101. [PDF]

关于我们

「Stata 连享会」 由中山大学连玉君老师团队创办，定期分享实证分析经验，公众号：StataChina。
公众号推文同步发布于 CSDN 、简书和知乎Stata专栏。可在百度中搜索关键词「Stata连享会」查看往期推文。
点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
欢迎赐稿： 欢迎赐稿。录用稿件达三篇以上，即可免费获得一期 Stata 现场培训资格。
E-mail： StataChina@163.com
往期推文：计量专题 || 精品课程 || 简书推文 || 公众号合集