空间回归模型综述

【摘要】 空间回归模型是按照空间区位研究变量之间关系的主要数学工具。根据回归模型是否存在空间同质特征（或反之空间异质性），通常可以将空间回归模型划分为 全局空间回归模型 和 局部空间回归模型。 【原文】 自编 【作者】 濮国梁，北京大学

空间回归模型是按照空间区位研究变量之间关系的主要数学工具之一。根据空间回归模型是否同质（或反之是否异质，可以简单理解为模型参数是否会随空间位置变化而变化 ），可以将空间回归模型划分为 全局空间回归模型 和 局部空间回归模型。其中：

全局空间回归模型以空间依赖性研究为主体，主要探究的是不同变量、误差项之间的空间交互效应；

局部回归模型则相对复杂，它不仅要研究变量、误差项之间的空间交互效应，还要研究模型本身（ 通常指 模型结构 和 模型参数 ）的空间变化规律，探究的重点是 空间异质性。

空间效应 (spatial effects) 包括 空间依赖性 与 空间异质性。而全局空间回归重点考虑的是 空间依赖性，主要表现为不同变量、误差项之间的各种 空间交互效应。

以标准线性回归模型为起点，全局空间回归统计模型中主要涉及三种交互效应：

内生交互效应：不同空间单位的响应变量（ YYY ）之间存在的交互效应，即 “空间单位 AAA 的响应变量 yAy_AyA​ 与空间单位 BBB 的响应变量 yBy_ByB​ 之间的空间交互关系”。 例如：在一个行政区的税收与相邻行政区的税收之间存在交互效应。

外生交互效应：一个空间单位的响应变量（ YYY ）与其他空间单位的解释变量（ XXX ）之间的交互效应，即 “空间单位 BBB 的解释变量 xBx_BxB​ 与空间单位 AAA 的响应变量 yAy_AyA​ 之间的空间交互关系”。例如：一个国家的经济增长不仅取决于其内部收入水平、储蓄率、人口增长、技术变革等因素，还与相邻国家的这些因素有关。

误差项交互效应 ：不同空间单位的误差项（ uuu ）之间的交互效应，即 “空间单位 AAA 的误差项 uAu_AuA​ 与空间单位 BBB 的误差项 uBu_BuB​ 之间的空间交互关系”，例如：行政指令对不可预测的财政政策变化的校正机制。

为了解释上述空间交互效应，历史上提出了大量模型来对这些空间交互效应进行建模，而其中线性模型的研究最为普遍和深入。

空间计量模型最初关注的基本焦点是 空间滞后模型（也称为空间自回归模型，Spatial Auto Regression Model，SAR) 以及 空间误差项模型 (Spatial Error Model，SEM) ，分别对应 内生交互效应 和 误差项交互效应 。Anselin (1988，1996) 提出了基于稳健拉格朗日乘子的空间滞后模型检验和空间误差模型检验程序。

在 200720072007 年，在首届空间计量经济学协会世界大会上，研究人员对于包含多种空间交互效应的模型，研究兴趣日益增加。Harry Kelejian (1998) 提倡使用同时包含 内生交互效应 和 误差项交互效应 的模型。 LeSage and Pace (2009) 将其称为 SAC 模型，尽管没有指出首字母缩略语代表什么。 Elhorst (2010) 则称这种模型为 Kelejian-Prucha 模型。

同样在 200720072007 年的第 545454 届地区科学协会国际会议 （北美洲） 上， LeSage 提倡使用同时包括 内生交互效应 和 外生交互效应 的模型，并称其为 空间杜宾模型（Spatial Durbin Model, SDM）。但 Gibbons 和 Overman 等（2012） 对 空间滞后模型 （SAR）、 空间误差项模型（ SEM ） 和 空间杜宾模型（ SDM ） 模型进行了批判，提倡使用仅含 外生交互效应 的 SLX 模型。

具有所有类型交互效应的完整模型形式如下：

Y=α⋅lN⏟截距项+X⋅β⏟解释变量+δ⋅W1⋅Y⏟内生交互（空间滞后）项+W2⋅X⋅θ⏟外生交互（空间杜宾）项+u⏟误差项u=λ⋅M⋅u+ε\begin{array}{l} Y=\underbrace{\alpha \cdot l_N}_{截距项} + \underbrace{X \cdot \beta}_{解释变量} +\underbrace{\delta \cdot W_1 \cdot Y}_{内生交互（空间滞后）项} + \underbrace{W_2 \cdot X \cdot \theta}_{外生交互（空间杜宾）项} + \underbrace{u}_{误差项} \\ u=\lambda \cdot M \cdot u + \varepsilon \end{array} Y=截距项α⋅lN​​​+解释变量X⋅β​​+内生交互（空间滞后）项δ⋅W1​⋅Y​​+外生交互（空间杜宾）项W2​⋅X⋅θ​​+误差项u​​u=λ⋅M⋅u+ε​

该模型被称为 “广义嵌套空间模型”，它包括了上述所有类型的空间交互效应，其中：

其中：

W1,W2,MW_1, W_2, MW1​,W2​,M 为空间权重矩阵，是大小为 N×NN \times NN×N 的非负矩阵， NNN 为空间单位的数量，三者可以相同也可以不同。

YYY 为所有空间单位的响应变量观测值构成的向量，是 N×1N \times 1N×1 的列向量。

XXX 为解释变量矩阵，由所有空间单位内解释变量的观测值组成，大小为 N×KN \times KN×K ，KKK 为解释变量数目。

