基于空间滤波的大型数据集空间变系数建模

〖摘 要〗 虽然空间变系数 (SVC) 建模在应用科学中很流行，但其计算负担很大。如果考虑空间变系数的多尺度属性，则尤其如此。鉴于此背景，本研究开发了一种基于 Moran 特征向量的空间变系数 (M-SVC) 建模方法，可有效地估计多尺度空间变系数模型。该估计通过 (1) 秩降低、(2) 预压缩和 (3) 顺序似然最大化来加速。步骤 (1) 和 (2) 从似然函数中消除样本大小 N；在这些步骤之后，似然最大化成本与 N 无关。步骤 (3) 进一步加速似然最大化，因此即使空间变系数的数量 K 很大，也可以估计多尺度空间变系数模型。通过蒙特卡罗模拟实验将 M-SVC 方法与地理加权回归 (GWR) 进行比较。这些模拟结果表明，当 N 很大时，本文方法比地理加权回归快得多，尽管数值估计了 2K 个参数，而地理加权回归仅数值估计了 1 个参数。然后，将所提出的方法应用于土地价格分析作为说明。开发的空间变系数估计方法在 R 包 “spmoran” 中实现

〖原 文〗 Murakami, D. and Griffith, D.A. (2019) ‘Spatially varying coefficient modeling for large datasets: Eliminating N from spatial regressions’, Spatial Statistics, 30, pp. 39–64. Available at: https://doi.org/10.1016/j.spasta.2019.02.003.

1. 导言

在地统计学（例如，Cressie，1993 ^[8]）和空间计量经济学（例如，LeSage 和 Pace，2009^[36]）中提出的空间统计方法，通常需要计算复杂度为 $\mathcal{O}(N^3)$ 的密集协方差矩阵的逆，其中 $N$ 是样本量。这些方法不适用于大样本。

为了克服这一局限性，大空间数据建模得到了迅速发展。

地统计学家主要针对空间插值问题研究大数据建模。稀疏方法和低秩方法是其中的典型代表（参见 Heaton 等 2017 年的评论^[32]）：

• 稀疏方法包括
  • 协方差锥化（Furrer 等，2006 年^[20]）
  • 使用随机偏微分方程（Lindgren 等，2011 年 ^[38]）
  • 最近邻高斯过程（Datta 等，2016 年^[10]）
• 低秩方法包括
  • 固定秩克里金法（Cressie 和 Johannesson，2008 年^[9]）
  • 预测过程建模（Banerjee 等，2008 年^[3]）
  • 多分辨率近似（Nychka 等，2015 年^[49]）。

空间计量学家主要针对存在空间依赖性的回归问题研究大数据建模。例如：

• 基于多项式的近似（Smirnov 和 Anselin，2001 ^[53]；Griffith，2004 ^[26]）
• 组合似然近似（Arbia，2014 ^[2]）
• 其他似然近似（例如，LeSage 和 Pace，2007 ^[35]）。
• Burden 等 (2015) ^[5] 提出了一种低秩方法。
• 空间两阶段最小二乘法（Kelejian 和 Prucha，1998 ^[34]）被认为是一种计算效率高的方法。

虽然空间依赖性和空间异质性是空间数据的两个普遍性质（Tobler, 1970 ^[56]; Goodchild, 2004 ^[23]; Anselin, 2010 ^[1]），但之前的这些研究都更关注前者，却忽视了后者。尽管如此，空间异质性对进行恰当地空间数据建模至关重要（参见 Fotheringham 等，2002 ^[18]；Geniaux 和 Martinetti，2018 ^[22]；Debarsy 和 Yang，2018 ^[11]）。

鉴于此背景，本研究重点关注空间变系数 (SVC) 建模，这是空间异质性建模的代表性方法。大量研究面向空间变系数建模（ 见 第 2 节）。最近，Fotheringham 等 (2017)^[19] 和 Murakami 等 (2018)^[43] 揭示了考虑空间变系数的多尺度特性，在准确估计它们时的重要性；例如，第一个空间变系数可能具有比其他空间变系数更大比例的地图模式。然而， 如果考虑多尺度性质，计算成本会进一步迅速增加，不仅取决于 $N$，还取决于空间变系数的数量，如何在计算上高效估计多尺度空间变系数尚不明确。

本研究的目的正是针对此局限性。具体来说，我们开发了一种方法，通过扩展 Moran 的特征向量空间过滤方法（Griffith 2003 ^[25]；Murakami 等，2017 年 ^[45]），能够以计算有效的方式估计多尺度空间变系数模型。

随后的部分组织如下。

• 第 2 节 简要回顾现有空间变系数建模方法，并比较了其计算特性。
• 第 3 节 介绍了我们的空间变系数模型。
• 第 4 节 开发了一种快速空间变系数估计方法。
• 第 5 节 将所提出方法的计算时间与现有方法进行比较。
• 第 6 节 介绍了一个实证应用。
• 第 7 节 总结了我们的讨论。

2. SVC 建模方法

假设响应变量 $y(\mathbf{s}_i)$ 在分布在研究区域 $D ⊂ \mathbb{R}^2$ 中的 $N$ 个样本点 $\mathbf{s}_i \mid i \in \{1,\ldots ,N\}$ 处被采样。线性空间变系数模型的公式如下：

$$\begin{align*} &y(\mathbf{s}_i) = \sum^K_{k=1} x_k(\mathbf{s}_i) \beta_k(\mathbf{s}_i) + \varepsilon(\mathbf{s}_i) \tag{1}\\ &\mathbb{E}[\varepsilon(\mathbf{s}_i)] = 0\\ &\mathbb{V}ar[\varepsilon(\mathbf{s}_i)] = \sigma^2 \end{align*}$$

其中，$x_k(\mathbf{s}_i)$ 表示第 $k$ 个解释变量，$\beta_k(\mathbf{s}_i)$ 表示第 $k$ 个变系数，$\sigma^2$ 为方差参数。空间变系数建模旨在估计和推断 $\beta_k(\mathbf{s}_i)$。

（1）地理加权回归

地理加权回归 (GWR; Brunsdon 等, 1996 ^[6]) 是空间变系数建模的代表性方法。地理加权回归通过为附近样本分配比远处测点 $\mathbf{s}_i$ 的样本更多的权重来估计 $\beta_k(\mathbf{s}_i)$。估计 $\boldsymbol{\beta}(\mathbf{s}_i) = [\beta_1(\mathbf{s}_i), \ldots, \beta_k(\mathbf{s}_i)]^T$ 由下式给出（其中 T 表示转置算子）：

$$\hat{\boldsymbol{\beta}}(\mathbf{s}_i) = [\mathbf{X}^T \mathbf{G}(\mathbf{s}_i)\mathbf{X}]^{-1} \mathbf{X}^T \mathbf{G}(\mathbf{s}_i) \mathbf{y}$$

其中 $\mathbf{X}$ 是解释变量的 $N × K$ 矩阵，$\mathbf{y}$ 是响应变量的 $N × 1$ 向量，$\mathbf{G}(\mathbf{s}_i)$ 是 $N × N$ 对角矩阵，其第 $j$ 个元素 $g(\mathbf{s}_i, \mathbf{s}_j)$ 表示分配到第 $j$ 个样本的权重。$g(\mathbf{s}_i, \mathbf{s}_j)$ 由距离衰减核给出。例如，Wheeler 和 Calder (2007) ^[58] 以及 Wheeler 和 Waller (2009) ^[60] 使用指数核 $g(\mathbf{s}_i, \mathbf{s}_j) = \exp(−d(\mathbf{s}_i, \mathbf{s}_j)/b)$，其中 $b$ 是确定尺度的带宽变系数；当 $b$ 大时，$\hat{\boldsymbol{\beta}}(\mathbf{s}_i)$ 为大尺度地图模式，$b$ 小时为小尺度地图模式。通过最小化修正后的 Akaike 信息准则 (AIC) 或最小化交叉验证 (CV) 的预测误差，我们可以估计带宽 $b$。地理加权回归已广泛应用于环境研究（例如，Harris 等，2010a ^[30]；Dong 等，2018 ^[12]）、犯罪学（例如，Cahill 和 Mulligan，2007 ^[7]）、流行病学（例如，Nakaya 等，2005 ^[48]），以及许多其他领域。

