地理加权回归模型

〖摘 要〗 空间异质性是地理学第二定律的核心。从地理信息科学角度，空间异质性主要包含两种类型，一是随空间变化，空间某些变量之间的关系发生了明显变化；二是随空间变化，空间某些变量的统计量（如：均值、方差）会出现平稳或者不平稳的变化。地理加权回归是空间计量学、地理空间统计学中为研究第一种空间异质性（即变量间关系的空间异质性）而提出的工具，在多元变量的空间插值或预测等方面具有重要作用。本文为相关原理的基本介绍。

〖原 文〗 Yamagata, Y. and Seya, H. (eds) (2020) Spatial analysis using big data: methods and urban applications. London, United Kingdom ; San Diego, CA: Academic Press, an imprint of Elsevier (Spatial econometrics and spatial statistics). Chapter 6

1 引言

1.1 全局空间最小二乘回归的问题

在地学空间分析中，$n$ 组观测数据通常是在 $n$ 个不同地理位置上获取的样本数据，全局空间回归模型就是假定回归参数与样本数据的地理位置无关，或者说在整个空间研究区域内保持稳定一致，那么在 $n$ 个不同地理位置上获取的样本数据，就等同于在同一地理位置上获取的 $n$ 个样本数据，其回归模型就是最小二乘法回归模型，采用最小二乘估计得到的回归参数户既是该点的最优无偏估计，也是研究区域内所有点上的最优无偏估计。

而在实际问题研究中我们发现：回归参数在不同地理位置上往往表现为不同，也就是说回归参数随地理位置变化而变化，存在空间非平稳性，此时如果仍采用全局空间回归模型，得到的回归参数将是在整个研究区域内的参数平均值，无法反映回归参数的真实空间特征。

1.2 一些分析空间非平稳性的方法

目前，空间非平稳性的分析方法已经有了很多。其中：

针对某些变量之间关系的空间异质性问题，有：

• 广义最小二乘法。该方法把回归模型的参数转化成地理空间位置的线性函数，利用模型参数的线性函数得到模型的参数估计，模型参数随空间位置的变化而变化。该方法在参数变化十分复杂的情况下的应用受到一定限制。

• 空间滤波方法。该方法借鉴时间序列分析中利用“预测–修正”原理来补偿参数漂移的思路，将其应用到空间数据中来分析回归参数的漂移，用迭代计算来估计参数。该方法局限性是其参数估计不能够做统计检验。

• 随机系数模型。假设回归模型的系数为随机变量，相互独立且为有限混合分布，基于贝叶斯原理和最大似然法来估计参数的概率。

• 多层次模型。假设回归模型的系数为随机变量，相互独立且呈正态分布，基于贝叶斯原理和最大似然法来估计参数的概率。

1.3 地理加权回归模型的提出

为解决该问题，国外有学者提出了空间变参数回归模型（Spatially Varying-Coefficient Regression Model） （Fosterand Gorr，1986； Gorrand Olligschlaeger，1994），将数据的空间结构嵌入回归模型中，使回归参数变成观测点地理位置的函数。

Fortheringham 等（Brunsdon 等，1996；Fortheringham 等，1997；Brunsdon 等，1998）在空间变系数回归模型基础上利用局部光滑思想，提出了地理加权回归模型（Geographieally Weighted Regression Model， GWR）。

2 基础模型

地理加权回归模型（GWR）是对全局线性回归模型的扩展。该模型将样本点的地理位置嵌入到回归参数之中，即：

$$\mathbf{y}_{i}=\beta_{0}\left(u_{i}, v_{i}\right)+\sum_{k=1}^{p} \beta_{k}\left(u_{i}, v_{i}\right)\mathrm{x}_{i k}+\varepsilon_{i} \quad i=1,2, \ldots, n$$