α⋅lN\alpha \cdot l_Nα⋅lN​ 为截距项，其中 lNl_NlN​ 是大小为 N×1N \times 1N×1 的列向量，α\alphaα 是大小为 N×NN \times NN×N 的常系数对角矩阵。

δ\deltaδ 为所有空间单位之间待估计的空间滞后系数，为 N×NN \times NN×N 的矩阵。

β\betaβ 为所有空间单位内待估计的解释变量系数， 是 K×1K \times 1K×1 的列向量。

θ\thetaθ 为所有空间单位之间待估计的外生交互系数，是 K×1K \times 1K×1 的列向量。

λ\lambdaλ 是空间误差的自相关系数，为 N×NN \times NN×N 的矩阵。

ϵ=(ϵ1，…，ϵN)T\epsilon =(\epsilon_1，…， \epsilon_N)^Tϵ=(ϵ1​，…，ϵN​)T 为 N×1N \times 1N×1 的干扰项。如果假设 ϵi\epsilon_iϵi​ 为独立同分布（对于所有 iii ），则其均值为零且方差为 σ2σ^2σ2 。

在上述统一模型框架下，下图总结了八种常见的线性回归模型图谱（假设 W1=W2=M=WW_1 =W_2 = M = WW1​=W2​=M=W ）。其中图右侧为 普通线性回归模型（ OLS ），不含任何交互效应；图左侧是为 广义嵌套空间模型（ GNS ），包含所有空间交互效应。GNS 右侧的每一个模型，都可以通过对 GNS 模型的一个或多个参数施加限制而得到，图中箭头旁边的标注表明了限制因素。读者也可以从右往左看，以了解在传统线性回归模型基础上，通过追加了哪些交互效应生成了哪些新模型。

需要注意的是：

该图同时表明，有很多可能的模型在计量经济学的理论与实证研究中使用较少。如：同时包含 外生交互效应 和 误差项交互效应 ，但不含 内生交互效应 的 空间杜宾误差模型（SDEM）。

理论研究者和实际操作者对不同类型交互效应的兴趣水平有很大差异。理论研究者最感兴趣的是 空间滞后模型（ SAR ） 和 空间误差模型（ SEM ），以及同时包含 内生交互效应 和 误差项交互效应 的 SAC 模型，因为计量经济学问题常伴随着对这些模型的估计方法。

注: GNS = general nesting spatial model（广义空间嵌套模型）, SAC = spatial autoregressive combined model（空间自回归组合模型）, SDM = spatial Durbin model（空间杜宾模型）, SDEM = spatial Durbin error model（空间杜宾误差模型）, SAR = spatial autoregressive model（空间自回归模型）, SLX = spatial lag of X model（外生空间滞后模型）, SEM = spatial error model（空间误差模型）, OLS = ordinary least squares model（普通最小二乘模型）.

在实证研究中，空间权重矩阵非常重要，起到了量化空间交互效应程度的作用。空间权重矩阵 WWW 是一个非负的已知常数矩阵，通常是对称矩阵，但也能是非对称矩阵。

经常采用的空间权重矩阵包括：

（1）ppp 阶邻接矩阵（如果 p=1p=1p=1 ，则只包含邻居；如果p=2p=2p=2 , 则不仅包含邻居，还包含邻居的邻居，以此类推），有时用 W1W^1W1 和 W2W^2W2 分别表示一阶和二阶；

（2）逆距离矩阵（带有或者不带有分界点）

（3）最近邻矩阵（需指定最近邻个数 qqq 作为超参数）

（4）分块对角矩阵，其中每一个分块代表一组空间单位，这些空间单位之间存在交互效应，但与其他分块内的空间单位之间不存在交互关系。

为了进行更简单的解释，常见做法是对 WWW 进行归一化，使其每一行元素之和等于 111。WWW 非负确保了所有权重介于 000 到 111 之间，且可以将加权操作解释为对其邻近空间单位取值的加权平均。经过行归一化的 WWW 被称为行随机矩阵。

每行元素之和均为 111，这意味着区域 iii 的所有邻居对 iii 的影响之和，等于区域 jjj 的邻居对 jjj 的影响之和。

此假定有时可能过强，Keleji 和 Prucha (2010) 证明：与全矩阵归一化不同，按行归一化很可能导致模型的错误设计，尤其是使用 逆距离矩阵 时，此问题更容易发生。此外，行归一化后的空间权重矩阵一般不再是对称矩阵，这也是其缺陷之一。

也可以按列进行归一化，即每列元素的和等于 111 。空间权重矩阵的列元素表示一个特定空间单位对所有空间的影响，而空间权重矩阵的行元素代表所有其他空间单位对某个特定空间单位的影响。因此，行归一化的效应是把其他空间单位对当前单元的影响进行归一化，而列归一化的效应是把当前空间单位对其他空间单位的影响进行归一化。

（1）模型的定义

对于样本容量为 NNN 的空间序列 $${ y_i }_{i=1}^N$$ ，即使假设空间自相关形式是线性的，待估计参数在理论上最多也可达 (n2−n)(n^2 - n)(n2−n) 个（每个区域最多可受 n−1n-1n−1 个区域影响，共有 nnn 个区域），大大超出样本容量。因为，必须对空间依赖性做出某种共同模式的假设，才能简化并计算参数。