单一带宽隐含地强加了一个强有力的假设，即所有空间变系数具有相同的空间变化尺度。但是，一些系数可能比其他系数具有更大的尺度模式。为了考虑这种多尺度属性，Yang (2014) ^[61]和 Fotheringham 等（2017）^[19] 提出了灵活带宽地理加权回归（FB-GWR，也称为多尺度地理加权回归），它用空间变系数特定带宽 $b_k$ 代替了常见的 $b$。多尺度地理加权回归通过反向拟合算法估算，该算法顺序迭代每个 $b_k$ 的 AIC 最大化（或 CV）直到收敛。最近，Lu 等 (2018)^[41] 提出了一种计算效率高的多尺度地理加权回归模型估计算法，计算时间减少了 $60\%$ 以上。但是，原始多尺度地理加权回归即便在 $N ≤ 5,000$ 时也很慢，因此对于更大的样本量可能需要更快的算法。

（2）贝叶斯空间变系数模型

贝叶斯空间变系数模型（ B-SVC ） 提供了另一种可用方法，参见 Gelfand 等，2003 年^[21]。该模型假设第 $k$ 个空间变系数 $\boldsymbol{\beta}_k = [\beta_k(\mathbf{s}_1)，\ldots, \beta_k(\mathbf{s}_N)]^T$ 服从以下多元高斯过程：

$$\beta_k \sim \mathcal{N}(b_k \mathbf{1}, \tau_k \mathbf{C}(r_k)) \tag{3}$$

其中 $b_k$、$\tau_k$ 和 $r_k$ 是参数，$\mathbf{1}$ 是一个 $N × 1$ 向量。 $\mathbf{C}(r_k)$ 是一个 $N × N$ 相关矩阵，其第 $(i, j)$ 个元素 $c_{ij}(r_k)$ 由距离衰减核给出，以捕获空间相关性，例如指数核 $c_{ij}(r_k) =\exp[-d(\mathbf{s}_i, \mathbf{s}_j)/r_k]$，其中 $d(\mathbf{s}_i, \mathbf{s}_j)$ 是测点 $\mathbf{s}_i$ 和 $\mathbf{s}_j$ 之间的欧氏距离。式 (1) 和 式 (2) 的贝叶斯空间变系数模型通过马尔可夫链蒙特卡罗 (MCMC) 技术估计。 Wheeler 和 Waller (2009)^[60]，以及 Finley 等 (2011a)^[16] 将贝叶斯空间变系数建模与地理加权回归进行了比较，并展示了贝叶斯空间变系数方法的稳健性。

不幸的是，当 $N$ 和/或 $K$ 很大时，贝叶斯空间变系数建模的计算成本令人失望（Finley 等，2011a ^[16]）。因为 MCMC 模拟需要在每次迭代中求密集矩阵 $\mathbf{C}(r_1),\ldots, \mathbf{C}(r_k)$ 的逆。为了减轻这种成本，Finley 等 (2011b)^[17] 对贝叶斯空间变系数模型应用了降阶。然而，当 $N = 8,774$ 、$K = 5$ 时，其 MCMC 模拟具有 $50,000$ 次迭代，仍然需要大约 72 小时。在实际工作中，即便更大的 $N$ 和 $K$ ，如果能够在几分钟内完成工作仍然大家愿意看到的。

（3）莫兰特征向量空间滤波法

Moran 特征向量空间变系数(M-SVC) 建模方法（Griffith，2008）^[27] 是 Moran 特征向量空间滤波方法 的扩展（参见 Griffith 2003 ^[25]；Griffith 和 Chun，2014 ^[29]），为多尺度空间变系数建模提供了第三种选择。虽然 Helbich 和 Griffith（2016 年）^[33] 以及 Oshan 和 Fotheringham（2018 年）^[50] 提到了这种方法的不稳定性，但 Murakami 等 (2017) ^[45] 将其扩展为近似贝叶斯空间变系数模型，Murakami 等 (2017 ^[45]; 2018 ^[43]) 通过蒙特卡罗实验证明了其估计准确性。

与贝叶斯空间变系数建模方法不同，该空间滤波方法不需要 MCMC 迭代，但当 $N$ 很大时仍然很慢，因为它需要特征分解（复杂度为 $\mathcal{O}(N^3)$ ）和 $2K$ 个参数的数值估计。

综上所述，如何在计算上有效地估计多尺度空间变系数仍然难以捉摸。但鉴于变系数模型的适用领域非常广泛且方法很流行，因此 不依赖高性能计算环境实现快速计算 仍然是大众所需的。

3 莫兰特征向量空间变系数(M-SVC) 建模

3.1 模型

该方法基于 Moran 系数（MC；Moran，1950 年）^[42]，它是空间依赖性的诊断统计量。$\mathbf{y} = [y(\mathbf{s}_1), \ldots, y(\mathbf{s}_N)]^T$ 的莫兰系数形式化为：

$$MC[\mathbf{y}] = \frac{N}{\mathbf{1}^T \mathbf{C} \mathbf{1}} \cdot \frac{\mathbf{y}^T (\mathbf{I} − \mathbf{1}\mathbf{1}^T/N) \mathbf{C}(\mathbf{I} − \mathbf{1}\mathbf{1}^T /N)\mathbf{y}}{ \mathbf{y}^T (\mathbf{I} − \mathbf{1}\mathbf{1}^T /N)\mathbf{y}} \tag{4}$$

其中 $\mathbf{C}$ 是对角线为零的对称空间邻接矩阵，$(\mathbf{I} − \mathbf{1}\mathbf{1}^T/N)$ 是实现居中化的矩阵。如果 $y(\mathbf{s}_1), \ldots, y(\mathbf{s}_N)$ 空间正相关时，$MC[\mathbf{y}] > −\frac{ 1 }{N−1}$，如果它们负相关，则 $MC[\mathbf{y}] < − \frac{1}{N−1}$ (Griffith, 2003)^[25]。可见邻接矩阵 $\mathbf{C}$ 对莫兰系数影响非常大。根据 Dray 等(2006) ^[13] 以及 Murakami 和 Griffith (2015) ^[44]，我们可以采用如下方法定义矩阵 $\mathbf{C}$：其第 $(i, j)$ 个元素为 $\exp(−d(\mathbf{s}_i, \mathbf{s}_j)/r)$，其中变程参数 $r$ 由所有样本点连接（或样本点对）之间的最小生成树中那个最大距离给出。与传统方法不同，本文中变程 $r$ 并不固定，我们引入了另一个参数来控制空间依赖性的 “尺度/变程” 比。这种替换是快速计算所必需的。

假设矩阵 $(\mathbf{I} − \mathbf{1}\mathbf{1}^T/N) \mathbf{C}(\mathbf{I} − \mathbf{1}\mathbf{1}^T /N)$ 的特征值为 $\boldsymbol{\lambda}=\lambda_1,\ldots,\lambda_N$ ，相应的特征向量为 $\mathbf{e}=\mathbf{e}_1,\ldots,\mathbf{e}_N$，令 $\mathbf{E} = [\mathbf{e}_1, \ldots , \mathbf{e}_L]$ 是由其中 $L (< N)$ 个特征向量构成的 $N × L$ 矩阵，而 $\boldsymbol{\Lambda}$ 是对角元素由相应的 $L$ 个特征值组成的 $L \times L$ 对角矩阵。那么下面的等式成立：

$$MC[\mathbf{e}_l] = \frac{N}{\mathbf{1}^T\mathbf{C} \mathbf{1}}λ_l$$

这意味着: 较大的正特征值对应的特征向量解释了较强的空间依赖性，而负空间依赖性则相反。也就是说， 特征向量提供了潜在空间依赖性的正交映射模式描述，其每一个量级都可以由莫兰系数值索引（Griffith，2003 年^[25]；Tiefelsdorf 和 Griffith，2007 年 ^[55]）。在现有空间变系数研究（例如，Fotheringham 等，2002 年^[18]）中，一般都隐含地假设了正空间依赖性（现实世界中负空间依赖性比较少见，一般出现在两个区域之间存在某种竞争的情况），因此本研究仅使用了其中 $L$ 个正特征值对应的 $(\mathbf{E}, \boldsymbol{\Lambda})$ 特征对。

Murakami 等 (2017) 的 M-SVC 模型指定如下：

