空间变系数模型的新旧方法

【摘要】 本文比较了空间异质性建模中空间变化参数模型的一些主要方法，其中包括：（1）21 世纪前提出的传统方法，包括：空间展开模型、空间自适应滤波模型和地理加权回归模型。（2）21 世纪初以来出现的一些新方法，包括：空间平滑过渡自回归模型、空间高斯过程模型、含自回归过程的随机参数模型。（3）一些通用变参数模型方法在空间异质性建模中的应用，包括：空间样条方法等。文中采用人工合成数据，以图形方式展示了不同方法之间的差异。注意：不知为何原因，本文未提及空间滤波方法。

【原文】 D. M. Lambert, “Old and new approaches for spatially varying coefficient models,” Review of Regional Studies, vol. 51, no. 2, pp. 113–128, 2021, doi: 10.52324/001c.27969.

---

《空间回归模型综述》 一文对面元数据的空间回归模型进行了综述，总体分为全局空间回归模型和局部空间回归模型，文中也谈到全局空间回归模型重点讨论邻近的目标变量、因变量、以及误差项之间的空间交互效应,探究的核心是空间依赖性。在此类模型中，空间异质性通常隐含存在于协方差函数的结构建模中。局部空间回归模型与全局模型不同，它显式地对模型参数随空间位置的变化建模，也就是说，它研究模型本身（ 模型结构 、 模型参数 等）的空间变化规律，因此探究的核心是一种 空间异质性 或 空间结构不稳定性 的显式建模方法。

1 简介

由于局部空间回归模型中的参数会随着空间位置而变化，因此其推断和预测面临较大的计算复杂度挑战，因此在过去的有段时间内发展比较缓慢。进入 21 世纪后，一些计算和软件等方面的新进展，重新激发了业内人士对 空间局部回归方法 (Local Regression Models, LRM) 的兴趣。该方法允许空间单元根据自身区位对控制做出不同的反应，这一点对区域科学家来说非常具有有吸引力。相对于全局回归模型而言，空间局部回归模型放宽了 “所有地点对局部条件、资源或其他变化因素具有相似响应” 的基本假设。毕竟，在看地图时很容易得出感性的结论：某些地点的迁移率高于地区平均水平、某地的企业进出速度具有局部特质、政府的政策干预在不同区域会产生不同影响。

研究人员已经提出了多种空间局部回归模型，包括：

• Casetti (1972) 的 空间展开模型 （ Spatial Expansion ）
• Foster 和 Gorr (1982) 的 空间自适应滤波器 ( Spatial Adaptive Filter, SAF)
• Cleveland 和 Devlin (1988) 的 地理加权回归 (Geographically Weighted Regression, GWR )
• Pede 等 (2014) 的 平滑过渡模型 ( Smooth Transition AutoRegression, STM )
• Finley 和 Banerjee (2020) 的 贝叶斯空间高斯过程模型 ( Spatial Gaussian Process Model, SGP)
• Congdon (2020) 的 空间自回归型随机参数模型 (Spatial Autogressive Random Parameter Model, SRP)

本文主要侧重于空间数据，未考虑时空结合的局部回归模型，不过本文结果大部分可以推广到时空局部回归模型。 局部回归模型具有直观的吸引力，并且大多数方法都相对容易估计，但研究人员在选择适合手头问题的模型时应谨慎，最好使用探索性策略对样本内性能作出模型评估和比较后再作出选择。

2 背景

（1）传统空间异质性建模方法

• Whittle (1954) 将空间自回归过程引入线性回归是对横截面单元之间的空间相互作用进行建模的早期努力。 Whittle 将局部结果的空间滞后作为内生变量包含在线性响应模型中。
• 大约 20 年后，Cliff 和 Ord (1973) 将 Whittle 的自回归模型重新定义为具有相关误差的模型。
• Whittle、 Cliff 和 Ord 的贡献现在分别被称为 空间自回归模型 (SAR) 和 空间误差模型 (SEM)。 Anselin (2003) 曾对相关的模型及变体进行了总结。
• SAR 和 SEM 方法通过自协方差结构对空间交互进行建模，但对局部变化的边缘响应是全局性的；即单个估计的斜率系数适用于所有空间单元。
• LeSage 和 Pace (2009) 将全局效应分解为直接和间接效应，并给出了使用 “Leontief 逆矩阵” 实现分解的框架程序。分解后的SAR 模型边缘响应，提供了关于反馈和溢出效应的详细信息，同时 SEM 的边际效应是全局性的。