空间滞后的概念引申自时间序列数据的时间滞后，即假设当前时刻的数据线性滞后于前一时刻的数据。不过与时间滞后不同的是，空间滞后更为复杂，因为空间滞后可以来自不同的方向，而且是双向的。

Y=δWY+ε\boldsymbol{Y}=\delta \boldsymbol{W} \boldsymbol{Y}+\boldsymbol{\varepsilon} Y=δWY+ε

其中，WWW 为已知的空间权重矩阵（ 一般是一阶空间邻接矩阵），WYWYWY 是对 YYY 的时间滞后，δ\deltaδ 是空间滞后系数。

该模型中，每个空间单位的空间依赖性仅来自于其空间滞后，因此被称为 空间滞后模型（Spatial Lag Model，SLM）。 相邻地区的响应变量可能相互依赖， 并最终形成一个均衡的结果。 例如，假设以地区税收为响应变量，则不同地区政府出于相互竞争或博弈的考虑，在制定本地区税收时，会考虑周边地区的税收水平。空间依赖性导致变量 $${y_i}_{i=1}^N$$ 之间互相影响，产生内生性。

但是更一般地，会在方程中加入解释变量：

Y=δWY+Xβ+ε\boldsymbol{Y}=\delta \boldsymbol{W} \boldsymbol{Y}+ \boldsymbol{X}\beta+\boldsymbol{\varepsilon} Y=δWY+Xβ+ε

其中，XXX 为 N×KN\times KN×K 数据矩阵，包括 KKK 列解释变量；βk×1β_{k\times1}βk×1​ 为对应的解释系数。

此方程也称为 空间自回归模型（Spatial Auto Regression，SAR 模型）。

如果 δ=0\delta=0δ=0 ，则 SAR 简化为 普通线性回归模型 OLS 。

（2）模型的估计

对于空间自回归模型，常使用 MLE 估计。

首先，假设扰动项 ε∼N(0,σ2I)\boldsymbol{\varepsilon} \sim N\left(\boldsymbol{0}, \sigma^{2} \boldsymbol{I}\right)ε∼N(0,σ2I)

其次，方程可写为$ \boldsymbol{A y} \equiv(\boldsymbol{I}-\lambda \boldsymbol{W}) \boldsymbol{y}=\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{\varepsilon} $

其中， $ \boldsymbol{A} \equiv(\boldsymbol{I}-\lambda \boldsymbol{W}) $ 。

由于雅可比行列式 $ J \equiv\left|\frac{\partial \varepsilon}{\partial \boldsymbol{y}}\right|=\left|\frac{\partial(\boldsymbol{A} \boldsymbol{y}-\boldsymbol{X} \boldsymbol{\beta})}{\partial \boldsymbol{y}}\right|=\left|\frac{\partial \boldsymbol{A} \boldsymbol{y}}{\partial \boldsymbol{y}}\right|=\left|\boldsymbol{A}^{\prime}\right|=|\boldsymbol{A}| $, 根据多维正态密度公式，可写出样本的似然函数：

L(y∣λ,σ2,β)=(2πσ2)−n/2(abs⁡∣A∣)exp⁡{−12σ2(Ay−Xβ)′(Ay−Xβ)}L\left(\boldsymbol{y} \mid \lambda, \sigma^{2}, \boldsymbol{\beta}\right)=\left(2 \pi \sigma^{2}\right)^{-n / 2}(\operatorname{abs}|\boldsymbol{A}|) \exp \left\{-\frac{1}{2 \sigma^{2}}(\boldsymbol{A} \boldsymbol{y}-\boldsymbol{X} \boldsymbol{\beta})^{\prime}(\boldsymbol{A} \boldsymbol{y}-\boldsymbol{X} \boldsymbol{\beta})\right\} L(y∣λ,σ2,β)=(2πσ2)−n/2(abs∣A∣)exp{−2σ21​(Ay−Xβ)′(Ay−Xβ)}

其中，$ \mathrm{abs}|\boldsymbol{A}| $ 表示行列式 $ |\boldsymbol{A}| $ 的绝对值。

对数似然函数为：

ln⁡L(y∣λ,σ2,β)=−n2ln⁡2π−n2ln⁡σ2+ln⁡(abs∣A∣)−12σ2(Ay−Xβ)′(Ay−Xβ)\ln L\left(\boldsymbol{y} \mid \lambda, \sigma^{2}, \boldsymbol{\beta}\right)=-\frac{n}{2} \ln 2 \pi-\frac{n}{2} \ln \sigma^{2}+\ln (\mathrm{abs}|\boldsymbol{A}|)-\frac{1}{2 \sigma^{2}}(\boldsymbol{A} \boldsymbol{y}-\boldsymbol{X} \boldsymbol{\beta})^{\prime}(\boldsymbol{A} \boldsymbol{y}-\boldsymbol{X} \boldsymbol{\beta}) lnL(y∣λ,σ2,β)=−2n​ln2π−2n​lnσ2+ln(abs∣A∣)−2σ21​(Ay−Xβ)′(Ay−Xβ)

（1）模型定义

空间外生滞后模型（ Spatial Lag of X model, SLX ）假设某空间单位的响应变量 yiy_iyi​ 依赖于其相邻空间单位的解释变量。例如：空间单位 iii 的犯罪率不仅依赖于本单元警力，还可能依赖于相邻空间单位的警力。

此时的模型为：

y=Xβ+WXθ+ε\boldsymbol{y}=\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{W} \boldsymbol{X} \boldsymbol{\theta}+\boldsymbol{\varepsilon} y=Xβ+WXθ+ε

