基本地理加权回归模型

对于有 \(n\) 个样本和 \(p\) 个自变量的数据集,样本 \(i\) 对应的基本地理加权回归模型定义为

\[y_i = \beta_{i0} + \sum_{k=1}^p \beta_{ik}x_{ik} + \epsilon_{i}\]

其中 \(y_i\) 是因变量, \(x_{ik}\) 是第 \(k\) 个自变量, \(\beta_{ik}\) 是第 \(k\) 个回归系数, \(\beta_{i0}\) 是截距,\(\epsilon_i\) 是随机误差且 \(\epsilon_i \sim N(0, \sigma^2)\)\(\sigma\) 是标准差。 \(\beta_i\) 通过如下方式进行估计

(1)\[\beta_i = \left(X' W X \right)^{-1} X' W y\]

其中 \(\beta_i=(\beta_{i0},\beta_{i1},\cdots,\beta_{ip})'\) 是回归系数向量, \(X\) 是自变量矩阵,\(y\) 是因变量向量, \(W\) 被称为空间权重矩阵,通过如下方式定义

\[\begin{split}W = \begin{pmatrix} w_{i1} & & & \\ & w_{i2} & & \\ & & \ddots & \\ & & & w_{in} \\ \end{pmatrix}\end{split}\]

\(W\) 中的每个 \(w_{ij}\) 都通过一个核函数 \(k\) 根据样本 \(i\) 到样本 \(j\) 的距离得到。距离越远,权重越小。

核函数

下面是一些常用的核函数:

Gaussian
\[k(d;b) = \exp\left\{- \frac{d^2}{2 b^2}\right\}\]
Exponential
\[k(d;b) = \exp\left\{- \frac{|d|}{b}\right\}\]
Bi-squared
\[\begin{split}k(d;b) = \left\{ \begin{array}{ll} \left[ 1 - \left( \frac{d}{b} \right)^2 \right]^2, & \mbox{if } d < b \\ 0, & \mbox{otherwise} \end{array} \right.\end{split}\]
Tri-cube
\[\begin{split}k(d;b) = \left\{ \begin{array}{ll} \left[ 1 - \left( \frac{d}{b} \right)^3 \right]^3, & \mbox{if } d < b \\ 0, & \mbox{otherwise} \end{array} \right.\end{split}\]
Box-car
\[\begin{split}k(d;b) = \left\{ \begin{array}{ll} 1, & \mbox{if } d < b \\ 0, & \mbox{otherwise} \end{array} \right.\end{split}\]

参数 \(b\) 被称为“带宽”。 其值通常 通过黄金分割算法从数据中根据某些指标进行优选。支持下列指标,

十字交叉验证(CV)

对于给定的带宽 \(b\) 十字交叉验证值通过如下方式定义

\[CV(b) = \sum_{i=1}^n \left( y - x_i \hat{\beta}_{-i} \right)^2\]

其中 \(x_i\) 是矩阵 \(X\) 的第 \(i\) 行,\(\hat{\beta}_{-i}\) 是刨除第 \(i\) 个样本的回归系数估计值。它同样也是根据 (1) 进行估计,但是 \(w_{ii} = 0\)

校正的赤池信息准则(AIC:sub:c

对于给定的带宽 \(b\),AIC 值通过如下方式定义

\[AIC(b) = 2n \ln \hat{\sigma} + n \ln 2pi + n \left\{ \frac{n+tr(S)}{n - 2 - tr(S)} \right\}\]

其中 \(\hat{\sigma}\) 是随机误差标准差的估计值,\(S\) 被称为帽子矩阵 “hat matrix” ,定义为

\[\begin{split}S = \begin{pmatrix} x_1 (X'W_1X)^{-1}X'W_1 \\ x_2 (X'W_2X)^{-1}X'W_2 \\ \vdots \\ x_n (X'W_nX)^{-1}X'W_n \end{pmatrix}\end{split}\]

它的效果是

\[\hat{y} = Sy\]

距离度量

除了欧氏距离,还有其他很多距离度量可以被用于GWR。目前支持如下两种距离度量。

坐标系距离

根据其空间参考(CRS)类型计算直线距离。 如果是投影坐标系,对于在位置 \((u_i,v_i)\)\((u_j,v_j)\) 处的两个样本,

\[d_{ij} = \sqrt{ (u_i - u_j)^2 + (v_i - v_j)^2 }\]

如果是地理坐标系,计算他们之间的大圆距离(测地线距离)。

Minkwoski 距离

该距离度量仅可应用于投影坐标系。定义为

\[d_{ij} = \sqrt[p]{ |u_i - u_j|^p + |v_i - v_j|^p }\]

未来,我们将支持通过距离矩阵设置距离。

案例

如果要拟合一个基础GWR模型,使用 gwm::GWRBasic

基本用法

#include <armadillo>
using namespace arma;

mat coords = randr(100, 2, distr_param(0, 25));
mat x = join_rows(ones(100, 1), randn(100, 2));
mat beta = join_rows(
    ones(100) * 3.0,
    1.0 + (coords.col(0) + coords.col(1)) / 12.0,
    1.0 + (36.0 - (6.0 - coords.col(0) / 2)) % (36.0 - (6.0 - coords.col(1) / 2)) / 324
);
vec y = sum(x % beta, 1);

CGwmCRSDistance distance(false);
CGwmBandwidthWeight bandwidth(25, true, CGwmBandwidthWeight::Gaussian);
CGwmSpatialWeight spatial(&bandwidth, &distance);

GWRBasic algorithm;
algorithm.setCoords(coords);
algorithm.setDependentVariable(y);
algorithm.setIndependentVariables(x);
algorithm.setSpatialWeight(spatial);
mat beta_hat = algorithm.fit();

带宽优选

如果不确定带宽值,可以按照如下方式让算法自动优选带宽:

GWRBasic algorithm;
algorithm.setCoords(coords);
algorithm.setDependentVariable(y);
algorithm.setIndependentVariables(x);
algorithm.setSpatialWeight(spatial);
algorithm.setIsAutoselectBandwidth(true);
algorithm.setBandwidthSelectionCriterion(GWRBasic::BandwidthSelectionCriterionType::AIC);
mat beta_hat = algorithm.fit();

传递给 gwm::GWRBasic::setBandwidthSelectionCriterion() 的参数可以是 gwm::GWRBasic::BandwidthSelectionCriterionType 中的任意值。

自变量优选

如果不希望所有自变量都被纳入模型,而是仅仅纳入一些显著的变量,可以通过下列方式零算法自动优选自变量:

GWRBasic algorithm;
algorithm.setCoords(coords);
algorithm.setDependentVariable(y);
algorithm.setIndependentVariables(x);
algorithm.setSpatialWeight(spatial);
algorithm.setIsAutoselectIndepVars(true);
algorithm.setIndepVarSelectionThreshold(3.0);
mat beta_hat = algorithm.fit();

传递给 gwm::GWRBasic::setIndepVarSelectionThreshold() 的参数是一个 AIC 变化阈值,用于确定一个模型与另一个模型相比拟合效果是否有显著改进。一般来说,该值地大小取决于样本地数量。样本数量越大,就需要一个更大的阈值。