GWmodel | 地理加权模型（Ⅱ-1）：地理加权主成分分析（GWPCA）

地理加权回归（GWR）相比于普通的回归模型能够考虑到回归系数的空间异质性。但实际上，地理加权并非GWR模型所独有，其他数据分析方法同样也能通过加入地理权重进行改进。本篇介绍的就是地理加权主成分分析（GWPCA）。

本部分使用的工具包主要是GWmodel工具包，开发者是武汉大学的卢宾宾。它相比于spgrw工具包，功能更加丰富，除了GWR外，还能用来运行GWPCA等多种地理加权模型。

本篇包含以下内容：

使用GWmodel工具包中的函数运行GWR模型；
GWPCA模型的理论基础；
使用GWmodel工具包中的函数运行GWPCA模型。

1 GWR模型

在第Ⅰ部分已经介绍了GWR模型的主要构成部分以及spgwr工具包中的相关函数，这里再介绍下GWmodel工具包中的相关函数。

函数	线性GWR	广义线性GWR
带宽函数	`bw.gwr`	`bw.ggwr`
模型函数	`gwr.basic`	`ggwr.basic`

以下仅以线性GWR模型为例。

工具包和原始的示例数据加载如下：

library(GWmodel)
data(Georgia)
data(GeorgiaCounties)

上述代码运行后可以得到Gedu.counties和Gedu.df两个数据，数据结构分别为SpatialPolygonsDataFrame和data.frame。通过以下代码将二者合并：

library(sf)
library(sp)
library(tidyverse)
Gedu.df$AreaKey <- factor(Gedu.df$AreaKey)
dta <- Gedu.counties %>%st_as_sf() %>%left_join(Gedu.df, by = c("AREAKEY" = "AreaKey")) %>%as("Spatial")

得到的dta是一个含属性的sp格式面状矢量数据。至此，示例数据的准备工作已经完成。

下面看下主要函数的语法结构。

带宽函数的语法结构：

bw.gwr(formula, data,approach = "CV",kernel = "bisquare",adaptive = FALSE,p = 2, theta = 0,longlat = F, dMat,parallel.method = F, parallel.arg = NULL)

GWmodel::bw.gwr()函数和spgwr::gwr.sel()函数中的参数基本上是对应的：

approach——method：确定“最优”带宽的方法；可选项有交叉验证（CV）和赤池信息准则（AIC）；

kernel——gweight：距离加权函数。GWmodel工具包提供了5种选择：gaussian、exponential、bisquare、tricube、boxcar，具体见函数说明文档；

adaptive——adapt：FALSE表示带宽为固定的距离；TRUE表示权重不为0的邻近单元的个数固定（k-邻近法），此时带宽是变化的；

longlat——longlat：当空间数据的坐标为经纬度格式时，需将改参数设置为TRUE；

dMat：距离矩阵；可忽略。

模型函数的语法结构：

gwr.basic(formula, data,regression.points,bw, kernel = "bisquare",adaptive = FALSE,p = 2, theta = 0,longlat = F, dMat,F123.test = F,cv = F, W.vect = NULL,parallel.method = FALSE, parallel.arg = NULL)

以下为示例代码。

定义模型表达式和带宽：

formula <- PctPov ~ PctBlack + PctBach
bw.1 <- bw.gwr(formula, data = dta,approach = "aic",adaptive = T,kernel = "gaussian")# “最优”带宽寻找过程
## Adaptive bandwidth (number of nearest neighbours): 105 AICc value: 910.8443
## Adaptive bandwidth (number of nearest neighbours): 73 AICc value: 902.3318
## Adaptive bandwidth (number of nearest neighbours): 51 AICc value: 891.5931
## Adaptive bandwidth (number of nearest neighbours): 40 AICc value: 884.2077
## Adaptive bandwidth (number of nearest neighbours): 30 AICc value: 875.4572
## Adaptive bandwidth (number of nearest neighbours): 27 AICc value: 873.3422
## Adaptive bandwidth (number of nearest neighbours): 22 AICc value: 867.0877
## Adaptive bandwidth (number of nearest neighbours): 22 AICc value: 867.0877

运行GWR模型：

model.1 <- gwr.basic(formula, data = dta,bw = bw.1,adaptive = T,kernel = "gaussian")

GWmodel和spgwr这两个工具包运行GWR模型的方法是类似的，有兴趣的读者可以尝试使用GWmodel工具包中的函数运行广义线性GWR模型。

2 GWPCA模型

2.1 理论基础

关于主成分分析（PCA）详见推文：数据降维之主成分分析（PCA）

地理加权主成分分析（GWPCA）与普通主成分分析的区别在于计算协方差矩阵时加入了地理权重矩阵。

记为行、列的矩阵，其中行表示样本，列表示变量，再将进行列中心化得到矩阵。

PCA计算协方差矩阵：

GWPCA计算协方差矩阵：

其余步骤相同。

2.2 相关函数

带宽函数为bw.gwpca()，模型函数为gwpca()。

gwpca()函数的语法结构如下：

gwpca(data, elocat,vars, k = 2,robust = FALSE,kernel = "bisquare",adaptive = FALSE,bw, p = 2, theta = 0,longlat = F, cv = T,scores=F, dMat)

vars：参与主成分分析的变量名；

k：保留的主成分个数。

# 变量名
varnames <- colnames(select(dta@data, starts_with("Pct")))
varnames
## [1] "PctRural" "PctBach"  "PctEld"   "PctFB"    "PctPov"   "PctBlack"# 带宽
bw.2 <- bw.gwpca(data = dta, vars = varnames, k = 3)
bw.2
## [1] 24450.79# 模型
model.2 <- gwpca(data = dta, vars = varnames,bw = bw.2, k = 3)

参考资料：

[1] Harris P, Clarke A, Juggins S, Brunsdon C, Charlton M (2015). Enhancements to a geographically weighted principal components analysis in the context of an application to an environmental data set. Geographical Analysis 47: 146-172