其中，WXθ\boldsymbol{W} \boldsymbol{X} \boldsymbol{\theta}WXθ 表示来自相邻单元解释变量的影响，而 θ\boldsymbol{\theta}θ 为相应的系数向量。

此模型被称为 空间杜宾模型（Spatial Durbin Model，SDM）。

（2）模型的估计

此方程不存在内生性，故可直接进行普通最小二乘估计。

解释变量 XXX 与 WXWXWX 之间可能存在多重共线性。如果 θ=0\theta=0θ=0，则简化为普通线性回归模型。

将空间杜宾模型与空间自回归模型相结合，可以得到：

Y=δWY+Xβ+WXθ+ε\boldsymbol{Y}=\delta \boldsymbol{W} \boldsymbol{Y}+\boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{W} \boldsymbol{X} \boldsymbol{\theta} + \boldsymbol{\varepsilon} Y=δWY+Xβ+WXθ+ε

（1）模型的定义

空间依赖性还可能通过误差项来体现。考虑以下 空间误差项模型(Spatial Errors Model，简记 SEM)：

y=Xβ+u\boldsymbol{y}=\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{u} y=Xβ+u

其中，扰动项 uuu 的生成过程为 ：

u=λMu+ε,  ε∼N(0,σ2I)\boldsymbol{u}=\lambda \boldsymbol{M} \boldsymbol{u}+\boldsymbol{\varepsilon},\,\, \boldsymbol{\varepsilon} \sim N(0,\sigma^2I) u=λMu+ε,ε∼N(0,σ2I)

扰动项 u\boldsymbol{u}u 之间存在空间依赖性，但不依赖于 X\boldsymbol{X}X ，而且该扰动项与某些对 Y\boldsymbol{Y}Y 有影响的未知变量存在空间相关性。

如果 λ=0\lambda=0λ=0 ，则简化为普通线性回归模型。

（2）模型的估计

对于 SEM 模型，尽管扰动项存在自相关，OLS 仍一致，但缺乏效率。最有效率的方法为 MLE。

（1）模型的定义

（Spatial Autoregressive Combined model，SAC）更一般的空间计量模型将空间自回归模型 ( SAR ) 与空间误差项模型 ( SEM ) 结合起来：

Y=δWY+Xβ+u\boldsymbol{Y}=\delta \boldsymbol{W} \boldsymbol{Y}+ \boldsymbol{X}\beta+\boldsymbol{u} Y=δWY+Xβ+u

u=λMu+ε,ε∼N(0,σ2I)\boldsymbol{u}=\lambda \boldsymbol{M u}+\boldsymbol{\varepsilon},\boldsymbol{\varepsilon} \sim \boldsymbol{N}(0,\sigma ^ 2I) u=λMu+ε,ε∼N(0,σ2I)

其中，WWW 与 MMM 分别为响应变量 YYY 与扰动项 uuu 的空间权重矩阵，二者可以相等。则此模型称为 带空间自回归误差项的空间自回归模型 (Spatial Auto Regressive Model with Spatial Autoregressive Disturbances，SARAR)。

上述 空间自回归模型（SAR）、空间误差模型（SEM）以及 SARAR 模型，被统称为 Cliff-Ord 模型 (Cliff and Ord, 1973, 1981; Ord, 1975)。

SAR 模型与 SEM 模型都是 SARAR 模型 的特例，分别对应于 λ=0\lambda=0λ=0 与 δ=0\delta=0δ=0的情形。

（2）模型的估计

对于 SARAR 模型，可进行 MLE 估计。 样本的对数似然函数为 :

ln⁡L(y∣λ,ρ,σ2,β)=−n2ln⁡2π−n2ln⁡σ2+ln⁡(abs⁡∣A∣)+ln⁡(abs⁡∣B∣)−12σ2(Ay−Xβ)′B′B(Ay−Xβ)\begin{aligned} \ln L\left(\boldsymbol{y} \mid \lambda, \rho, \sigma^{2}, \boldsymbol{\beta}\right)=&-\frac{n}{2} \ln 2 \pi-\frac{n}{2} \ln \sigma^{2}+\ln (\operatorname{abs}|\boldsymbol{A}|)+\ln (\operatorname{abs}|\boldsymbol{B}|) \\ &-\frac{1}{2 \sigma^{2}}(\boldsymbol{A} \boldsymbol{y}-\boldsymbol{X} \boldsymbol{\beta})^{\prime} \boldsymbol{B}^{\prime} \boldsymbol{B}(\boldsymbol{A} \boldsymbol{y}-\boldsymbol{X} \boldsymbol{\beta}) \end{aligned} lnL(y∣λ,ρ,σ2,β)=​−2n​ln2π−2n​lnσ2+ln(abs∣A∣)+ln(abs∣B∣)−2σ21​(Ay−Xβ)′B′B(Ay−Xβ)​

最大化仍可分两步进行，并进行迭代，与 SEM 模型 的情形类似。

但 MLE 估计量 存在三个缺陷：

首先，MLE 估计量 可能计算不方便，特别当空间权重矩阵 WWW 与 MMM 的维度较高时。

其次，MLE 估计量 的大样本理论尚不健全。

再次，如果 扰动项 ε 不服从独立同分布的正态分布，则 QMLE 估计量 可能不一致。