（2）局部回归模型（空间变系数模型）

局部回归模型明确地对空间异质性进行模，而无需空间自相关或交互作用（Cho et al. (2010)）。也就是说，局部回归模型解决的空间异质性类型与 SAR 和 SEM 等模型的异质性类型有本质区别。

• 对于 SAR 和 SEM 模型，异质性通过协方差结构间接建模，或者对于 SAR，通过 LeSage 和 Pace 的边际效应分解进行建模。
• 局部回归模型中，位置直接与局部变化的生成或需求函数有关，而这可能是截距和斜率等系数随空间变化的根本原因 (McMillen (1996))，这意味着每一个空间单元的斜率和截距项都是可估计的。

一种实现空间异质斜率和截距系数的直接方法，是人为地对空间单元进行分类，在每一种类型中，将响应模型视为一种特定于该类型的固定效应（McGranahan 等（2011 年））。固定效应方法需要利用外部信息对空间单元进行分类，例如 “北”与“南”、“大都市”与“非大都市”等。此方法的优点是简单，而且对于横截面数据而言，固定效应的 “空间机制” 模型易于估计，非常适合做跨机制的响应比较。此方法的问题在于：空间单元的分类方案无法代表数据生成过程，或者不是数据生成过程的组成部分；此外，该方法还保持这样的假设，即聚集到同一类别组中的空间单元共享相同的响应机制。

本文考虑的局部回归模型以线性横截面模型为主，可以概括为，

$$\begin{array}{r} y_i=\sum_{k=1}^K \beta_{i k}\left(\mathbf{s}_i\right) \cdot x_{i k}+\varepsilon_i, i \in 1 \ldots N, \\ \varepsilon_i \sim \operatorname{iid}\left(0, \sigma^2\right) \end{array}$$

$y_i$ 是在位置 $i$ 观测到的结果，$N$ 是横截面中的位置数，$\mathbf{x}_i$ 是在位置 $i$ 观测到的控制向量，$\beta_i(\mathbf{s}_i)$ 代表 $N$ 个长度为 $k$ 的边际效应向量，$\mathbf{s}_i$ 是定位位置的坐标向量，$\sigma^2$ 是误差的方差。本文所考虑的局部回归模型，其显着特征在于系数 $\beta_i(\mathbf{s}_i)$ 的参数化方式。

3 局部回归模型

3.1 空间展开模型

Casetti (1972) 的空间展开模型将 $\beta_i(\mathbf{s}_i)$ 指定为 “位置坐标 × 参数” 交互形式的多项式函数。例如，如果使用二阶多项式展开来参数化 $\beta_i(\mathbf{s}_i)$，则：

$$\begin{array}{r} \beta_{i k}\left(\mathbf{s}_i\right)=\alpha_{0 k}+\alpha_{1 k} \cdot x c_i+\alpha_{2 k} \cdot y c_i+ \\ \alpha_{3 k} \cdot x c_i^2+\alpha_{4 k} \cdot y c_i^2+\alpha_{5 k} \cdot x c_i \cdot y c_i \end{array}$$

其中 $\alpha$ 是全局参数，$\mathbf{s}_i \in (yc_{i}, xc_{i})$ 是位置坐标。高阶多项式是可以的，但会引入共线性。可以省略交互作用和高阶项以减少共线性引起的问题。该方法受到的主要批评是：其与异质性来源之间的关系没有明显的理论依据。该特点（连同所引入的共线性）使得空间展开模型的推断极具挑战性。 Anselin (2003) 为 Casetti 的空间展开模型提供了另外的视角。

2.2.空间自适应滤波 (SAF)

Foster 和 Gorr (1982) 使用空间滤波过程来估计特定位置的系数。该方法借鉴了时间序列文献中使用的、估计“滚动系数”的方法，即系数是随时间演变的移动窗口均值。空间近邻类似于时间序列中使用的、将序列划分为数据子集的局部窗口。

局部系数估计为：