$$\begin{align*} \mathbf{y} &= b_1 \mathbf{1} + \sum^K_{k=2}{\mathbf{x}_k} \circ \boldsymbol{\beta}_k + \mathbf{E} \boldsymbol{\gamma}_1 + \boldsymbol{\varepsilon} \tag{5}\\ \boldsymbol{\gamma}_1 &\sim \mathcal{N}(0, τ_1^2 \boldsymbol{\Lambda}^{α_1})\\ \varepsilon &\sim \mathcal{N}(0, \sigma^2I) \end{align*}$$

其中 “$\circ$” 表示列乘积运算符，$\mathbf{y}$ 是响应变量的 $N × 1$ 向量，$\mathbf{x}_k$ 是第 $k$ 个协变量的 $N × 1$ 向量。 $\mathbf{E} \boldsymbol{\gamma}_1$ 是捕获残差空间依赖性的项。 Murakami 和 Griffith (2018)^[43] 表明，$L = 200$ 的 $\mathbf{E} \boldsymbol{\gamma}_1$ 项可以大大降低残差的空间依赖性（意即可以捕获大部分残差的空间依赖性）。

第 $k$ 个变系数可以参数化如下：

$$\boldsymbol{\beta}_k = b_k \mathbf{1} + \mathbf{E} \boldsymbol{\gamma}_k,\, \boldsymbol{\gamma}_k \sim \mathcal{N}(0, \tau_k^2 \boldsymbol{\Lambda}^{α_k}) \tag{6}$$

其中 $b_k$ 是一个系数。空间变系数由常数项 $b_k \mathbf{1}$ 和空间变化项 $\mathbf{E} \boldsymbol{\gamma}_k$ 组成。$\mathbf{E} \boldsymbol{\gamma}_k$ 项具有以下性质：

• 性质（1）：可以根据莫兰系数进行解释；
• 性质（2）：特征向量均值为零。

性质(1) 对于准确建模空间依赖性的变化很重要，而 性质(2) 能使均值项 $b_k \mathbf{1}$ 与变化项 $\mathbf{E} \boldsymbol{\gamma}_k$ 之间可识别（参见 Gelfand 等，2003 年 ^[21]）。

向量 $\boldsymbol{\gamma}_k$ 中包含受收缩参数 $\tau_k^2$ 和 $α_k$ 控制的随机系数。

• 参数 $\tau_k^2$ 控制空间变系数 $\boldsymbol{\beta}_k$ 的方差
• 参数 $\alpha_k$ 控制空间尺度
  • 较大的 $\alpha_k$ 值会收缩具有小特征值的次要特征对上的系数，由此产生具有大尺度映射模式的空间变系数。
  • 较小的 $\alpha_k$ 则相反。

采用变化的 $\alpha$ 值可以极大地加速模型估计，同时保持估计尺度的灵活性。 Murakami 和 Griffith (2018)^[43] 表明，当估计 $\alpha_k$ 而不是 $r$ 时，所产生的估计误差非常小。

给定 $\mathbf{x}_1 = \mathbf{1}$ 和 $\widetilde{\mathbf{E}}_1 = (\mathbf{x}_1 \circ \mathbf{E}) = \mathbf{E}$，M-SVC 模型表示为:

$$\begin{align*} \mathbf{y} &= \mathbf{Xb} + \sum^{K}_{k=1} \widetilde{\mathbf{E}}_k \mathbf{V}(\boldsymbol{\theta}_k) \mathbf{u}_k + \boldsymbol{\varepsilon} \tag{7}\\ \mathbf{u}_k &\sim \mathcal{N}(0, \sigma^2 \mathbf{I}) \\ \boldsymbol{\varepsilon} &\sim \mathcal{N}(0, \sigma^2 \mathbf{I}) \end{align*}$$

其中 $\mathbf{X}=[\mathbf{x}_1, \ldots, \mathbf{x}_k]$, $\mathbf{b}=[b_1, \ldots, b_k]^T$, $\boldsymbol{\theta}_k \in \{ \tau_k , \alpha_k \}$, $\mathbf{V}(\boldsymbol{\theta}_k) = (\tau_k / \sigma) \boldsymbol{\Lambda}^{\alpha_k/2}$。式 (7) 表明 M-SVC 模型与线性混合效应模型相同（例如，Bates，2010）^[4]。

3.2.参数估计

第一类最大似然和第二类最大似然准则

（1）第一类最大似然

在统计与概率论中，估计统计模型参数的一个常用方法是使用最大似然估计，具体定义如下

$$\hat{\theta} \triangleq \mathop{\arg\max}\limits_{\theta}log \ p(\mathcal D |\theta)$$

上述最大似然估计对应的似然通常称之为第一类最大似然。

（2）第二类最大似然

当使用贝叶斯方法时，可以通过数值优化方法先获得超参数的经验最优值，具体如下：

$$\lambda^{*}=\mathop{\arg\max} \limits_{\lambda}p(D| \lambda)$$

其中 $\lambda$ 是超参数，此类方法也被称为 经验贝叶斯方法 或者 第二类最大似然估计方法，与交叉验证法评估不同超参数的方法相比，此方法更有效。

M-SVC 模型通过最大化第二类似然（经验贝叶斯）方法进行估计，即最大化

$$\text{loglik}_R(\boldsymbol{\Theta})=\int \text{loglik} (\mathbf{y} \mid \mathbf{b},\boldsymbol{\Theta}) d \mathbf{b}$$

其中 $\boldsymbol{\Theta} \in \{ \boldsymbol{\theta}_1, \ldots,\boldsymbol{\theta}_K \}$， $\text{loglik} (\mathbf{y} \mid \mathbf{b}, \boldsymbol{\Theta}) = \int \ldots \int p(\mathbf{y}, \mathbf{u}_1, \ldots, \mathbf{u}_K \mid \mathbf{b},\boldsymbol{\Theta}) p(\mathbf{u}_1, \ldots, \mathbf{u}_K)d \mathbf{u}_1, \ldots, d \mathbf{u}_K$ 。

如果将 $p(\mathbf{y}, \mathbf{u}_1, \ldots, \mathbf{u}_K \mid \mathbf{b}, \boldsymbol{\Theta})$ 和 $p(\mathbf{u}_1, \ldots, \mathbf{u}_K)$ 定义为高斯密度函数，则此似然具有解析表达式。使用此性质，可以推导出 式 (7) 的对数似然为（参见 Bates，2010 ^[4]；Murakami 等，2017 ^[45]）：

$$\text{loglik}_R(\boldsymbol{\Theta}) = − \frac{1}{2} \ln |P| − \frac{N-K}{2} \left(1 + \ln \left( \frac{2πd(\boldsymbol{\Theta})}{N − K} \right)\right) \tag{8}$$

$$\mathbf{P}=\left[\begin{array}{cccc} \mathbf{X}^{T} \mathbf{X} & \mathbf{X}^{T} \widetilde{\mathbf{E}}_1 \mathbf{V}\left(\boldsymbol{\theta}_1\right) & \cdots & \mathbf{X}^{T} \widetilde{\mathbf{E}}_K \mathbf{V}\left(\boldsymbol{\theta}_K\right) \\ \mathbf{V}\left(\boldsymbol{\theta}_1\right) \widetilde{\mathbf{E}}_1^{T} \mathbf{X} & \mathbf{V}\left(\boldsymbol{\theta}_1\right) \widetilde{\mathbf{E}}_1^{T} \widetilde{\mathbf{E}}_1 \mathbf{V}\left(\boldsymbol{\theta}_1\right)+\mathbf{I} & \cdots & \mathbf{V}\left(\boldsymbol{\theta}_1\right) \widetilde{\mathbf{E}}_1^{T} \widetilde{\mathbf{E}}_K \mathbf{V}\left(\boldsymbol{\theta}_K\right) \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{V}\left(\boldsymbol{\theta}_K\right) \widetilde{\mathbf{E}}_K^{T} \mathbf{X} & \mathbf{V}\left(\boldsymbol{\theta}_K\right) \widetilde{\mathbf{E}}_K^{T} \widetilde{\mathbf{E}}_1 \mathbf{V}\left(\boldsymbol{\theta}_1\right) & \cdots & \mathbf{V}\left(\boldsymbol{\theta}_K\right) \widetilde{\mathbf{E}}_K^{T} \widetilde{\mathbf{E}}_K \mathbf{V}\left(\boldsymbol{\theta}_K\right)+\mathbf{I} \end{array}\right] \tag{9}$$

