其他分享
首页 > 其他分享> > 地理加权归回模型 (GWR) 参数估计

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

作者：互联网

 

作者：陈凤 (西安交通大学)
 
Stata连享会   计量专题 || 精品课程 || 简书推文 || 公众号合集

点击查看完整推文列表

连享会：内生性问题及估计方法专题

连享会-内生性专题现场班-2019.11.14-17

 

 

连享会 空间计量专题 2019.10.24-27

2019.10-杨海生-空间计量专题-连享会

 

文章目录

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

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

(1)Yi=β0(ui,vi)+∑j=1pβj(ui,vi)Xij+εiY_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 为例， 如果不考虑房屋的地理位置信息，可以建立一个简单的线性回归模型：

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

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

(3)Yi=β(ui,vi)Xi+εiY_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)，即

(4)min⁡∑i=1n[yi−∑j=1pβj(u0,v0)xij]2wi(u0,v0)\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​) 处的局部最小二乘估计值为

(5)β^(u0,v0)=(XTW(u0,v0)X)−1XTW(u0,v0)Y\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 估计值和如下线性模型的最小二乘估计是等价的：

(6)wi(u0,v0)Yi=∑j=1pwi(u0,v0)xijβj(u0,v0)+εˉi,i=1,2,⋯ ,n.\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∑p​wi​(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π​1​exp(−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∣&lt;h;0,∣dij∣&gt;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, &amp; \left|d_{ij}\right|&lt;h; \\ 0,&amp; \left|d_{ij}\right|&gt;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​)如下所示。
image.png

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∣&lt;dik;0,∣dij∣&gt;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, &amp; \left|d_{ij}\right|&lt;d_{ik}; \\ 0,&amp; \left|d_{ij}\right|&gt;d_{ik}. \end{array}\right. wi​(uj​,vj​)=Kh​(dij​)=⎩⎨⎧​[1−(dik​dij​​)2]2,0,​∣dij​∣<dik​;∣dij​∣>dik​.​

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

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)=h1​i=1∑n​(Yi​−Y^(−i)​(h))2,

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

h0=arg⁡min⁡h&gt;0CV(h)h_0=\arg\min\limits_{h&gt;0} {\rm CV}(h)h0​=argh>0min​CV(h)

 

连享会计量方法专题……

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

设

Y^(h)=(Y^1(h),Y^2(h),⋯&ThinSpace;,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)))2n​i=1∑n​(Yi​−Y^i​(h))2,

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

h0=arg⁡min⁡h&gt;0GCV(h)h_0=\arg\min\limits_{h&gt;0} {\rm GCV}(h)h0​=argh>0min​GCV(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&gt;0AICc(h)h_0=\arg\min\limits_{h&gt;0}{\rm AICc}(h)h0​=argh>0min​AICc(h)

AICc\rm AICcAICc 准则选择最优窗宽如下所示：
image.png
注：模拟实验以及经验表明，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 的回归系数的局部估计值如下：
Dublin.jpg

5. 参考文献

 

关于我们

点击此处-查看完整推文列表

在这里插入图片描述

标签：归回,boldsymbol,u0,v0,beta,参数估计,GWR,rm,left
来源： https://blog.csdn.net/arlionn/article/details/100641980