$$\beta_{i k}\left(\mathbf{s}_i\right)=\frac{1}{n_i} \sum_{j \in n_i} \beta_{i j}\left(\mathbf{s}_i\right)$$

其中 $n_i$ 是位置 $i$ 的窗口中包含的邻居数。换句话说，$\beta_{ik}(\mathbf{s}_i)$ 是基于其 $\beta_{ij}(\mathbf{s}_j)$ 邻居的加权平均值的预测。邻居是使用具有预定阈值的连续距离测量或通过邻接来定义的。迭代过程用于最小化响应变量的误差方差。

2.3 地理加权回归

Cleveland 和 Devlin (1988) 引入了地理加权回归 (GWR)。 McMillen (1996) 和 Fotheringham 等 (2003) 将地理加权回归扩展到了广义线性模型，包括 Poisson、logit 和 probit 的地理加权回归估计。

位置 $i$ 处的响应为：

$$y_i=\sum_k \beta_{i k}\left(w\left(\mathbf{s}_i\right)\right) \cdot x_{i k}+\varepsilon_i \tag{4}$$

其中 $\mathbf{s}_i$ 是位置 $i$ 与其 $N - 1$ 个邻居之间的距离向量， $w(\mathbf{s}_i)$ 是特定于每个位置的核权重函数，以及先前定义的其他符号。核函数确定了位置 $i$ 的邻居对其参数估计的影响范围和强度。

典型的核包括：

$$\begin{array}{r} \text { Gaussian : } w\left(\mathbf{s}_i\right)=\phi\left(\rho \cdot\left[\left\|\mathbf{d}_i\right\|^2\right] / \widehat{\sigma}_{d i}\right) \\ \text { Exponential : } w\left(\mathbf{s}_i\right)=\exp \left(-\rho \cdot\left\|\mathbf{d}_i\right\|^2\right) \\ \text { Bi-square }: w\left(\mathbf{s}_i\right)=\left(\mathbf{d}_i<d_{\max }\right) \cdot\left[1-\left(\mathbf{d}_i / d_{\max }\right)^2\right]^2 \end{array} \tag{5}$$

其中 $\phi$ 是标准正态概率密度函数，̂$\widehat{\sigma}_{d i}$ 是距离向量的标准偏差，$\mathbf{d}_i$ 是位置 $i$ 和其邻居之间的距离向量，$\rho$ 是一个空间参数。

双平方核是局部核，因为它通过生成观测子集来工作，而不是像 指数核 和 高斯核 那样使用所有 $N -1$ 个邻居。双平方参数 $\mathbf{d}-{max}$ 决定了用于估计局部向量的邻居数量。核参数使用交叉验证程序确定（Fotheringham 等（2003 年）；Bivand 等人（2008 年））。一旦核参数被恢复，地理加权回归通过为每个位置估计局部参数向量展开后续任务。

地理加权回归系数的矩阵是:

$$\boldsymbol{\beta}_i=\left(\mathbf{X}^{\prime} \mathbf{W}(i)^{\prime} \mathbf{X}\right)^{-1} \mathbf{X}^{\prime} \mathbf{W}(i) \mathbf{y} \tag{6}$$

与其他局部回归模型一样，地理加权回归具有直观的吸引力，现在有多种软件环境可供研究人员使用，但这并非没有问题。

• 问题一： 推断问题。随着样本量增加，模型参数的数量以相同的速度增加，造成推断困难。地理加权回归的这一特点与标准的渐近理论不一致，后者认为随着样本量的增加，参数会收敛到其真实值。此外，相同的观测值会被重复用于估计特定位置的参数，这违反了关于模型误差项的独立同分布假设。尽管计算特定位置的标准差和 $t$ 检验在机制上很简单，但很难确定假设检验的合理自由度。

• 问题二： 共线性问题（Farber 和 Ṕaez (2007)；Wheeler 和 Tiefelsdorf (2005)），在极端情况下，人们理解的映射系数中的结构化地理模式，实际上有可能是地理加权回归估计的共线性问题造成的假象。为了解决此问题，新的研究将套索和岭回归技术扩展到 GWR（Wulandari 等人（2017 年）；He 等人（2021 年）），在此不作赘述。