导数 $d(\boldsymbol{\Theta})$ 等于:

$$\boldsymbol{d}(\boldsymbol{\theta})=\left\|\mathbf{y}-\mathbf{X} \hat{\mathbf{b}}-\sum_{k=1}^K \widetilde{\mathbf{E}}_k \mathbf{V}\left(\boldsymbol{\theta}_k\right) \widehat{\mathbf{u}}_k\right\|^2+\sum_{k=1}^K\left\|\widehat{\mathbf{u}}_k\right\|^2 \tag{10}$$

其中

$$\left[\begin{array}{c} \hat{\mathbf{b}} \\ \widehat{\mathbf{u}}_1 \\ \vdots \\ \widehat{\mathbf{u}}_K \end{array}\right]=\mathbf{P}^{-1}\left[\begin{array}{c} \mathbf{X}^{T} \mathbf{y} \\ \mathbf{V}\left(\boldsymbol{\theta}_1\right) \widetilde{\mathbf{E}}_1^{T} \mathbf{y} \\ \vdots \\ \mathbf{V}\left(\boldsymbol{\theta}_K\right) \widetilde{\mathbf{E}}^{T}{ }_K \mathbf{y} \end{array}\right] \tag{11}$$

式 (10) 平衡了模型准确性和复杂性之间的权衡。可以通过以下步骤估计参数：

• 步骤(1)： 通过最大化 式（8） 估计 $\boldsymbol{\Theta} \in \{ \boldsymbol{\theta}_1, \ldots, \boldsymbol{\theta}_K \}$
• 步骤(2)： 通过 式（11） 估计 $\{\mathbf{b}, \mathbf{u}_1, \ldots, \mathbf{u}_K \}$
• 步骤(3)： 通过$\hat{\beta}_k =\mathbf{b}_k \mathbf{1} + \mathbf{E} \hat{\boldsymbol{\gamma}}_k =\hat{b}_k \mathbf{1} + \mathbf{EV}(\hat{\boldsymbol{\theta}}_k) \hat{\mathbf{u}}_k$ 估计空间变系数。

4. 快速 M-SVC 建模

不幸的是，Moran 特征向量方法在建模和估计方面的计算效率都很低，在本节中将讨论加速的方法。第 4.1 节 介绍了如何降低建模成本，而 第 4.2 节 介绍了如何降低估计成本。

4.1 建模

$(\mathbf{I} − \mathbf{1}\mathbf{1}^T/N) \mathbf{C}(\mathbf{I} − \mathbf{1}\mathbf{1}^T /N)$ 的特征分解成本为 $\mathcal{O}(N^3)$ ，这对于大 N 来说是难以处理的。此外，如果使用距离核给出 $\mathbf{C}$，则必须在分解之前存储 $N × N$ 的矩阵。就计算复杂性和内存使用而言，效率过低。

为了解决此问题，我们使用了基于 Nystrom 扩展（Drineas 和 Mahoney，2005 年 ^[14]）的近似特征分解方法，求得 Moran 特征向量和特征值的近似（Murakami 和 Griffith 2018 ^[43]），公式如下：

$$\hat{\mathbf{E}}=\left[\mathbf{C}_{N L}-\mathbf{1} \otimes \frac{\mathbf{1}_L^{\prime}\left(\mathbf{C}_L+\mathbf{I}_L\right)}{L}\right] \mathbf{E}_L\left(\boldsymbol{\Lambda}_L+\mathbf{I}_L\right)^{-1} \tag{12}$$

$$\hat{\boldsymbol{\Lambda}}=\frac{L+N}{L}\left(\boldsymbol{\Lambda}_L+\mathbf{I}_L\right)-\mathbf{I}_L \tag{13}$$

其中 “$\otimes$” 表示 Kronecker 积运算符。$\mathbf{C}_L$ 是一个 $L × L$ 空间邻接矩阵 ($L \ll N$) ，其中 $L$ 个结点分布在整个研究区域。

结点的坐标由 k 均值聚类 的簇中心坐标定义，可以簇内样本的空间坐标计算得出（参见 Zhang 和 Kwok，2010 ^[62]）。 $\mathbf{C}_{NL}$ 是 $L$ 个结点和 $N$ 个样本点之间的 $N × L$ 邻近矩阵。集合 $\{ \mathbf{E}_L, \boldsymbol{\Lambda}_L, \mathbf{I}_L, \mathbf{1}_L\}$ 的定义就像是所有结点的 $\{ \mathbf{E}, \boldsymbol{\Lambda}, \mathbf{I}, \mathbf{1} \}$ 集合。

根据继 Murakami 和 Griffith (2018) ^[43] ，采用 $200$ 个特征对足以模拟空间依赖性。为此，如果 $L$ 低于 $200$ 时，则 $L$ 固定为 $\hat{\boldsymbol{\Gamma}}$ 中正特征值的数量；否则将 $L = 200$。不管哪种情况，$L$ 都是固定的，特征近似的计算复杂度减少为 $\mathcal{O}(N)$，而内存使用量为 $N × L$。关于空间变系数的数量，考虑到识别难度较大（参见 Wheeler 和 Tiefelsdorf，2005 年^[59]；Farber 和 Páez，2007 年^[15]；Helbich 和 Griffith，2016 年^[33]），通常至多设置 $10$ 个空间变系数可能是合理的。

根据上述内容，我们假设 $L$ 是固定的且 $K \ll N$。当必须考虑大量解释变量时，只能针对某些重点选定的解释变量给出空间变系数，而其他变量仍然假定为常数系数。此时，第 $k$ 个系数通过用 $b_k \mathbf{1}$ 替换掉 $\beta_k = b_k \mathbf{1} + \mathbf{E}\boldsymbol{\gamma}_k$ 来转换成常数系数。此性质有助于避免增加 $K$ 并确保空间变系数的可识别性。

我们的近似 M-SVC 模型是使用 式(14) 和 式 (15)：

$$\begin{align*} \mathbf{y} &=b_1 \mathbf{1}+\sum_{k=2}^K \mathbf{x}_k \circ \boldsymbol{\beta}_k + \hat{\mathbf{E}} \boldsymbol{y}_1+\boldsymbol{\varepsilon} \tag{14}\\ \boldsymbol{\gamma}_1 &\sim N\left(\mathbf{0}, \tau_1^2 \widehat{\boldsymbol{\Lambda}}^{\alpha_1}\right)\\ \boldsymbol{\varepsilon} &\sim N\left(\mathbf{0}, \sigma^2 \mathbf{I}\right) \end{align*}$$

$$\begin{align*} \boldsymbol{\beta}_k &= b_k \mathbf{1}+\hat{\mathbf{E}} \boldsymbol{\gamma}_k \tag{15}\\ \quad \boldsymbol{\gamma}_k &\sim N\left(\mathbf{0}, \tau_k^2 \widehat{\boldsymbol{\Lambda}}^{\alpha_k}\right) \end{align*}$$

下一节将解释如何计算有效地对模型进行估计。

4.2 估计

似然最大化的计算过程可能非常慢，因为在 $\boldsymbol{\Theta} \in \{ \boldsymbol{\theta}_1, \ldots, \boldsymbol{\theta}_K \}$ 中包含 $2K$ 个收缩参数，其中 $\boldsymbol{\theta}_k \in \{ \tau_k , \alpha_k \}$ 没有封闭形式解。例如，如果假设有 $10$ 个空间变系数，则必须对 $20$ 个收缩参数进行数值估计。 M-SVC 模型的数值估计比地理加权回归困难很多，因为地理加权回归只有一个带宽参数或 $K$ 个不同尺度下的带宽参数。

本节表明，通过采用选定的矩阵技巧，M-SVC 模型估计的计算成本大大降低。第 4.2.1 节 从似然函数中消除了 $N$，第 4.2.2 节 推导出一种顺序方法来最大化似然。

4.2.1 从似然函数中消除“N”

要估计 $\boldsymbol{\Theta} \in \{\boldsymbol{\theta}_1, \ldots,\boldsymbol{\theta}_K\}$，必须重复计算受限对数似然。为了减轻这种负担，本研究在估计之前，消去了那些大小取决于 $N$ 的矩阵和向量。具体来说，似然的计算有以下几个步骤：

