spgwr | R语言与地理加权回归（Ⅰ-1）：线性地理加权回归[通俗易懂]

地理加权回归（Geographically Weighted Regression, GWR）经过多年发展，已经具备了多种形式，在R语言中也对应着多个工具包，其中spgwr是一个开发较早、比较经典的工具包，功能也相对基础。

library(spgwr)

在该包中，运行线性地理加权回归的函数是gwr()。语法结构如下：

gwr(formula, data = list(), coords,
    bandwidth, gweight = gwr.Gauss, 
    adapt = NULL, hatmatrix = FALSE, fit.points,
    longlat = NULL, se.fit = FALSE, weights,
    cl = NULL, predictions = FALSE, 
    fittedGWRobject = NULL, se.fit.CCT = TRUE)

本篇先介绍它的几个主要的参数：

该包目前的版本号是0.6-34，还不支持sf格式的矢量对象。

数据源

library(rgdal)
NY8 <- readOGR(system.file("shapes/NY8_utm18.shp",
                           package = "spData"))

模型形式

form <- Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME

带宽

带宽（bandwidth）确定了局部的范围，该包的gwr.sel()函数提供了两种确定带宽的方法：交叉验证法和AIC信息准则法。语法结构如下：

gwr.sel(formula, data = list(),
        coords, adapt = FALSE, gweight = gwr.Gauss,
        method = "cv", verbose = TRUE,
        longlat = NULL, RMSE = FALSE, weights,
        tol = .Machine$double.eps^0.25,
        show.error.messages = FALSE)

• method：可选项有cv（交叉验证）、aic（AIC准则）。

bw <- gwr.sel(formula = form, data = NY8,
              gweight = gwr.Gauss, method = "cv")
bw
## [1] 179942.6

formula、data、gweight参数需要与gwr()函数的对应参数保持一致。

距离加权函数

距离加权函数是一个随距离增加而逐渐衰减的函数，该包提供了4种地理加权函数：gwr.gauss、gwr.Gauss（默认）、gwr.bisquare、gwr.tricube。

以b表示带宽，d表示距离。以d = 100为例：

完整形式

线性回归：

model.lm <- lm(formula = form, data = NY8@data)
summary(model.lm)
## 
## Call:
## lm(formula = form, data = NY8@data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7417 -0.3957 -0.0326  0.3353  4.1398 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.51728    0.15856  -3.262  0.00124 ** 
## PEXPOSURE    0.04884    0.03506   1.393  0.16480    
## PCTAGE65P    3.95089    0.60550   6.525 3.22e-10 ***
## PCTOWNHOME  -0.56004    0.17031  -3.288  0.00114 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6571 on 277 degrees of freedom
## Multiple R-squared:  0.1932, Adjusted R-squared:  0.1844 
## F-statistic:  22.1 on 3 and 277 DF,  p-value: 7.306e-13

线性地理加权回归：

library(rgdal)
NY8 <- readOGR(system.file("shapes/NY8_utm18.shp",
                           package = "spData"))
form <- Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME
bw <- gwr.sel(formula = form, data = NY8,
              gweight = gwr.Gauss, method = "cv")

model <- gwr(formula = form, data = NY8,
             bandwidth = bw, gweight = gwr.Gauss)
model
## Call:
## gwr(formula = form, data = NY8, bandwidth = bw, gweight = gwr.Gauss)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 179942.6 
## Summary of GWR coefficient estimates at data points:
##                   Min.   1st Qu.    Median   3rd Qu.      Max.  Global
## X.Intercept. -0.522172 -0.520740 -0.520154 -0.514439 -0.511092 -0.5173
## PEXPOSURE     0.047176  0.048032  0.049527  0.049722  0.050477  0.0488
## PCTAGE65P     3.911526  3.933832  3.959192  3.962334  3.979552  3.9509
## PCTOWNHOME   -0.559358 -0.557968 -0.557682 -0.555498 -0.554563 -0.5600

模型结果的数据结构是list，模型的主要结果在下图红框所示的位置：

以截距为例进行可视化：

library(sf)
NY8_sf <- st_as_sf(NY8)
NY8_sf$Intercept <- model$SDF@data$X.Intercept.

plot(NY8_sf["Intercept"])

参考文献：
https://mirrors.tuna.tsinghua.edu.cn/CRAN/web/packages/spgwr/spgwr.pdf