• 问题三： 核选择问题。与从业者常选择的空间加权矩阵 $W$ 不同（LeSage and Pace (2009)），地理加权回归估计方法对所使用的核（用于重新加权相邻位置对位置 $i$ 估计的影响）非常敏感。图 1 使用模拟数据展示了不同核选择对局部估计的影响。该模型包括一个常数项和一个协变量，该协变量服从均值为 $0$、方差为 $1$ 的标准正态分布；误差项也是一个均值为 $0$、方差为 $1$ 的标准正态随机变量。图中的网格大小为 $40 × 40$，共产生了 $1,600$ 个位置。局部自适应双平方核与高斯核及指数核（两者均为连续核）之间的区别非常明显。连续核使用所有 $N - 1$ 个观测值估计位置 $i$ 处的参数，而在双平方核估计中则包含更少的观测值。从双平方核生成的数据集的截断性质导致一种子弹模式，而由连续核函数生成的模式更平滑。这些机制差异表明，从业者应该仔细考虑数据生成过程（Data Geraration Process，DGP），因为它与局部及全局的随机空间场有关，但也会涉及工作经验问题。在某些情况下，数据生成过程被被表示为主要源自局部事件（而不是面向整个区域的平滑、连续的全局过程）可能更好。

由于上述原因，一些人认为使用地理加权回归和其他局部回归方法最适合于探索性空间分析（例如，Wheeler 和 Calder（2007）；Farber 和 Ṕaez（2007））。事实上，Roger Bivand 在 R 中的地理加权回归例程包括这个警告（Bivand (2020)）。

研究人员已经开发了用于检验地理加权回归参数平稳性和拟合优度的自由度调整方法（ Leung 等 (2000)），但该方法似乎并未扩展到特定位置的假设检验。

图 1：地理加权回归核的比较

2.4 平滑过渡回归

Pede 等 (2014) 将时间序列文献中开发的 状态切换 或 平滑过渡自回归模型 (STM) 扩展到了空间数据。空间空间平滑过渡模型 包含 SAR(1)、SEM(1) 和 ARAR(1) 空间过程模型。不失一般性，这里考虑的空间平滑过渡模型排除了空间自回归过程，只关注局部回归模型参数的估计。

响应方程是非线性的：

$$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + g(\mathbf{s}; \gamma, c) \mathbf{X} \boldsymbol{\delta} + \boldsymbol{\varepsilon} \tag{7}$$

其中 $g(\mathbf{s}; \gamma, c)$ 是连续累积密度函数 (cdf)，它将观察结果分类为两种空间状态，$\mathbf{s}$ 是距离向量或其他特定于位置的数据，$\gamma$ 是 cdf 的斜率，$c$ 是指示 cdf 中地理位置中值的参数。模型参数随位置变化为 $\beta_{ik} = \beta_{k} + g(\mathbf{s}_i; \gamma, c) · \delta_k$，其中 $\delta_k$ 是斜率（或截距）偏移。

Pede 等使用逻辑斯谛 cdf 来指定转换函数：

$$g(\mathbf{s}; \gamma, c) = \frac{1}{1 + exp(−\gamma · [s − c])} tag{8}$$

当 $\gamma = \delta = 0$ 时，模型简化为空间线性回归。随着 $\gamma$ 的增加，$g(\mathbf{s}; \gamma, c) \rightarrow 1$，线性模型表现得好像是一个虚拟变量回归，其中 $\delta_k$ 是状态之间的差异（Pede 等人（2014 年）；Brown 和 Lambert（2016 年））。对于 $\gamma$ 的中间值和 $\delta \neq 0$，公式 7 是一个局部回归模型。Pede 等将转换函数的斜率和位置参数限制为在所有协变量中相同。一个直接的扩展是允许 $g(·)$ 的斜率和位置参数在 $\mathbf{X}$ 中的协变量之间变化（即 $g(\mathbf{s}; \gamma_k, c_k)$ ）。

2.5 空间高斯过程模型

贝叶斯计量经济学的进步，以及随之而来的软件可用性和计算资源的进步，使得对空间变化参数（被视为相对于总体平均响应 $\boldsymbol{\beta}$ 的随机变化）建模成为可能。某个位置处参数随空间的变化（相对于总体均值）被构造为空间高斯过程 (SGP)，通常带有 半变异函数（Semi-variogram）。贝叶斯克里金法 和 空间自回归型空间高斯过程 都可以在 R 或 R-Stan 环境中进行估计。