Arraiz 等人 (2010) 通过蒙特卡罗模拟证明，在异方差的情况下，QMLE 估计量 存在不一致。 Kelejian and Prucha (1999, 1998, 2004, 2010) 提出利用工具变量，通过广义矩估计来估计 SARAR 模型，称为 广义空间两阶段最小二乘法 (Generalized Spatial Two-stage Least Square，GS2SLS)。GS2SLS 的计算方便，不会因为样本容量太大而产生困难。

（1）模型的定义

（2）模型的估计

（1）模型的定义

（2）模型的估计

主要局限性：

思维导图：

处理区域间互动关系的另一方法是直接将地区间的距离作为变量引入回归模型，即将空间坐标或空间距离直接纳入模型，实现显式的空间依赖模型。例如：在国际贸易领域，常使用以下模型来考察两国间贸易量的决定因素：

Exportij=αYiβYjγDijδ\text{Export}_{i j}=\frac{\alpha Y_{i}^{\beta} Y_{j}^{\gamma}}{D_{i j}^{\delta}} Exportij​=Dijδ​αYiβ​Yjγ​​

其中，$$\text{Export}_{ij}$$ 代表 iii 国对 jjj 国的出口额，$$Y_i,Y_j$$ 分别为两国的 GDP，而 DijδD^{\delta}_{ij}Dijδ​ 是两国之间的地理距离。

此式很像牛顿万有引力公式，故称为 引力模型。对该方程两边取对数，在加上误差项，即可构成线性依赖模型：

ln⁡Exportij=ln⁡α+βln⁡Yi+γln⁡Yj−δln⁡Dij+ε\ln Export_{i j}=\ln \alpha+\beta \ln Y_{i}+\gamma \ln Y_{j}-\delta \ln D_{i j}+\varepsilon lnExportij​=lnα+βlnYi​+γlnYj​−δlnDij​+ε

一个基本的 GWR 模型可以指定为：

yi=βi0+∑k=1pβikxik+ϵi,i=1,…,n,y_{i}=\beta_{i 0}+\sum_{k=1}^{p} \beta_{i k} x_{i k}+\epsilon_{i}, \quad i=1, \ldots, n, yi​=βi0​+k=1∑p​βik​xik​+ϵi​,i=1,…,n,

其中 yiy_{i}yi​ 是位置 iii 处的响应变量，βi0\beta_{i0}βi0​ 是位置 iii 处的截距系数，xikx_{ik}xik​ 是位置 iii 处的第 kkk 个解释变量，

βik\beta_{ik}βik​ 是位置 iii 处第 kkk 个解释变量局部回归系数，ϵi\epsilon_{i}ϵi​ 是与位置 iii 相关的随机误差项。使用加权最小二乘估计，可以为每个位置的每个解释变量生成一组系数。在矩阵形式中，这由下式给出：

β^(i)=[X′W(i)X]−1X′W(i)y\hat{\boldsymbol{\beta}}(i)=\left[\boldsymbol{X}^{\prime} \boldsymbol{W}(i) \boldsymbol{X}\right]^{-1} \boldsymbol {X}^{\prime} \boldsymbol{W}(i) \boldsymbol{y} β^​(i)=[X′W(i)X]−1X′W(i)y

Li 2019 年的《快速地理加权回归（FastGWR）》： 地理加权回归 (GWR) 是一种广泛使用的工具，用于探索地理空间过程的空间异质性。 GWR 计算特定位置的参数估计值，这使得其校准过程需要大量计算。当前开源 GWR 软件可以处理的最大数据点数是标准桌面上的大约 15,00015,00015,000 个观测值。在大数据时代，这严重限制了 GWR 的使用。为了克服这一限制，我们提出了一种高度可扩展的开源 FastGWR 实现，该实现基于 Python 和消息传递接口 (MPI)，可扩展到数百万个观测值的数量级。 FastGWR 优化内存使用以及并行化以显着提高性能。为了说明 FastGWR 的性能，对来自洛杉矶市 Zillow 数据集的大约 130130130 万个单户住宅物业进行了特征房价模型校准，这是将 GWR 应用于这种规模的数据集的首次尝试。结果表明，随着高性能计算 (HPC) 环境中内核数量的增加，FastGWR 呈线性扩展。它还优于当前可用的开源 GWR 软件包，在标准桌面上速度大幅降低——最高可达数千倍。