式中： $(u_i,v_i)$ 为第 $i$ 个样本点的坐标；$\mathbf{y}_i$ 为第 $i$ 个样本点的因变量观测值；$\mathbf{x}_{ik}$ 为第 $k$ 个自变量，$\beta_k(u_i,v_i)$ 是在第 $i$ 个样本处第 $k$ 个自变量的权重（即模型参数）；$p$ 为自变量个数； $n$ 为样本总数；$\epsilon_i$ 为第 $i$ 个样本点的随机误差。为表述方便，将上式简写为：

$$\mathbf{y}_{\mathrm{i}}=\beta_{0}+\sum_{k=1}^{p} \beta_{i k} \mathbf{x}_{i k}+\varepsilon_{i} \quad i=1,2, \ldots, n$$

该式意味着：对于样本集合中的每个样本点，地理加权回归模型都将利用所有样本为其构建一个线性回归模型，所以 $GWR$ 模型的参数总量将是全局加权回归模型的 $n$ 倍 。若全局加权回归的模型参数为 $(p+1)$ 个，则 $GWR$ 模型的参数为 $n*(p+1)$ 个。

当所有样本点的模型参数一致时（ $\beta_{1k}=\beta_{2k}=\ldots=\beta_{nk}$ ），地理加权回归模型就退化为全局加权回归模型。

3 空间自相关性和异质性的体现

依据地理加权回归模型的基本原理，使用所有样本来为每一个样本点生成一套权重系数，而根据地理学第一定律，距离目标点越远的样本，对目标因变量的影响越小，因此可以在每个样本点的地理加权回归模型参数中嵌入空间结构要素，来体现该空间自相关性。而不同样本点具有不同的回归模型参数，则体现了地理学第二定律的空间异质性原则。

由此可见，地理加权回归模型的关键是如何将空间结构要素嵌入到回归模型的参数中。Fotheringham 等（1996）在传统最小二乘法基础上，通过加入空间权重矩阵 $W$ 的方式，实现空间结构要素的嵌入。每个目标点的 $W$ 不同，体现了空间异质性，而矩阵内权重的大小（根据地理学第一定律，通常受距离影响），体现了自相关程度。

$$\begin{equation} \begin{aligned} 传统最小二乘：\hat\beta(u_{i} , v_{i})&=(X^TX ) ^ { - 1 } X ^ { T } Y\\ \downarrow\\ 地理加权最小二乘：\hat\beta(u_{i} , v_{i})&=(X^TW(u_i,v_i) X ) ^ { - 1 } X ^ { T } W ( u _ { i } , v _ { i } ) Y \end{aligned} \end{equation}$$

式中：

$$X=\left[\begin{array}{cccc}1 & x_{11} & \ldots & X_{1 k} \\ 1 & x_{21} & \ldots & X_{2 k} \\ \ldots & \ldots & \ldots & \ldots \\ 1 & X_{n 1} & \ldots & X_{n k}\end{array}\right]$$

$$W\left(u*{i}, v*{i}\right)=W(i)=\left[\begin{array}{cccc}W_{i 1} & 0 & \ldots & 0 \\ 0 & W_{i 2} & \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots \\ 0 & 0 & \ldots & W_{i n}\end{array}\right]$$

$$\beta=\left[\begin{array}{cccc}\beta_{0}\left(u_{1}, v_{1}\right) & \beta_{1}\left(u_{1}, v_{1}\right) & \ldots & \beta_{p}\left(u_{1}, v_{1}\right) \\ \beta_{0}\left(u_{2}, v_{2}\right) & \beta_{1}\left(u_{2}, v_{2}\right) & \ldots & \beta_{p}\left(u_{2}, v_{2}\right) \\ \ldots & \ldots & \ldots & \ldots \\ \beta_{0}\left(u_{n}, v_{n}\right) & \beta_{1}\left(u_{n}, v_{n}\right) & \ldots & \beta_{p}\left(u_{n}, v_{n}\right)\end{array}\right]$$

$$Y=\left[\begin{array}{c}y_{1} \\ y_{2} \\ \cdots \\ y_{n}\end{array}\right]$$