2.5.1 贝叶斯克里金法和空间高斯过程

贝叶斯克里金法使用 重采样方法 生成 空间变化估计 与 空间协方差参数 的联合后验分布。贝叶斯克里金法的早期例子包括 Handcock 和 Stein (1993)、Gelfand 等（2003）、Wheeler 和 Calder（2007）。最新的应用包括 Park 等(2019)、Finley 和 Banerjee (2020) 、McElreath (2020) 。

贝叶斯克里金法（和更普遍的贝叶斯方法）对自由度问题具有鲁棒性，而这恰恰是用于估计空间变参数模型的频率派方法存在的问题。尽管如此，贝叶斯克里金法是计算密集型的，估计空间变化参数所需的时间随着位置单元数呈指数增长。出于这个原因，明智地选择先验会有助于收敛。 Finley 和 Banerjee (2020) 提供了一个 R 包 spSVC，用于估计贝叶斯克里金空间高斯过程。

使用贝叶斯克里金法的空间高斯过程的广义模型结构是分层的（Finley 和 Banerjee (2020)）：

$$\begin{array}{r} y_i \sim \mathcal{N}\left(\mu_i, \sigma\right) \\ \mu_i=\mathbf{x}_i \boldsymbol{\beta}_i \\ p\left(\beta_{i k}\right) \sim \mathcal{M V N}\left(\bar{\beta}_k, \mathbf{V}_k\right) \end{array}$$

其中 $y_i$ 是正态分布的响应变量，均值函数为 $\mu_i$，标量方差为 $\sigma$。平均响应是线性的，并且是 $x_i$ 中特定于位置的协变量与相应参数向量 $\boldsymbol{\beta}_i$ 的函数。

总体均值（ $\bar{\beta}$ ）的先验是正态分布，其均值为零，标准差为 $10$（一个方差相对分散的先验）。 $\sigma$ 的一个合理先验是衰减率为 $1$ 的指数分布（McElreath (2020)）。换言之，响应与总体均值之间的平均偏差预计在 $±1$ 标准偏差范围内。

主效应参数的 $k$ 长度向量的先验是多元正态，其总体均值为 $\bar{\beta}_k$)，协方差结构为 $\mathbf{V}_k$。空间协方差（控制一个位置与其 $N - 1$ 个邻居之间效应偏差的相关性）使用半变异函数进行建模，其中 $g (d_{ij}, \rho_k)$ 是指数函数，可以使用高斯、球面或其他可接受的半变异函数来指定 $g(·)$。参数 $\sigma^2_{\beta_k}$ 和 $\eta^2_{\beta_k}$（分别）是半变异函数的块金效应和基台效应，而 $\rho_k$ 是其范围（Cressie (1993)）。空间协方差先验也是指数的，衰减率为 $1$。

2.5.2  空间自回归型高斯过程模型

贝叶斯克里金法的一个有趣扩展是 Congdon (2020) 引入了对空间自回归型结构的空间高斯过程估计。具有空间变参数的 SAR 型结构对于区域科学家可能更熟悉。Congdon 的方法也是一种分层贝叶斯模型，除了总体水平主效应及其标准差的先验之外，平均响应函数类似于贝叶斯克里金空间高斯过程。 Congdon 提供了 R-Stan 代码来估计 SAR 型空间高斯过程的版本。

SAR 型空间高斯过程模型为：

$$\begin{array}{r} y_i \sim \mathcal{N}\left(\mu_i, \sigma\right) \\ \mu_i=\mathbf{x}_i \boldsymbol{\beta}_i \\ p\left(\beta_{i k}\right) \sim \mathcal{M V \mathcal { N }}\left(\beta_k, \sigma_k^2 \cdot\left[\mathbf{B}_k^{\prime} \mathbf{B}_k\right]^{-1}\right) \\ p\left(\bar{\beta}_k\right) \sim \mathcal{N}(0,10) \\ p(\sigma), p\left(\sigma_k\right) \sim \operatorname{Exponential}(1) \\ \rho_k=2 \cdot\left(z_k-0.5\right) \\ p\left(z_k\right) \sim \operatorname{Beta}(2,2) \end{array}$$

