基本地理加权回归模型
对于有 \(n\) 个样本和 \(p\) 个自变量的数据集,样本 \(i\) 对应的基本地理加权回归模型定义为
其中 \(y_i\) 是因变量, \(x_{ik}\) 是第 \(k\) 个自变量, \(\beta_{ik}\) 是第 \(k\) 个回归系数, \(\beta_{i0}\) 是截距,\(\epsilon_i\) 是随机误差且 \(\epsilon_i \sim N(0, \sigma^2)\) , \(\sigma\) 是标准差。 \(\beta_i\) 通过如下方式进行估计
其中 \(\beta_i=(\beta_{i0},\beta_{i1},\cdots,\beta_{ip})'\) 是回归系数向量, \(X\) 是自变量矩阵,\(y\) 是因变量向量, \(W\) 被称为空间权重矩阵,通过如下方式定义
在 \(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 变化阈值,用于确定一个模型与另一个模型相比拟合效果是否有显著改进。一般来说,该值地大小取决于样本地数量。样本数量越大,就需要一个更大的阈值。