其中：$\hat{\beta}$ 是模型参数 $β$ 的估计值，$n$ 是空间样本点数，$p$ 是自变量个数，$W_{ik}$ 是对位置 $i$ 刻画模型时赋予第 $k$ 个样本的权重。由此可见，对于地理加权回归模型来说，关键是如何计算空间权重矩阵 $W$ （空间权重矩阵量化描述了 $n$ 个空间位置之间相互的关系，可百度）。

由于地理加权回归模型中的回归参数在各采样点上均不同，因此其未知参数的个数为 $n×（p + 1）$，远大于观测样本个数 $n$ ，采用传统参数估计方法同时计算未知参数比较困难，因此有学者提出了一些方法来拟合该模型。

Foste & Gorr（1986）和 Gorr & Olligsehiaeger（1994）利用广义阻尼负反馈（generalized damped negative feedback）方法估计未知参数在各地理位置的值，该估计方法只是在很直观的意义上考虑数据的空间结构，加之估计方法较为复杂，很难对估计量作深入的统计推断方面的研究。

Brunsdon 等（1996）在局部多项式光滑思想上提出了偏差和方差折衷（Bias-Variance Trade-off）的解决思路：假设回归参数为一连续表面，位置相邻则回归参数非常相似，在估计采样点 $i$ 的回归参数时，以采样点 $i$ 及其邻域采样点上的观测值构成局域子样本，建立全局线性回归模型，然后采用最小二乘方法得到回归参数估计 $\hat\beta_{ik}（k=0，1，2，…，p）$ 。对于另一个采样点 $i+1$ ，采用另一个相应的局域子样来估计，以此类推。由于在回归过程中，由其它样点上的观测值来估计 $i$ 点上的模型参数，因此 $i$ 点上的参数估计不可避免存在偏差，为有偏估计。而且，参与回归估计的子样规模越大，参数估计的偏差就越大，反之，参与回归估计的子样规模越小，参数估计的偏差就越小。从降低偏差角度考虑应当尽量减少子样规模，但子样规模的减少必然导致回归参数估计值的方差增加，精度降低。

4 空间权函数的选择

空间权重矩阵是地理加权回归模型（$GWR$）的核心（Brunsdonetal， 2000），而权重是由空间坐标的某种函数（即空间权函数）形式体现的，因此，如何选取空间权函数对地理加权回归模型（$GWR$）的参数估计影响很大。常见方法有：

（1）距离阈值法

距离阈值法是最简单的空间权函数，它的关键是选取合适的距离阈值 $D$ ，然后与每一个数据点 $j$ 与回归点 $i$ 之间的距离 $d_{ij}$ 进行比较，若 $d_{ij}$ 大于该阈值则权重为 0 ，否则为 1 ，即