反复。为了减轻这种负担，本研究在估计之前消除了大小取决于 N 的矩阵和向量。具体来说，似然度的评估有以下几个步骤：

• 步骤 1: 计算 $\mathbf{M}_{0,0}=\mathbf{X}^TX$、$\mathbf{M}_{0,k} = \mathbf{X}^T (\mathbf{x}_k \circ \mathbf{E})$, $\mathbf{M}_{k,\widetilde{i}} = (\mathbf{x}_k \circ \mathbf{E})^T (\mathbf{x}_{\widetilde{k}} \circ \mathbf{E})$、 $m_0=\mathbf{X}^T \mathbf{y}$、 $\mathbf{m}_k = ( \mathbf{x}_k \circ \mathbf{E})^T \mathbf{y}$ 和 $\mathbf{m}_{y,y} = \mathbf{y}^T \mathbf{y}$；

• 步骤 2: 将矩阵和向量代入 式(8) 来重写受限对数似然，如下所示：

$$\log \operatorname{lik}_R(\boldsymbol{\Theta})=-\frac{1}{2} \ln |\mathbf{P}|-\frac{N-K}{2}\left(1+\ln \left(\frac{2 \pi}{N-K} \boldsymbol{d}(\boldsymbol{\Theta})\right)\right) \tag{16}$$

其中

$$\mathbf{P}=\left[\begin{array}{cccc} \mathbf{M}_{0,0} & \mathbf{M}_{0,1} \mathbf{V}\left(\boldsymbol{\theta}_1\right) & \cdots & \mathbf{M}_{0, K} \mathbf{V}\left(\boldsymbol{\theta}_K\right) \\ \mathbf{V}\left(\boldsymbol{\theta}_1\right) \mathbf{M}_{1,0} & \mathbf{V}\left(\boldsymbol{\theta}_1\right) \mathbf{M}_{1,1} \mathbf{V}\left(\boldsymbol{\theta}_1\right)+\mathbf{I} & \cdots & \mathbf{V}\left(\boldsymbol{\theta}_1\right) \mathbf{M}_{1, K} \mathbf{V}\left(\boldsymbol{\theta}_K\right) \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{V}\left(\boldsymbol{\theta}_K\right) \mathbf{M}_{K, 0} & \mathbf{V}\left(\boldsymbol{\theta}_K\right) \mathbf{M}_{K, 1} \mathbf{V}\left(\boldsymbol{\theta}_1\right) & \cdots & \mathbf{V}\left(\boldsymbol{\theta}_K\right) \mathbf{M}_{K, K} \mathbf{V}\left(\boldsymbol{\theta}_K\right)+\mathbf{I} \end{array}\right] \tag{17}$$

并且 $\boldsymbol{d}(\boldsymbol{\Theta})=\|\hat{\boldsymbol{\varepsilon}}\|^2+\sum_{k=1}^K\left\|\hat{\mathbf{u}}_k\right\|^2$，其中：

$$\|\hat{\boldsymbol{\varepsilon}}\|^2=m_{y y}-2\left[\hat{\mathbf{b}}^{T}, \widehat{\mathbf{u}}_1^{T}, \cdots \widehat{\mathbf{u}}_K^{T}\right]\left[\begin{array}{c} \mathbf{m}_0 \\ \mathbf{V}\left(\boldsymbol{\theta}_1\right) \mathbf{m}_1 \\ \vdots \\ \mathbf{V}\left(\boldsymbol{\theta}_K\right) \mathbf{m}_K \end{array}\right]+\left[\hat{\mathbf{b}}^{T}, \widehat{\mathbf{u}}_1^{T}, \cdots \widehat{\mathbf{u}}_K^{T}\right] \mathbf{P}_0\left[\begin{array}{c} \hat{\mathbf{b}} \\ \widehat{\mathbf{u}}_1 \\ \vdots \\ \widehat{\mathbf{u}}_K \end{array}\right] \tag{18}$$

$$\left[\begin{array}{c} \hat{\mathbf{b}}\\ \hat{\mathbf{u}}_1 \\ \vdots \\ \widehat{\mathbf{u}}_K \end{array}\right]=\mathbf{P}^{-1}\left[\begin{array}{c} \mathbf{m}_0 \\ \mathbf{V}\left(\boldsymbol{\theta}_1\right) \mathbf{m}_1 \\ \vdots \\ \mathbf{V}\left(\boldsymbol{\theta}_K\right) \mathbf{m}_K \end{array}\right] \tag{19}$$

$\mathbf{P}_0=\mathbf{P}$, 其中 $\mathbf{V}\left(\boldsymbol{\theta}_k\right) \mathbf{M}_{k, k} \mathbf{V}\left(\boldsymbol{\theta}_k\right)+\mathbf{I}$ 被替换为 $\mathbf{V}\left(\boldsymbol{\theta}_k\right) \mathbf{M}_{k, k} \mathbf{V}\left(\boldsymbol{\theta}_k\right)$。

• 步骤 3: 通过最大化 式 (16) 来估计$\boldsymbol{\Theta}$

因为 式(16) 不包括任何大小取决于 $N$ 的矩阵，用于估计 $\boldsymbol{\Theta}$ 的计算复杂度与 $N$ 无关。然而，如果 $K$ 很大，则 步骤 (c) 中的最大化速度太慢。为了减少这种计算成本，下一节开发了一种方法来按顺序估计每个 $\boldsymbol{\theta}_k$。

4.2.2 顺序似然最大化

本节讨论通过对每个 $k$ 依次求解 式(20) 直到收敛使得似然值 $\text{loglike}_R(\boldsymbol{\Theta})$ 最大化：

$$\boldsymbol{\theta}_k = \arg \max_{\boldsymbol{\theta}_k} \text{loglike}_R(\boldsymbol{\theta}_k \mid \boldsymbol{\Theta}_{−k}) \tag{20}$$

其中 $\boldsymbol{\Theta}_{-k} \in \{\boldsymbol{\theta}_1, \ldots, \boldsymbol{\theta}_{k-1}, \boldsymbol{\theta}_{k+1}, \ldots \boldsymbol{\theta}_K \}$。

这种最大化需要评估 $\mathbf{P}^{-1}$ 和 $|\mathbf{P} |$，它们的复杂度都是 $\mathcal{O}((K + KL)^3)$；随着 $K$ 的增长，复杂性迅速增加，而 $L$ 被限制为最多 $200$。由于必须多次评估 $\text{loglike}_R(\boldsymbol{\theta}_k \mid \boldsymbol{\Theta}^{−k})$ 以使其最大化，因此计算负担不可忽略。

为了降低计算成本，第 4.2.2.1 节 和 第 4.2.2.2 节 分别介绍了计算 $\mathbf{P}^{-1}$ 和 $|\mathbf{P}|$ 的快速方法。尽管此部分假设以估计 $\boldsymbol{\theta}_K$ 为例，但同样的方法也可用于估计 $\{\boldsymbol{\theta}_1, \ldots, \boldsymbol{\theta}_{K-1}\}$。

（1）快速计算包含 $\mathbf{P}^{-1}$ 的项

假设 $\mathbf{m}_{-K}=\left[\begin{array}{llll}\mathbf{m}_0^{T} & \mathbf{m}_1^{T} & \cdots & \mathbf{m}_{K-1}^{T}\end{array}\right]^{T}$ 。根据附录 1，包括逆 $\mathbf{P}^{-1}$ 的 式（9） 可以表示为：

$$\begin{aligned} & {\left[\begin{array}{c} \widehat{\mathbf{b}} \\ \widehat{\mathbf{u}}_1 \\ \vdots \\ \widehat{\mathbf{u}}_K \end{array}\right]=\left[\begin{array}{cc} \widetilde{\mathbf{v}}_{-K}^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{V}\left(\boldsymbol{\theta}_K\right)^{-1} \end{array}\right] \mathbf{Q}^{-1}\left[\begin{array}{c} \mathbf{m}_{-K} \\ \mathbf{m}_K \end{array}\right]} \\ & -\left[\begin{array}{c} \widetilde{\mathbf{V}}_{-K}^{-1} \mathbf{Q}_{-K, K}^* \\ \mathbf{V}\left(\boldsymbol{\theta}_K\right)^{-1} \mathbf{Q}_{K, K}^* \end{array}\right]\left(\mathbf{V}\left(\boldsymbol{\theta}_K\right)^2+\mathbf{Q}_{K, K}^*\right)^{-1}\left[\mathbf{Q}_{K,-K}^* \mathbf{m}_{-K}+\mathbf{Q}_{K, K}^* \mathbf{m}_K\right] \tag{21} \\ \end{aligned}$$