其中 $\mathbf{B}_k = I - \rho_k \mathbf{W}$，$\mathbf{W}$ 是一个行归一化的空间权重矩阵，$\mathbf{I}$ 是一致的单位矩阵，$\rho_k \in (-1, 1)$ 是 $\beta_{ik}$ 的空间自回归参数。在这里，空间自回归系数被参数化为一个 $\text{Beta}$ 随机变量的函数，但其他先验也是可以接受的，例如 逻辑斯谛 或 均匀分布。

5 结论

鉴于空间数据的固有性质，使得区域科学家处于一个独特的位置，他们可以为空间变参数模型的发展做出贡献。本演示的目的是对不同类型的局部回归模型进行定性比较，同时尽量减少可归因于规模、复杂性和其他模型制品的差异。每种局部回归模型都有自己的优势和劣势，本文考虑的局部回归模型之间，差异主要体现在于定义和构建空间过程的方式。

从计算效率的角度来看，Casetti 的空间展开和地理加权回归优于其他局部回归模型。空间展开和空间自适应滤波的特殊定义，使其很难对参数异质性的模式作出合理解释。空间展开过程将参数变化建模为多项式趋势。如果实际的数据生成过程与展开函数无关，那么经验模型可能会被严重错误地指定。同样的警告也适用于空间自适应滤波局部回归模型。

如果理论表明某些控制变量会导致跨空间的结构中断，那么空间平滑过渡模型与此概念是一致的。在此处提供的合成数据示例中，包含在状态函数中的空间变量旨在强调状态转换潜力。在实践中，可以包括状态分裂函数中协变量的空间滞后（Pede 等（2014 年）），或某个位置到重要要素之间的距离。或者，协变量的水平可以包含在函数 $g(·)$ 中，例如人口密度。

计算最复杂的局部回归模型是空间高斯过程和空间自回归类型的估计器。与地理加权回归一样，这两个估计器需要指定协方差过程。是否应该在空间高斯过程的协方差项中包含块金效应的决定并不简单，正如其对参数异质性的平滑效应的强度所证明的那样。忽略块金效应可能会导致参数过度平滑。相反，空间自回归型局部回归模型通过参数的协方差项间接包含块金效应。

Casetti 的空间展开模型以及 Gorr 和 Foster 的空间自适应滤波最先提出了一种使用空间变参数对异质性进行建模的方法。空间展开模型可以使用普通最小二乘法估计。但多项式展开项引入了共线性。此外，展开程序与研究人员面对的数据生成过程可能相关，也可能不相关。空间自适应滤波更难估计，需要一个优化程序来分别调整每个位置的参数向量，与空间展开方法一样，空间自适应滤波可能与真实数据生成过程脱节。

地理加权回归是一种流行的替代方法，用于对具有空间变参数的异质性进行建模。与空间展开和空间自适应滤波方法不同，地理加权回归不依赖特定定义来估计空间变化的参数。当前算法在使用大型空间数据集估计地理加权回归模型方面是有效的。但研究人员不应在公共领域软件中使用默认选项，而应仔细考虑用于异质性建模的核类型，以及核选择对数据生成过程假设的影响。

空间平滑过渡模型是将异质性建模为空间变参数问题的另一种替代方法。空间平滑过渡模型的一个优点是能够提供一种方法来测试跨空间参数变化的结构中断。研究区域增长动态中的集聚经济作用的研究人员可能会对此特点感兴趣。未来的研究可以考虑将状态函数修改为与回归模型中包含的每个协变量相关的函数。

软件可用性的进步、对开源存储库的访问增加以及计算能力已经普及了用于解决复杂协方差结构和模型设计的分层贝叶斯建模方法。除了空间变参数之外，这些进步为建模复杂的空间协方差结构打开了大门。但与地理加权回归一样，必须注意其核以及半变异函数（在贝叶斯克里金模型的情况下）的选择。此外，尽管计算的进步使空间贝叶斯分层模型的估计更容易处理，但更大的空间数据集仍然在收敛时间和参数稳定性方面带来计算挑战。