$$w_{i j}=\left\{\begin{array}{ll}1 & d_{i j} \leq D \\ 0 & d_{i j}>D\end{array}\right.$$

该权重函数实质就是滑动窗口，计算简单但函数不连续，因此较少采用。

（2）距离反比法

Tobler（1970）地理学第一定律认为空间相近的地物比相远的地物具有更强的相关性，因此在估计回归点 $i$ 的参数时，应对回归点的邻域给予更多的关注。根据该思路，人们自然想到用距离来衡量该空间关系：

$$w_ { i j } = 1 / d _ { i j } ^ { \alpha }$$

这里 $a$ 为合适的常数，当 $a$ 取值为 1 或 2 时，对应的是距离倒数和距离倒数的平方。该方法简洁明了，但对于回归点本身也是样本数据点的情况，就会出现回归点观测值权重无穷大的情况，若从样本中剔除自身又会降低参数估计精度，所以直接采用距离反比法的模型也不多，需要对其修正。

（3）高斯（高斯）函数法

高斯（高斯）函数法就是表示 $w_{ij}$ 与 $d_{ij}$ 之间的连续单调递减函数，可以克服上述空间权函数不连续的缺点。其函数形式如下：

$$w_{i j}=\frac {1}{\exp \left(0.5*\left(d_{i j} / b\right)^{2}\right)}$$

图 1 高斯空间权函数

式中 $b$ 是描述权重与距离之间函数关系的非负衰减参数，称之为带宽（Bandwidth）。带宽越大，权重随距离的增加而衰减越慢，带宽越小，权重随距离增加而衰减越快。

（3） bi-square 函数法

实际应用中往往会将对回归参数估计几乎没有影响的样点截掉，不参与计算，从而形成以有限高斯函数来代替高斯函数的模型，最常采用的便是 bi-square 函数（Bmndonetal，1997； Fotheringham et al， 1998）：

$$w_{i j}=\left\{\begin{array}{cl}{\left[1-\left(d_{i j} / b\right)^{2}\right]^{2}} & d_{i j} \leq b \\ 0 & d_{i j}>b\end{array}\right.$$

图 2 bi-square 空间权函数

从上式可以看出，bi-square 函数法可以看成是距离阈值法和高斯（高斯）函数法的结合。带宽范围内的回归点，可以通过有限高斯函数来计算数据点的权重，而带宽之外的数据点权重为 0。

5 带宽的确定与优化

在实际应用中发现，地理加权回归分析对高斯权函数和 bi-Square 权函数的选择并不是很敏感，但对特定权函数的带宽却很敏感（如图 3），带宽过大回归参数估计的偏差过大，带宽过小又会导致回归参数估计的方差过大。最小二乘平方和是最常采用的优化原则之一，但对于地理加权回归分析中的带宽选择却失去了作用，这是因为对于 $ \sum*{i=1}^{n}\left[y*{i}-\hat{y}{i}(b)\right]^{2}=\min $ 而言，带宽 $b$ 越小，参与回归分析的数据点权重越小，预测值 $\hat y_i(b)$ 越接近实际观测值 $y_i$ ，从而 $ \sum{i=1}^{n}\left[y_{i}-\hat{y}_{i}(b)\right]^{2} \approx 0 $ ，也就是说最优带宽是只包含一个样本点的狭小区域。为此，科学家给出了多种确定或者优化带宽的方法。

图 3 不同权函数与带宽选择对参数估计的影响

（1）交叉验证方法

Cleveland （1979）、Bowman （1984）建议采用用于局域回归分析的交叉验证（cross-validation, CV）方法，其公式表达为：

$$\mathrm{CV}=\sum_{i=1}^{n}\left[y_{i}-\hat{y}_{\neq \mathrm{i}}(b)\right]^2$$

其中，$\hat y_{\neq i}(b)$ 是 $y_i$ 的估计值。该方法在生成 $\hat y$ 时排除了点 $i$ 的观测值。这样当 $b$ 变得很小时，模型仅仅刻画点 $i$ 附近样点而没有包括 $i$ 本身。把不同的带宽 $b$ 和相应的 $CV$ 绘制成趋势线，那么就可以找出 $CV$ 值最小时，对应的最佳带宽 $b$ 。

在实际应用中为减少计算量，Loader 于 1999 年提出了一种近似交叉验证统计量的方法，称为广义交叉验证方法（generalized cross validation，GCV）：

$$G C V=\frac{1}{n} \frac{\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}(b)\right)^{2}}{(1-\operatorname{tr}(\boldsymbol{S}(b) / n))^{2}}=\frac{n \sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}(b)\right)^{2}}{\left(n-\operatorname{tr}(\boldsymbol{S}(b))^{2}\right.}$$

由 $KT$ 变换矩阵 $S$ 的构成可知，当带宽很小时，地理加权回归分析的有效参数个数趋近样本数量 $n$ ，上式中的分母趋于零，这样即便预测值 $\hat y_i(b)$ 趋向 $y_i$ ，GCV 也不会等于 0。

（2） AIC 准则

Akaike 通过对基于极大似然原理的参数估计方法加以修正，提出了一种较为通用的模型选择准则，称为 Akaike 信息量准则（Akaike Information Criterion，AIC）。AIC 定义为（Akaike，1974）：

$$A I C = - 2 \ln L ( \hat\theta _ { L } , x ) + 2 q$$

其中， $\hat \theta_L$ 为 $\theta$ 的极大似然估计，$q$ 为参数的个数。

AIC 准则应用比较广泛，Hurvich et al 将 AIC 准则扩展到非参数回归分析中的光滑参数选择（Hurvich et al， 1998），Brunsdon 和 Fotheringham 则在 Hurvich 等研究基础上将其进一步用于地理加权回归分析中的权函数带宽选择（Brunsdon et al，2002； Fotheringham et al， 2002），其公式为：

$$AIC_c=2 n \ln (\hat{\sigma})+n \ln (2 \pi)+n \frac{n+\operatorname{tr}(S)}{n-2-\operatorname{tr}(S)}$$

其中，下标 $c$ 表示“修正后的” AIC 估计值，$n$ 是样点大小，$\hat \sigma$ 是误差项估计的标准差，$tr(S)$ 是$GWR$ 的 K-T 变换矩阵 S 的迹，它是带宽的函数。AIC 有利于评价 $GWR$ 模型是否比传统最小二乘模型更好地模拟数据。

其简单形式表示为：

$$AIC=2 n \ln (\hat{\sigma})+n \ln (2 \pi)+n+\operatorname{tr}(S)$$

（3） 贝叶斯信息准则

1978 年 SehwartZ 提出了贝叶斯信息准则（Bayesian Information Criterion，BIC），该准则可以使自回归模型的阶数适中，故常被用来确定回归模型中的最优阶数，2002 年 Nakaya 将其用于地理加权回归分析中的权函数带宽选择。BIC 准则与 AIC 准则非常相似，只是惩罚因子不同，其公式为

$$BIC= -2 \ln L\left(\hat{\theta}_{L}, x\right)+q \ln n$$

式中 $\hat \theta_L$ 为 $\theta$ 的极大似然估计，$q$ 为参数的个数，$n$ 为样本个数，使 BIC 最小的模型为“最优”模型。式中可以看出，BIC 准则对于具有相同未知参数个数的模型，样本数越多，惩罚度越大，对于具有相同样本的情况，则趋于选择具有更少参数的模型为最优。

与 AIC 不同的是，BIC 准则要求模型为 Bayesian 模型，即每个候选模型都必须具有相同的先验概率，而实际上模型参数的先验分布通常是不知道的，另外如何将 BIC 准则扩展到可变带宽的非参数模型，用有效参数个数来代替全局参数个数还不是很清楚。

6 空间非平稳性的显著性检验

7 大 N 问题

有一些 GWR 方法适用于大样本。

一种方法是使用样本进行基于 CV 的带宽优化，这是 GWR 中最重的部分。例如，Yu（2007）在确定样本量的 68606 个样本中随机抽取 3437 个样本，使估计误差在 2% 的范围内。 Feuillet 等 (2018) 将他们的研究区域分成更小的单元，并在每个单元中应用 GWR。

Sachdeva 2022 年的 《采用 “空间变化过程建模” 还是 “非线性建模” ？》: 在各种类型的局部统计模型的校准中获得的空间变化参数估计的研究是司空见惯的。这种估计的变化通常用空间变化过程来解释。本文强调，在将这种变化与空间变化过程相关联之前，应先检查非线性方面的空间变化参数估计的另一种解释。这可以通过描述和演示的简单筛选程序来实现，并且可以轻松应用于任何局部模型的结果。突出显示问题并展示解决方案，使用一组模拟数据，然后使用真实世界的数据集。该论文还强调了相反的情况，即当实际关系是线性但空间变化时，GAM 的不当应用会产生虚假的非线性结果。

8 应用案例

9 总结

参考文献

• [1] Brunsdon, C., Fotheringham, S., Charlton, M., 1998. Geographically weighted regression. Journal of the Royal Statistical Society: Series D (The Statistician) 47 (3), 431e443.
• [2] Brunsdon, C., Charlton, M., Harris, P., 2012. Living with collinearity in local regression models. In: 10th International Symposium on Spatial Accuracy Assessment in Natural Resources and Environmental Sciences, Florianopolis, Brazil.
• [3] da Silva, A.R., Fotheringham, A.S., 2016. The multiple testing issue in geographically weighted regression. Geographical Analysis 48 (3), 233e247.
• [4] Feuillet, T., Commenges, H., Menai, M., Salze, P., Perchoux, C., Reuillon, R., Kesse-Guyot, E., Enaux, C., Nazare, J.-A., Hercberg, S., Simon, C., Charreire, H., Oppert, J.M., 2018. A massive geographically weighted regression model of walking-environment relationships. Journal of transport geography 68, 118e129.
• [5] Fotheringham, A.S., Brunsdon, C., Charlton, M., 2002. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Wiley, New York.
• [6] Fotheringham, A.S., Yang, W., Kang, W., 2017. Multiscale geographically weighted regression (mgwr). Annals of the Association of American Geographers 107 (6), 1247e1265.
• [7] Getis, A., 2008. A history of the concept of spatial autocorrelation: A geographer’s perspective. Geographical Analysis 40 (3), 297e309.
• [8] Getis, A., Griffith, D.A., 2002. Comparative spatial filtering in regression analysis. Geographical Analysis 34 (2), 130e140.
• [9] Harris, R., Singleton, A., Grose, D., Brunsdon, C., Longley, P., 2010. Grid-enabling geographically weighted regression: a case study of participation in higher education in England. Transactions in GIS 14 (1), 43e61.
• [10] Leong, Y.Y., Yue, J.C., 2017. A modification to geographically weighted regression. International Journal of Health Geographics 16 (1), 11. Li, Z., Fotheringham, A.S.,
• [11] Li, W., Oshan, T., 2019. Fast Geographically Weighted Regression (FastGWR): a scalable algorithm to investigate spatial process heterogeneity in millions of observations. International Journal of Geographical. Information Science 33 (1), 155e175.
• [12] Lu, B., Brunsdon, C., Charlton, M., Harris, P., 2017. Geographically weighted regression with parameter-specific distance metrics. International Journal of Geographical Information Science 31 (5), 982e998.
• [13] Lu, B., Yang, W., Ge, Y., Harris, P., 2018. Improvements to the calibration of a geographically weighted regression with parameter-specific distance metrics and bandwidths. Computers, Environment and Urban Systems 71, 41e57.
• [14] Mei, C.L., He, S.Y., Fang, K.T., 2004. A note on the mixed geographically weighted regression model. Journal of Regional Science 44 (1), 143e157.
• [15] Nakaya, T., Fotheringham, A.S., Brunsdon, C., Charlton, M., 2005. Geographically weighted Poisson regression for disease association mapping. Statistics in Medicine 24 (17), 2695e2717.
• [16] Tran, H.T., Nguyen, H.T., Tran, V.T., 2016, October. Large-scale geographically weighted regression on Spark. In: 2016 Eighth International Conference on Knowledge and Systems Engineering (KSE), pp. 127e132. IEEE.
• [17] Wheeler, D., Tiefelsdorf, M., 2005. Multicollinearity and correlation among local regression coefficients in geographically weighted regression. Journal of Geographical Systems 7 (2), 161e187.
• [18] Yang, W., 2014. An Extension of Geographically Weighted Regression with Flexible Bandwidths. Doctoral dissertation. University of St Andrews.
• [19] Yu, D., 2007. Modeling owner-occupied single-family house values in the city of Milwaukee: a geographically weighted regression approach. GIScience & Remote Sensing 44 (3), 267e282.