wj(i)={[1−(dijb)2]2, if j∈{Ni}0, if j∉{Ni}w_{j}(i)= \begin{cases}{\left[1-\left(\frac{d_{ij}}{b}\right)^{2}\right]^{2},} & \text { if } j \in\left\{N_{i}\right\} \\ 0, & \text { if } j \notin\left\{N_{i}\right\}\end{cases} wj​(i)=⎩⎨⎧​[1−(bdij​​)2]2,0,​ if j∈{Ni​} if j∈/{Ni​}​

其中 dijd_{ij}dij​ 是 iii 和 jjj 之间的距离，bbb 是到第 NNN 个最近邻的距离，{Ni}\{N_{i}\}{Ni​} 是集合与 iii 的距离不超过 bbb 的所有观测值。位置 iii 的观测值的权重为 111 ，而那些比距离 bbb 更远的观测值的权重则为零。在拟合 GWR 模型之前，必须设定最近邻的数量，这进而又在上述平方核中定义了 bbb，也称为带宽。在整个研究过程中，可以使用 AIC 信息准则最小化来确定内核带宽，因为它提供了 拟合优度 和 自由度（或模型复杂度）之间的权衡（ Fotheringham、Brunsdon 和 Charlton 2002 ）。 GWR 特定的 AIC 计算为：

AICc=2nlog⁡e(σ^)+nlog⁡e(2π)+n{n+tr⁡(S)n−2−tr⁡(S)}AI C_{c}=2 n \log _{e}(\hat{\sigma})+n \log _{e}(2 \pi)+n\left\{\frac{n+\operatorname{tr} (\boldsymbol{S})}{n-2-\operatorname{tr}(\boldsymbol{S})}\right\} AICc​=2nloge​(σ^)+nloge​(2π)+n{n−2−tr(S)n+tr(S)​}

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

建模局部空间关系的另一种方法是使用贝叶斯层次模型来指定空间结构。Gelfand 等人（2003 年）提出的具有内生空间尺度的 贝叶斯空间变化系数模型（ Bayesian Spatially Varying Coefficient models，SVC ） 已很常见。为了对 SVC 的可能概括进行完整讨论，Finley、Banerjee 和 Carlin（2007 年）或 Banerjee、Carlin 和 Gelfand（2014 年）提供了非常出色的综述文章。 SVC 的核心是一个相关混合效应模型，其中 随机效应协方差 采用 空间相关核 来构建。

对于 SVC 模型，首先考虑响应变量 YYY 的模型，将 NNN 个位置观察到的 ppp 维解释变量编组为 N×pN\times pN×p 的设计矩阵 XXX 。 SVC 模型 将响应变量建模为全局效应 μ\muμ 和两个随机效应 ζ\zetaζ 和 ϵ\epsilonϵ 的函数。对于位置 iii，响应变量为：

yi=∑j=1pXij(μj+ζij)+ϵiy_{i}=\sum_{j=1}^{p} X_{i j}\left(\mu_{j}+\zeta_{i j}\right)+\boldsymbol{\epsilon}_{i} yi​=j=1∑p​Xij​(μj​+ζij​)+ϵi​

其中 ϵi\epsilon_{i}ϵi​ 是一个独立同分布正态误差项，均值为 000 ，方差为 τ2\tau^{2}τ2；ζij\zeta_{ij}ζij​ 是每个过程 jjj 的空间相关误差项。在这个模型中，μj\mu_jμj​ 代表过程 jjj 的平均全局效应。然后，通过计算每个单元的 μj+ζij\mu_{j}+\zeta_{i j}μj​+ζij​ 来恢复单元特定的参数估计值 βij\beta_{i j}βij​。局部扰动项中的空间结构是通过其分层先验协方差在 ζ\zetaζ 中引入的。为了以向量形式说明这个模型，Gelfand 等人 (2003) 将 XXX “拉伸”成一个N×NpN \times NpN×Np 的矩阵 XsX_{s}Xs​，其中每一行对应一个观测值，每组 ppp 列包含观测到的 ppp 个解释变量值。这种平铺将 ppp 长度的解释变量向量放在对角线上，在其他地方放置零。然后，ζ\zetaζ 可以成为 Np×1N p \times 1Np×1 的随机效应集合，其中第 iii 个块包含 ppp 个随机效应，用于每个 jjj 过程，μ\muμ 大小为 Np×1N p \times 1Np×1 并且包含过程的 ppp 长度向量，意味着重复了 NNN 次。通过这种平铺，模型可以被表述为：

Y=Xsμ+Xsζ+ϵY=X_{s} \mu+X_{s} \zeta+\epsilon Y=Xs​μ+Xs​ζ+ϵ

在此之后，ζ\zetaζ 的分层先验由一个 p×pp\times pp×p 的过程间协方差矩阵 TTT 和一个 N×NN\times NN×N 的单元间相关矩阵 H(ϕ)H(\phi)H(ϕ) （ 一个单独关于 ϕ\phiϕ 的函数）来构造。 Np×1N p \times 1Np×1 空间随机效应的完整协方差矩阵为：

ζ∼N(0,H(ϕ)⊗T)\zeta \sim \mathcal{N}(0, H(\phi) \otimes T) ζ∼N(0,H(ϕ)⊗T)

通常，H(ϕ)H(\phi)H(ϕ) 可能与 GWR 中的空间核函数类型相同。但这里的核直接影响随机效应的协方差，而不是对用于估计局部模型的数据进行加权。这是因为完整的 H(ϕ)⊗TH(\phi) \otimes TH(ϕ)⊗T 协方差矩阵定义了在所有 NNN 个单元处的过程之间的相关结构，体现了完整的 Np×NpN p \times N pNp×Np 单元-过程关系。与 GWR 一样，该核也可能更复杂，用于建模各向异性或替代距离度量。但是，核必须提供一个有效的协方差矩阵，以确保得到的协方差仍然有效。这排除了 GWR 中经常使用的 “自适应” 带宽，因为这会导致不对称的 HHH。对于其他参数，SVC 模型的常规先验被挑选用于条件共轭，并且通常通过 Gibbs 抽样 进行估计（Finley、Banerjee 和 Carlin 2007）。

为了将该模型推广到可与 MGWR 相媲美，需要一种对具有潜在不同尺度的多过程进行建模的方法。虽然 Finley（2011）提出了一种方法，但必须对其进行详细说明以使其可操作。这里重点考虑所有解释变量都存在局部变化的情况，但可以很容易推广到包括非变化项。首先，我们定义了一个新的 XXX 平铺，称为 XcX_{c}Xc​，它不同于 Gelfand 等人 (2003) 的单带宽 SVC 模型中使用的平铺。 堆叠发生在 ppp 的 NNN 个集合中。对于多过程模型，我们改为堆叠 NNN 的 ppp 个集合。

令 XcX_{c}Xc​ 是一个 N×NpN\times N pN×Np 矩阵，它水平堆叠对角化解释变量 XjX_{j}Xj​。这意味着 XcX_{c}Xc​ 是 pN×Np N \times NpN×N 对角矩阵：

Xc=[Diag⁡[X1]Diag⁡[X2]…Diag⁡[Xp]]X_{c}=\left[\begin{array}{llll} \operatorname{Diag}\left[X_{1}\right] & \operatorname{Diag}\left[X_{2}\right] \quad \ldots & \left.\operatorname{Diag}\left[X_{p}\right]\right] \end{array}\right. Xc​=[Diag[X1​]​Diag[X2​]…​Diag[Xp​]]​

然后，让 ζ\zetaζ 是一个 Np×1N p \times 1Np×1 的列向量，它通过堆叠每个ppp 过程的单元特定空间效应的 NNN-向量 来创建。由于每个过程都有自己的带宽，j∈{1,2,…,p}j \in\{1,2, \ldots, p\}j∈{1,2,…,p} 的每个 ζj\zeta_{j}ζj​ 都有自己的分层先验。这些被建模为可分离的以匹配 MGWR 规范，该方法将子过程视为有条件地相互独立。同时，XcζX_{c} \zetaXc​ζ 是 N×1N \times 1N×1 向量，它将正确的 “过程-单元” 随机效应 ζij\zeta_{ij}ζij​ 匹配到每个单元的正确解释变量 xijx_ {ij}xij​ ：

Xcζ=[∑jpx1jζ1j∑jpx2jζ2j⋮∑jpxnjζnj]X_{c} \zeta=\left[\begin{array}{c} \sum_{j}^{p} x_{1 j} \zeta_{1 j} \\ \sum_{j}^{p} x_{2 j} \zeta_{2 j} \\ \vdots \\ \sum_{j}^{p} x_{n j} \zeta_{n j} \end{array}\right] Xc​ζ=​∑jp​x1j​ζ1j​∑jp​x2j​ζ2j​⋮∑jp​xnj​ζnj​​​

有了这个，我们可以以向量形式为响应变量声明一个局部模型：

Y=Xμβ+Xcζ+ϵζ∼N(0,Diag⁡[H(ϕj)σj2]jp)ϵ∼N(0,τ2)\begin{gathered} Y=X \mu_{\beta}+X_{c} \zeta+\epsilon \\ \zeta \sim \mathcal{N}\left(0, \operatorname{Diag}\left[H\left(\phi_{j}\right) \sigma_{j}^{2}\right]_{j}^{p}\right) \\ \epsilon \sim \mathcal{N}\left(0, \tau^{2}\right) \end{gathered} Y=Xμβ​+Xc​ζ+ϵζ∼N(0,Diag[H(ϕj​)σj2​]jp​)ϵ∼N(0,τ2)​

其中 Diag⁡[M]ik\operatorname{Diag}[M]_{i}^{k}Diag[M]ik​ 是一个块对角矩阵，其中 MMM 沿对角线重复 kkk 次。 ζ\zetaζ 的协方差矩阵是块对角线 Np×NpN p \times N pNp×Np 矩阵，因此 ζ\zetaζ 可以通过分别收集每个过程的空间效应分布图来构造：

ζj∼N(0,H(ϕj)σj2)\zeta_{j} \sim \mathcal{N}\left(0, H\left(\phi_{j}\right) \sigma_{j}^{2}\right) ζj​∼N(0,H(ϕj​)σj2​)

使用这种块，可以通过将平均效应和该过程的随机扰动相加，再次恢复每个过程的局部效应：

βij=(μβ)j+ζij\beta_{i j}=\left(\mu_{\beta}\right)_{j}+\zeta_{i j} βij​=(μβ​)j​+ζij​

每个过程的空间随机效应分解成它们自己的分布这一事实，也确保了在边缘化或调节全局参数之后，每个过程的后验分布是独立的。我们可以定义一个在全局模型参数 {μ,τ,ζ}\{\mu, \tau, \zeta\}{μ,τ,ζ} 和过程特定参数 {ϕj,σj2}\{ \phi_{j}, \sigma_{j}^{2} \}{ϕj​,σj2​} 之间迭代抽样的 Gibbs 采样器。与单带宽 SVC 中使用的块不同，这允许多尺度 SVC 处理 N×NN\times NN×N 过程特定协方差矩阵中的 ppp，而不是处理完整的 Np×NpN p\times N pNp×Np 随机效应协方差单带宽模型。

特征向量空间滤波器基于这样的解释：修改后的连接矩阵的特征向量是一组可能正交和不相关的地图模式（即空间自相关程度）（Griffith 1996, 2011）。其第一个特征向量 E1E_1E1​ 为能够生成具有最大可实现 Moran's I 相关系数的地图模式的实数集；第二个特征向量 E2E_2E2​ 同样为能够生成具有最大可实现 Moran's I 相关系数的地图模式的实数集，不过与 E1E_1E1​ 之间不相关；这样继续直至与前面的 n−1n-1n−1 个特征向量不相关的 EnE_nEn​ 。改进后的连接矩阵 WWW 常被定义为：

(I−ll′/n)Cn(I−ll′/n)(5)(\mathbb{I}-ll^\prime/n)C_n(\mathbb{I}-ll^\prime/n) \tag{5} (I−ll′/n)Cn​(I−ll′/n)(5)

其中 III 是 n×nn \times nn×n 单元矩阵，lll 是 l′l^\primel′ 的 n×1n \times 1n×1 向量，CnC_nCn​ 是将研究空间互斥且完全划分的 nnn 个空间单位之间的二值连接矩阵。

通过选择从 WWW 派生的特征向量的子集并创建一个线性组合，可以生成一个能够解释残差中空间自相关的合成变量，并将其合并到线性回归模型：

y=β01+∑p=1PXpβp+∑k=1KEkβEk+ϵ\boldsymbol{y}=\beta_{0} \mathbf{1}+\sum_{p=1}^{P} \boldsymbol{X}_{p} \beta_{p}+\sum_{k=1}^{K} \boldsymbol{E}_{k} \beta_{E_{k}}+\boldsymbol{\epsilon} y=β0​1+p=1∑P​Xp​βp​+k=1∑K​Ek​βEk​​+ϵ

其中前两项给出了传统的线性回归规范（即截距以及系数和解释变量对的组合），第三项表示空间滤波器。这里，Ek\boldsymbol{E}_{k}Ek​ 是 KKK 个选定特征向量的总集合中的一个潜在特征向量。 Griffith 假设（6）中的规范还应该包括每个解释变量和每个选定特征向量之间的交互项，以便创建局部系数。他进一步假设添加新术语将提供一个近似等于 GWR 模型的规范，使得

y^gwr≈(β01+∑k0=1K0Ek0βk0)+∑p=1P(βp1+∑kp=1KpEkpβkp)⋅Xp+ϵ\hat{\boldsymbol{y}}_{g w r} \approx\left(\beta_{0} \mathbf{1}+\sum_{k_{0}=1}^{K_{0}} \boldsymbol{E}_{k_{0}} \beta_{k_{0}}\right)+\sum_{p=1}^{P}\left(\beta_{p} \mathbf{1}+\sum_{k_{p}=1}^{K_{p}} \boldsymbol{E}_{k_{p}} \beta_{k_{p}}\right) \cdot \boldsymbol{X}_{p}+\boldsymbol{\epsilon} y^​gwr​≈(β0​1+k0​=1∑K0​​Ek0​​βk0​​)+p=1∑P​​βp​1+kp​=1∑Kp​​Ekp​​βkp​​​⋅Xp​+ϵ

其中 y^\hat{y}y^​ 是 GWR 模型的预测值向量，⋅\cdot⋅ 表示逐元素矩阵乘法（即 Hadamard 矩阵乘法），kpk_{p}kp​ 表示来自集合的特征向量描述解释变量 ppp 的 KpK_{p}Kp​ 个特征向量。根据 Griffith (2008)，当 (6) 中的第一项和第三项在 (7) 中组合时，它们成为 GWR 截距系数，将在后续章节中更详细地探讨。最后，通过将 XpX_{p}Xp​ 分布在等式 (7) 中它们的预乘总和并重新排列各项，可以实现以下规范

y=β01+∑p=1PXp⋅1βp+∑k=1KEkβEk+∑p=1P∑k=1KXp⋅EkβpEk+ϵ\boldsymbol{y}=\beta_{0} \mathbf{1}+\sum_{p=1}^{P} \boldsymbol{X}_{p} \cdot 1 \boldsymbol{\beta}_{p}+\sum_{k=1}^{K} \boldsymbol{E}_{k} \beta_{E_{k}}+\sum_{p=1}^{P} \sum_{k=1}^{K} \boldsymbol{X}_{p} \cdot \boldsymbol{E}_{k} \beta_{p E_{k}}+\boldsymbol{\epsilon} y=β0​1+p=1∑P​Xp​⋅1βp​+k=1∑K​Ek​βEk​​+p=1∑P​k=1∑K​Xp​⋅Ek​βpEk​​+ϵ

其中第一项是回归截距，第二项表示解释变量及其全局系数，第三项和第四项分别表示全局特征向量和解释变量与特征向量之间的交互项。重要的是，全局解释变量是本规范要求构建交互项，最终允许构建局部系数。一旦选择了特征向量和交互项并估计了方程 (8)，每个解释变量的局部系数可以通过与每个全局解释变量相关的系数 ppp 和与与 GWR 框架不同，交互项基于相应的全局解释变量。然而，目前尚不清楚如何执行有关来自空间滤波的局部系数的假设检验，也不清楚如何应用多假设检验的校正。

基于空间滤波器的框架中的一个重要步骤是选择特征向量和交互项的子集。该方法需要基于统计显着性在所有交互项和全局特征向量中的前向逐步回归变量选择算法。首先，仅选择显示中度至高度空间自相关（即，莫兰系数 (MC) 大约大于 0.250.250.25）的特征向量。然后，对于逐步例程的每次迭代，在模型中测试每个候选变量。选择产生最小 P 值的变量，如果任何其他变量变得不显着（P 值大于 0.10.10.1），则将它们从模型中删除。该算法继续进行，直到没有变量可以添加到模型中，该模型也产生小于 0.10010.10010.1001 的关联 P 值，此时例程终止并且可以处理局部系数 (Griffith 2008)。在后面的部分中，本规范与作为本研究的一部分开发的这种逐步选择例程的两种替代均方误差最小化 (MSE) 选择变体进行比较。