其中:

$$\widetilde{\mathbf{V}}_{-K}=\left[\begin{array}{llll}\mathbf{I} & & & \\ & \mathbf{V}_1 & & \\ & & \ddots & \\ & & & \mathbf{V}_{K-1}\end{array}\right]$$

$$\mathbf{Q}=\left[\begin{array}{cc}\widetilde{\mathbf{M}}_{-K,-K}+\widetilde{\mathbf{V}}_{-K}^{-2} & \widetilde{\mathbf{M}}_{-K, K} \\ \widetilde{\mathbf{M}}_{K,-K} & \mathbf{M}_{K, K}\end{array}\right]$$

并且：

$$\widetilde{\mathbf{M}}_{-K,-K}= \left[\begin{array}{ccccc}\mathbf{M}_{0,0} & \mathbf{M}_{0,1} & & \cdots & \mathbf{M}_{0, K-1} \\ \mathbf{M}_{1,0} & \mathbf{M}_{1,1}+\mathbf{V}_1^{-2} & & \cdots & \mathbf{M}_{1, K-1} \\ \vdots & \vdots & \ddots & & \vdots \\ \mathbf{M}_{K-1,0} & \mathbf{M}_{K-1,1} & \cdots & \mathbf{M}_{K-1, K-1}+\mathbf{V}_{K-1}^{-2}\end{array}\right]$$

$$\widetilde{\mathbf{M}}_{-K, K}=\left[\begin{array}{llll}\mathbf{M}_{K, 0} & \mathbf{M}_{K, 1} & \cdots & \mathbf{M}_{K, K-1}\end{array}\right]^{\prime}$$

$$\left[\begin{array}{cc}\mathbf{Q}_{-K,-K}^* & \mathbf{Q}_{-K, K}^* \\ \mathbf{Q}_{K,-K}^* & \mathbf{Q}_{K, K}^*\end{array}\right]=\mathbf{Q}^{-1}$$

不同于 式（19），式(21) 不包含对大矩阵 $\boldsymbol{\theta}_K$ 的求逆运算。式(21) 依然有 $\mathbf{Q}^{-1}$，但其复杂度为 $\mathcal{O}((K+KL)^3)$。由于 $\mathbf{Q}$ 矩阵不包含$\boldsymbol{\theta}_K$，如果只对 $\mathbf{Q}^{-1}$ 进行一次求值，则可以在迭代计算中将其固定，使 $\text{loglike}_R(\boldsymbol{\theta}_k \mid \boldsymbol{\Theta}_{−k})$ 最大化。迭代计算中计算量最大的部分是 $(\mathbf{V}(\boldsymbol{\theta}_K)^2 + \mathbf{Q}^*_{K,K})^{−1}$， 其计算成本为 $\mathcal{O}(L^3)$，在 $L ≤ 200$ 约束下这是微不足道的。

（2）快速计算 $|\mathbf{P}|$

如 附录 2 所示，$|\mathbf{P}|$ 可以扩展如下：

$$\begin{gathered} |\mathbf{P}|=\left|\widetilde{\mathbf{V}}_{-K}\right|^2\left|\mathbf{V}\left(\boldsymbol{\theta}_K\right)\right|^2\left|\widetilde{\mathbf{M}}_{-K,-K}+\widetilde{\mathbf{V}}_{-K}^{-2}\right| \mid \mathbf{V}\left(\boldsymbol{\theta}_K\right)^{-2}+\mathbf{M}_{K, K} \\ -\widetilde{\mathbf{M}}_{K,-K}\left(\widetilde{\mathbf{M}}_{-K,-K}+\widetilde{\mathbf{V}}_{-K}^{-2}\right)^{-1} \widetilde{\mathbf{M}}_{-K, K} \mid . \end{gathered} \tag{22}$$

$|\widetilde{\mathbf{V}}_{-K}|^2$ 和 $|\mathbf{V}\left(\boldsymbol{\theta}_K\right)|^2$ 的计算是即时的，因为 $\widetilde{\mathbf{V}}_{-K}$ 和 $\mathbf{V}\left(\boldsymbol{\theta}_K\right)$ 是对角矩阵。计算 $\left|\widetilde{\mathbf{M}}_{-K,-K}+\widetilde{\mathbf{V}}_{-K}^{-2}\right|$ 和 $\left(\widetilde{\mathbf{M}}_{-K,-K}+\widetilde{\mathbf{V}}_{-K}^{-2}\right)^{-1}$ 的复杂度是 $\mathcal{O}\left((K+(K-1) L)^3\right)$， 比较大。然而，由于这些矩阵不包括 $\boldsymbol{\theta}_K$，因此可以在参数估计之前对其进行评估。毕竟似然最大化中必须迭代的最苛刻的计算是 $|\mathbf{V}(\boldsymbol{\theta}_K)^{−2} + \widetilde{\mathbf{M}}_{K,K}|$, 其中 $\widetilde{\mathbf{M}}_{K,K} = \mathbf{M}_{K,K} − \widetilde{\mathbf{M}}_{K,-K}(\widetilde{\mathbf{M}}_{-K,-K} + \widetilde{\mathbf{V}}_{−K}^{−2})^{−1} \widetilde{\mathbf{M}}_{-K,K}$ 提前给出。此行列式的计算复杂度为 $O(L^3)$。

总之，如果提前处理不包含 $\boldsymbol{\theta}_K$ 的矩阵和向量，则包含在函数 $\text{loglike}_R(\boldsymbol{\theta}_k\mid\boldsymbol{\Theta}_{−k})$ 中的 $\left|\widetilde{\mathbf{M}}_{-K,-K}+\widetilde{\mathbf{V}}_{-K}^{-2}\right|$ 和 $\left(\widetilde{\mathbf{M}}_{-K,-K}+\widetilde{\mathbf{V}}_{-K}^{-2}\right)^{-1}$ 的计算复杂度都会变为 $O(L^3)$。

4.3 总结

本文方法总结在 图 1 中。在建模步骤中，我们应用

• (1) 降阶，每个空间变系数表示为 $L$ 个近似特征向量 $\hat{\mathbf{E}}$ （$N × L$）的线性组合。在估计步骤中，我们首先应用

• (2) 预压缩，空间变系数模型被压缩成 $K^2(L + 1)^2 + K +1$ 个内积，见 第 4.2.1 节 的步骤 (1) 中。由于这一步，似然评估的计算成本与 $N$ 无关。

• (3) 顺序估计加速了似然最大化。 步骤 (3) 中的估计迭代处理 $\mathcal{O}(L^3)$ 。因为成本与 $N$ 无关，所以即使 $N$ 非常大，估计也很快。

  然而，计算成本随着 $K$ 的增长而增加。这个结果是因为大 $K$ 增加了迭代次数来依次估计 $\boldsymbol{\theta}_1, \ldots, \boldsymbol{\theta}_K$。如前所述，可以通过仅将空间变系数应用于某些重点解释变量来降低计算成本。

图 1：本文所提出方法的概述

5 蒙特卡洛模拟实验

本节将本文方法与现有方法在计算时间和估计精度方面进行比较。我们部署了三个模拟实验。第一个实验将本文方法与假设 N ≤ 4,000 的原始 M-SVC、地理加权回归和 FBGWR 方法进行比较。第二个模拟比较（相对）快速方法，包括地理加权回归和 M-SVC，假设 6,000 ≤ N ≤ 12,000。第三个实验检查本文方法在更大样本（20,000 ≤ N ≤ 100,000）下的性能。在本节中，我们使用 Mac Pro（3.5 GHz、6 核 Intel Xeon E5 处理器和 64 GB 内存）。 R（版本 3.6.2；https://cran.r-project.org/）用于模型估计。 GWmodel 包（版本 2-0.5；参见 Lu et al., 2014）用于估计地理加权回归和多尺度地理加权回归模型。

5.1 使用相对较小的样本进行模拟 (N ≤ 4,000)

6 实证研究

作为说明，地理加权回归和 M-SVC(3) 应用于 2010 年东京都市区官方评估的住宅地价数据（$N = 7,679$）。响应变量是对数土地价格 [JPY/m2]。协变量是到最近的火车站的距离 (Station_d) [km]、到东京火车站的距离 (Tokyo_d) [km]、1 公里网格中的绿地份额 (Green) 和预期的洪水深度 (Flood) [m]。这些数据来自日本 MLIT 提供的国家土地数值信息下载服务 (http://nlftp.mlit.go.jp/ksj-e/index.html)。

图 15 绘制了估计的系数（空间变化截距除外）。地理加权回归的估计用了 $721$ 秒，而 M-SVC 方法用了 $241$ 秒。使用地理加权回归估计的空间变系数具有相似的空间模式尺度。相比之下，使用 M-SVC 方法估计的空间变系数的空间模式具有相当不同的尺度；具体而言，Tokyo_d 的空间变系数具有大尺度空间格局，Forest 的空间变系数具有中等尺度的格局，Station_d 和 Flood 的空间变系数具有小尺度的格局。 Tokyo_d 的空间变系数在东京市中心和大都市区北部有很大的负值；这表明东京的可达性促进了城市化，尤其是在北部。 Station_d 的系数表明，车站可达性与郊区住宅区的土地价格具有更强的负相关关系。由于 Tokyo_d 和 Station_d 可被视为可访问性的全局和局部度量，因此这些空间变系数估计中的大尺度和小尺度空间模式在直觉上是合理的。森林的空间变系数表明绿地在郊区非城市地区具有很强的负相关关系。最后，洪水的空间变系数在荒川河周围有更强的负相关关系，荒川河是该地区的主要河流。这表明洪水风险已适当反映在荒川沿海地区的土地价格中。相比之下，地理加权回归估计似乎更难解释，因为它们在空间变系数中的规模相同。

总之，M-SVC 方法提供了多尺度空间变系数的合理估计。

图 15：估计的 SVC。洪水子图中显示了主要河道。

7 总结

虽然大空间数据建模是近期的热门话题，但涉及变系数建模的相关讨论却相当有限。鉴于此背景，本研究开发了一种快速 M-SVC 方法来估计多尺度空间变系数模型。我们使用 (1) 降阶、(2) 预压缩和 (3) 顺序估计来实现计算效率。蒙特卡洛模拟实验证实了我们方法的估计精度和计算效率。本文方法在相对较短的时间内估计空间变系数模型，无需应用任何并行计算。换句话说，本文方法不需要高性能计算环境。

（1）可并行性

然而，步骤 (1) 很容易并行化，类似于其他基于 Nystrom 扩展的特征近似方法（参见 Li 等，2015 ^[37]）。步骤（2） 也是如此，它只是简单地计算内积。尽管并行化顺序估计 步骤 (3) 可能比较困难，但其计算成本已经与 $N$ 无关。因此，我们方法的并行化允许更快地估计多尺度空间变系数模型，并将此处总结的结果扩展到数百万个样本。

（2）缺点

本文方法的一个缺点是它无法对精细尺度的空间变化进行建模，如 图 13 所示；我们最多只考虑了 200 个对应于大特征值的特征向量，这意味着忽略了其他解释小尺度变化的特征向量。地统计学中报告了低秩空间建模的这种局限性（Stein，2014 ^[54]）。解决该问题是未来一个重要的研究课题。幸运的是，已经针对这个问题提出了许多多尺度方法（例如，Sang 和 Huang，2012 年^[51]；Nychika 等，2015 年 ^[49]）。

局部方法比降秩方法（包括本文方法）更适合捕获精细尺度的变化（Stein，2014 ^[54]）。在这方面，加速的地理加权回归是估计具有大样本的精细尺度空间变系数的明智方法。考虑到 M-SVC 和多尺度地理加权回归方法的相似估计精度（参见 Murakami 等，2018 年^[43]），加速多尺度地理加权回归、具有参数特定距离矩阵的地理加权回归（PSDMGWR；Lu 等，2017 年 ^[40]），以及其他扩展的地理加权回归，也是开发快速灵活的空间变系数建模方法的有前途的方法。关于快速地理加权回归的研究包括 Harris (2010b) 等^[31]，Tran 等 (2016) ^[57]，和 Lu 等 （2017 年 ^[40]；2018 年 ^[41]）。

虽然我们开发了一种快速的空间变系数估计方法，但通过写成 式 (7)，同样的方法可能适用于其他混合效应模型，如空间（或非空间）加性模型和多级模型等。将我们的快速估计框架扩展到各种空间回归模型在大空间数据时代将很有价值。

M-SVC 建模方法在 R 包“spmoran”中实现（https://cran.r-project.org/web/packages/spmoran/index.html；参见 Murakami，2018 ^[43]）

参考文献

• [1] Anselin, L., 2010. Thirty years of spatial econometrics. Pap. Reg. Sci. 89(1), 3–25.
• [2] Arbia, G., 2014. Pairwise likelihood inference for spatial regressions estimated on very large datasets. Spat. Stat. 7, 21–39.
• [3] Banerjee, S., Gelfand, A. E., Finley, A. O., Sang, H., 2008. Gaussian predictive process models for large spatial data sets. J. R. Stat. Soc. Series B Stat. Methodol. 70 (4), 825–848.
• [4] Bates, D. M., 2010. lme4: Mixed-effects modeling with R. http://lme4.r-forge.rproject.org/book.
• [5] Burden, S., Cressie, N., Steel, D.G., 2015. The SAR model for very large datasets: a reduced rank approach. Econom. 3 (2), 317–338.
• [6] Brunsdon, C., Fotheringham, A.S., Charlton, M.E., 1996. Geographically weighted regression: a method for exploring spatial nonstationarity. Geogr. Anal. 28 (4), 281298.
• [7] Cahill, M., Mulligan, G., 2007. Using geographically weighted regression to explore local crime patterns. Soc. Sci. Comput. Rev. 25 (2), 174–193.
• [8] Cressie, N., 1993. Statistics for Spatial Data. John Wiley & Sons, New York.
• [9] Cressie, N., Johannesson, G., 2008. Fixed rank kriging for very large spatial data sets. J. R. Stat. Soc. Series B Stat. Methodol. 70 (1), 209–226.
• [10] Datta, A., Banerjee, S., Finley, A.O., Gelfand, A.E., 2016. Hierarchical nearestneighbor Gaussian process models for large geostatistical datasets. J. Am. Stat. Assoc. 111 (514), 800–812.
• [11] Debarsy. N. Yang, Y., 2018. Editorial for the special issue entitled: New advances in spatial econometrics: Interactions matter. Reg. Sci. Urban Econ. DOI: 10.1016/j.regsciurbeco.2018.02.004.
• [12] Dong, G., Nakaya, T., Brunsdon, C., 2018. Geographically weighted regression models for ordinal categorical response variables: An application to geo-referenced life satisfaction data. Comput. Environ. Urban Syst. 70, 35–42.
• [13] Dray, S., Legendre, P., PeresNeto, P.R., 2006. Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecol. Model. 196 (3–4), 483–493.
• [14] Drineas, P., Mahoney, M.W., 2005. On the Nyström method for approximating a Gram matrix for improved kernel-based learning. J. Mach. Learn. Res. 6, 2153–2175.
• [15] Farber, S., Páez, A., 2007. A systematic investigation of cross-validation in GWR model estimation: empirical analysis and Monte Carlo simulations. J. Geog. Sci. 9 (4), 371–396.
• [16] Finley, A.O., 2011. Comparing spatially-varying coefficients models for analysis of ecological data with non–stationary and anisotropic residual dependence. Methods Ecol. Evol. 2 (2), 143–154.
• [17] Finley, A.O., Banerjee, S., MacFarlane, D.W., 2011. A hierarchical model for quantifying forest variables over large heterogeneous landscapes with uncertain forest areas. J. Am. Stat. Assoc. 106 (493), 31–48.
• [18] Fotheringham, A.S., Brunsdon, C., Charlton, M., 2002. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. John Wiley & Sons Chichester, UK.
• [19] Fotheringham, A.S., Yang, W., Kang, W., 2017. Multiscale Geographically Weighted Regression (MGWR). Ann. Am. Assoc. Geogr. 107 (6), 1247–1265.
• [20] Furrer, R., Genton, M.G., Nychka, D., 2006. Covariance tapering for interpolation of large spatial datasets. J. Comput. Graph. Stat. 15 (3), 502–523.
• [21] Gelfand, A.E., Kim, H.J., Sirmans, C.F., Banerjee, S., 2003. Spatial modeling with spatially varying coefficient processes. J. Am. Stat. Assoc. 98 (462), 378–396.
• [22] Geniaux, G., Martinetti, D., 2018. A new method for dealing simultaneously with spatial autocorrelation and spatial heterogeneity in regression models. Reg Sci Urban Econ. DOI: 10.1016/j.regsciurbeco.2017.04.001.
• [23] Goodchild, M.F., 2004. The validity and usefulness of laws in geographic information science and geography. Ann. Assoc. Am. Geogr. 94 (2), 300–303.
• [24] Griffith, D.A., 2000. Eigenfunction properties and approximations of selected incidence matrices employed in spatial analyses. Linear Algebra Appl. 321 (1-3), 95112.
• [25] Griffith, D.A., 2003. Spatial Autocorrelation and Spatial Filtering: Gaining Understanding through Theory and Scientific Visualization. Springer, Berlin.
• [26] Griffith, D.A., 2004. Extreme eigenfunctions of adjacency matrices for planar graph employed in spatial analyses. Linear Algebra Appl. 388, 201–219.
• [27] Griffith, D.A., 2008. Spatial-filtering-based contributions to a critique of geographically weighted regression (GWR). Environ. Plann. A 40 (11), 2751–2769.
• [28] Griffith, D.A., 2015. Approximation of Gaussian spatial autoregressive models for massive regular square tessellation data. Int. J. Geogr. Inf. Sci. 29 (12), 2143–2173.
• [29] Griffith, D.A., Chun, Y., 2014. Spatial autocorrelation and spatial filtering. In: Fischer, M.M. and Nijkamp, N (eds). Handbook of Regional Science. Springer, Berlin, Heidelberg, pp. 1477–1507.
• [30] Harris, P., Fotheringham, A. S., Juggins, S., 2010a. Robust geographically weighted regression: a technique for quantifying spatial relationships between freshwater acidification critical loads and catchment attributes. Ann. Assoc. Am. Geogr. 100 (2), 286–306.
• [31] Harris, R., Singleton, A., Grose, D., Brunsdon, C., Longley, P., 2010B. Grid enabling geographically weighted regression: A case study of participation in higher education in England. Tran. GIS, 14(1), 43–61.
• [32] Heaton, M.J., Datta, A., Finley, A., Furrer, R., Guhaniyogi, R., Gerber, F., Gramacy, R.B., Hammerling, D., Katzfuss, M., Lindgren F., Nychka, D.W., Sun, F., ZammitMangion, A., 2017. A case study competition among methods for analyzing larg spatial data. Arxiv, 1710.05013.
• [33] Helbich, M., Griffith, D.A., 2016. Spatially varying coefficient models in real estate: Eigenvector spatial filtering and alternative approaches. Comput. Environ. Urban Syst. 57, 1–11.
• [34] Kelejian, H.H., Prucha, I.R., 1998. A generalized spatial two-stage least squares procedure for estimating a spatial autoregressive model with autoregressive disturbances. J. Real Estate Finance Econ. 17 (1), 99–121.
• [35] LeSage, J.P., Pace, R.K., 2007. A matrix exponential spatial specification. J. Econom. 140 (1), 190–214.
• [36] LeSage, J.P., Pace, R.K., 2009. Introduction to Spatial Econometrics. Chapman and Hall/CRC, New York.
• [37] Li, M., Bi, W., Kwok, J.T., Lu, B.L., 2015. Large-scale Nyström kernel matrix approximation using randomized SVD. IEEE Trans Neural Netw. Learn. 26 (1), 152164.
• [38] Lindgren, F., Rue, H., Lindström, J., 2011. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. J. R. Stat. Soc. Series B Stat. Methodol. 73 (4), 423–498.
• [39] Lu, B., Harris, P., Charlton, M., Brunsdon, C., 2014. The GWmodel R package further topics for exploring spatial heterogeneity using geographically weighted models. Geo-spatial Inform. Sci. 17 (2), 85–101.
• [40] Lu, B., Brunsdon, C., Charlton, M., Harris, P., 2017. Geographically weighted regression with parameter-specific distance metrics. Int. J. Geogr. Inf. Sci. 31 (5), 982–998.
• [41] 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. Comput. Environ. Urban Syst, DOI: 10.1016/j.compenvurbsys.2018.03.012.
• [42] Moran, P. A., 1950. Notes on continuous stochastic phenomena. Biometrika 37 (1/2), 17–23.
• [43] Murakami, D., 2018. spmoran: An R package for Moran's eigenvector-based spatial regression analysis. Arxiv, 1703.04467.
• [44] Murakami, D., Griffith, D.A., 2015. Random effects specifications in eigenvector spatial filtering: a simulation study. J. Geog. Sci. 17 (4), 311–331.
• [45] Murakami, D., Lu, B., Harris, P., Brunsdon, C., Charlton, M., Nakaya, T., Griffith, D.A., 2017. The importance of scale in spatially varying coefficient modeling. Arxiv, 1709.08764
• [46] Murakami, D., Griffith, D.A., 2018. Eigenvector spatial filtering for large data sets: fixed and random effects approaches. Geogr. Anal. DIO: 10.1111/gean.12156.
• [47] Murakami, D., Yoshida, T., Seya, H., Griffith, D.A., Yamagata, Y., 2017. A Moran coefficient-based mixed effects approach to investigate spatially varying relationships. Spat. Stat. 19, 68–89.
• [48] Nakaya, T., Fotheringham, A.S., Charlton, M., Brunsdon, C., 2005. Geographically weighted Poisson regression for disease associative mapping. Stat. Med. 24 (16), 2695–2717.
• [49] Nychka, D., Bandyopadhyay, S., Hammerling, D., Lindgren, F., Sain, S., 2015. A multiresolution Gaussian process model for the analysis of large spatial datasets. J. Comput. Graph. Stat. 24 (2), 579–599.
• [50] Oshan, T.M., Fotheringham, A.S., 2018. A comparison of spatially varying regression coefficient estimates using geographically weighted and spatial‐filter‐based techniques. Geogr. Anal. 50 (1), 53–75.
• [51] Sang, H., Huang, J. Z., 2012. A full scale approximation of covariance functions for large spatial data sets. J. R. Stat. Soc. Series B Stat. Methodol. 74 (1), 111–132.
• [52] Silvester, J.R., 2000. Determinants of block matrices. Math. Gazette 84 (501), 460467
• [53] Smirnov, O., Anselin, L., 2001. Fast maximum likelihood estimation of very large spatial autoregressive models: a characteristic polynomial approach. Comput. Stat. Data Anal. 35 (3), 301–319.
• [54] Stein, M.L., 2014. Limitations on low rank approximations for covariance matrices of spatial data. Spat. Stat. 8, 1–19.
• [55] Tiefelsdorf, M., Griffith, D.A., 2007. Semiparametric filtering of spatial autocorrelation: the eigenvector approach. Environ. Plann. A 39 (5), 1193.
• [56] Tobler, W.R., 1970. A computer movie simulating urban growth in the Detroit region. Econ. Geogr. 46, 234–240.
• [57] Tran, H.T., Nguyen, H.T., Tran, V.T., 2016. Large-scale geographically weighted regression on Spark. Proceedings of the 2016 International Conference on Knowledge and Systems Engineering (KSE), 127–132.
• [58] Wheeler, D.C., Calder, C.A., 2007. An assessment of coefficient accuracy in linear regression models with spatially varying coefficients. J. Geog. Sci. 9 (2), 145–166.
• [59] Wheeler, D.C., Tiefelsdorf, M., 2005. Multicollinearity and correlation among local regression coefficients in geographically weighted regression. J. Geog. Sci. 7 (2), 161–187.
• [60] Wheeler, D.C., Waller, L., 2009. Comparing spatially varying coefficient models: case study examining violent crime rates and their relationships to alcohol outlets and illegal drug arrests. J. Geog. Sci. 11 (1), 1–22.
• [61] Yang, W., 2014. An extension of geographically weighted regression with flexible bandwidths. PhD Thesis. University of St Andrews.
• [62] Zhang, K., Kwok, J.T., 2010. Clustered Nyström method for large scale manifold learning and dimension reduction. IEEE Trans. Neural Netw. 21 (10), 1576–1587.