空间变系数过程模型
【摘 要】许多应用的目标是建立一套回归模型，以便在空间相关性假设下对感兴趣区域上的响应变量作出解释。在几乎所有这些工作中，回归系数都假定为在该区域内恒定。但在某些应用中，预测系数会在局部或子区域水平上有所不同，而这种情形正是本文的重点。尽管空间表面（ Surface ）的参数化建模是可能的（如多项式表面建模、样条建模等），但我们认为将其视为空间随机过程的一次实现更为自然和灵活。在本文中，我们展示了在高斯响应背景下，如何对这种建模方法进行形式化，使其能够在随机效应和残差分析方面提供更有吸引力的解释。我们还提供了广义线性模型和时空场景的扩展。文中将在单户住宅售价数据集上展示静态和动态建模和解释能力。
【参 考】
- Gelfand, A. E., Kim, H.-J., Sirmans, C. F., & Banerjee, S. (2003). Spatial modeling with spatially varying coefficient processes. Journal of the American Statistical Association, 98(462), 387–396. https://doi.org/10.1198/016214503000170 ）
- Finley, a. O., & Banerjee, S. (2020). Bayesian spatially varying coefficient models in the spBayes R package. Environmental Modelling & Software, 125, 104608. https://doi.org/10.1016/j.envsoft.2019.104608
1 引言
快速、廉价计算的广泛可用性以及功能强大、用户友好的地理信息系统 (GIS) 软件导致在房地产/金融、流行病学、环境计量学/生态、通讯等领域产生了大量空间和时空数据。反过来，这也推动统计社区增加了大量空间建模和数据分析的活动。
（1）问题的提出
在许多应用中，目标是建立回归模型来解释在感兴趣区域（例如 $D$）上观测到的响应变量，并且假设响应是空间相关的。也就是说，虽然可以通过均值完成一些空间建模，但我们还是可以预期响应之间存在依赖性，并且随着两个响应变量之间在空间中变得更接近，这种依赖会变得更强。对于以点为参考的连续响应场景，如果作出高斯分布假设（可能在一个变换后的尺度上），则通常直接使用高斯过程对这种依赖性进行建模。已经有大量这方面的文献资料， Cressie (1993) 的书也许是一个起点。例如，在二值响应或计数响应情况下，通常会采用分层建模方法，在第一阶段（数据模型阶段）使用指数族模型，然后将高斯分布的空间随机效应引入到（与链接函数相关的）变换后尺度的均值结构中（这里指采用广义线性回归模型对数据建模，而过程模型则是高斯响应的；参见如 Diggle、Tawn 和 Moyeed 1998）。
在几乎所有这些工作中，回归系数被假定在整个区域内是恒定的。在某些应用中，这并不合适。回归系数有可能会在局部或次区域水平上有所不同。例如：
- Assunçao、Gamerman (1999) 引入了
贝叶斯空间变化参数模型来检验巴西农业中的微观区域要素生产率和要素可替代性。 - Agarwal、Gelfand、Sirmans 和 Thibadeau（2003 年）在局部平稳空间建模背景下，引入了影响学区（和子学区）层面房价因素的局部回归模型。
- Gamerman、Moreira 和 Rue (2001) 通过模拟研究开发了一种用于空间变化回归模型的灵活建模方法。
这些作者都做出了相当严格的假设，即系数在特定空间单元（面元）上保持恒定。这些空间单元上的表面（Surface）使用独立或有条件的自回归规范进行建模。因此，会产生新的担忧，如空间单元分辨尺度的随意性问题、表面缺乏平滑度的问题、表面是否所有位置都能内插的问题等。
当使用点参考数据时，如果能够允许系数随位置变化，形成一个特定系数的空间表面似乎更具吸引力。例如，在我们的应用中，模拟了单户住宅的（对数）售价，解释变量包括房屋年龄、居住面积、其他区域面积、浴室数量。如果感兴趣的地区是一个城市或更大的大都市区，那么很明显，资本化率（如年龄）将在整个地区有所不同。该地区中，某些区域的二手房价值将高（或低）于其他地区。如果允许年龄的回归系数随位置变化，我们就可以解决上述问题。考虑到实际价值（如房地产评估），我们可以预测任意属性的回顾系数，而不仅仅是在调查期间出售的那些。在对特定污染物的环境暴露进行建模时也会出现类似的问题，其中协变量可能包括温度和降水。
（2）一些解决办法及存在问题
一种可能的方法是对参数化系数的空间表面建模。在最简单的情况下，可以采用任意多项式表面函数来指定；这可能会导致一系列表面过于有限或过于灵活。
在二维空间上使用样条曲面也可以引入更多灵活性（参见例如，Luo 和 Wahba 1998），但需要在空间中选择样条函数并确定结点数量和位置。此外，对于多个系数，样条曲面的多变量指定是必需的。
（3）本文思路
本文采用的方法更自然，至少同样灵活。我们将空间变系数的表面建模为空间过程的实现。对于多个回归系数，则使用多元空间过程模型。事实上，我们使用了一个平稳性假设，可以通过选择协方差函数（例如 Matèrn 类）来模拟所需过程进而实现平滑度。 Kent (1989) 和 Stein (1999a) 提供了关于单变量过程的讨论； Banerjee 和 Gelfand (2003) 提供了有关多元过程的讨论。
我们的建模框架采用贝叶斯方法。这在建议的设置中很有吸引力，因为我们对随机空间效应的推断特别感兴趣。特别是，我们可以在观测到和未观测到的位置获得 空间系数过程 的完整后验，以及所有模型参数的后验。在任何其他框架中似乎都无法对既没有观测值也没有残差解释的过程进行插值。对于高斯型响应，虽然可以通过似然方法进行一些推断，但它更有限，特别是区间估计依赖于可能不合适的渐近。对于非高斯数据，几乎肯定需要近似值，可以采用我们在 第 7 节 中提供的形式，但推断仍然会受到限制。
为了清晰地解释和实现：
- 首先，我们在
单个协变量的情况下开发一般性方法，因此有两个空间变化的系数过程，一个用于 “截距”，一个用于 “斜率”。 - 然后，我们转向多个协变量的情况。因为即使在基本多元回归设置中，系数估计通常也会显示出一些强相关性，因此预计空间变系数过程的集合是相互依赖的。因此，我们使用多变量过程模型。实际上，我们提出了进一步的概括，以建立多元回归模型的空间模拟（例如，参见 Goldstein 1995）。
- 我们还考虑了灵活的时空可能性。
上文提到的房地产场景提供了测点级别的协变量，其系数具有相当大的实际意义，来自洛杉矶巴吞鲁日的单户住宅销售数据集可以进行说明。除了表现出特殊地形的区域外，我们预计空间变化系数模型将比趋势面模型更有用。也就是说，将纬度和经度的多项式合并到均值结构中这种做法，不会被期望用于允许跨区域变化的代理模型。
文章的结构安排如下：
- 第 2 节详细介绍了单个协变量的高斯建模方法。
- 第 3 节讨论了多协变量的情况。
- 第 4 节提出了时空环境中一系列不同的系数指定。
- 第 5 节评论简要介绍了模型比较。
- 第 6 节展示了一个使用上述单户住宅销售数据的示例。
- 第 7 节以对广义线性模型设置的一些讨论结束。
- 第 8 节讨论大 n 问题j
- 第 9 节介绍可用的软件包。
2 单自变量的空间变系数模型
2.1 传统非变系数的空间过程模型
回想一下普通的单变量高斯平稳空间过程模型（ Cressie, 1993)：
$$
Y(\mathbf{s})=\mu(\mathbf{s})+W(\mathbf{s})+\epsilon(\mathbf{s}) \tag{1}
$$
其中：
$\mu(\mathbf{s})=\mathbf{x}(\mathbf{s})^T \boldsymbol{\beta}$ ，含截距项 $\beta_0$
$\epsilon(\mathbf{s})$ 是白噪声过程，其均值为 $\mathbb{E}(\epsilon(\mathbf{s}))=0$，方差为 $\operatorname{Var}(\epsilon(\mathbf{s}))=\tau^2$，且任意两点之间的噪声相互独立，即有 $\operatorname{Cov}( \epsilon(\mathbf{s})，\epsilon(\mathbf{s}^{\prime}))=0$
$W(\mathbf{s})$ 为空间随机效应项，是二阶平稳的零均值空间过程，该过程独立于噪声过程。
- 均值 $\mathbb{E}(W(\mathbf{s}))=0$
- 方差 $\operatorname{Var}(W(\mathbf{s}))=\sigma^2$
- 任意两点之间的协方差 $\operatorname{Cov} (W(\mathbf{s}), W(\mathbf{s}^{\prime}))=\sigma^2 \rho(\mathbf{s}, \mathbf { s}^{\prime}; \phi)$，即协方差可以由方差和相关函数 $\rho$ 计算得到，其中相关函数 $\rho$ 通常是两点之间某种空间距离的函数。
式 (1) 隐含地定义了一个分层模型。
2.2 截距的空间随机过程
如果令 $\mu(\mathbf{s})=\beta_0+\beta_1 x(\mathbf{s})$，并将 $W(\mathbf{s})$ 改写为 $W(\mathbf{s})=\beta_0(\mathbf{s})$ 的形式，则可以定义一个辅助符号 $\tilde{\boldsymbol{\beta}}_0(\mathbf{s})=\beta_0+\beta_0(\mathbf{s})$。也就是说，$\beta_0(\mathbf{s})$ 可以被解释为在位置 $\mathbf{s}$ 处对截距 $\beta_0$ 的随机空间调整，而对应的 $\tilde{\boldsymbol{\beta}}_0(\mathbf{s})$ 可以被视为关于截距的随机过程。
现在再来看 式（1），对于观测到的一组位置 $\mathbf{s}_1, \mathbf{s}_2, \ldots, \mathbf{s}n$，当给定 $\beta_0$， $\beta_1$，${\beta_0(\mathbf {s}i)}$、$\tau^2$ 时，响应变量 $Y(\mathbf{s}i)=\beta_0+\beta_1 x(\mathbf{s} {\mathbf{i}})+\beta_0(\mathbf{s}{\mathbf{i}})+\epsilon(\mathbf{s}{\mathbf{i} })$ 条件独立（ $i=1, \ldots, n$ ）。因此第一阶段的似然可以形式化为：
$$
L (\beta_0, \beta_1,{\beta_0 (\mathbf{s}_i)}, \tau^2; \mathbf{y})= (\tau^2)^{-\frac{n}{2}} \exp {-\frac{1}{2 \tau^2} \sum [Y(\mathbf{s}_i).. ..-(\beta_0+\beta_1 x(\mathbf{s}_i)+\beta_0(\mathbf{s}_i))^2]} \tag{2}
$$
根据前面的定义，向量 $\boldsymbol{\beta_0}=(\beta_0(\mathbf{s}_1), \ldots, \beta_0(\mathbf{s} _n))^T$ 显然服从零均值的多元高斯分布：
$$
f(\boldsymbol{\beta}_0 \mid \sigma_0^2, \phi_0)=N(\mathbf{0}, \sigma_0^2 H_0(\phi_0)) \tag{3}
$$
其中任意两个测点之间的相关性函数表示为 $(H_0(\phi_0))_{i j}=\rho_0(\mathbf{s}_i-\mathbf{s}_j ; \phi_0)$，很明显和两点之间的距离定义 $| \mathbf{s}-\mathbf{s^\prime} |$ 有关。对于接下来的所有讨论和示例，我们都将采用 Matèrn 相关函数 $\rho(\mathbf{h}, \phi)$ 来表示，其中 $\rho(\mathbf{h}, \phi) \propto(\gamma(|\mathbf{h}|))^v K_v (\gamma(|\mathbf{h}|))$。该相关函数中，$K_v$ 是一个修改后的 Bessel 函数，参数 $\boldsymbol{\phi}=(\gamma, \nu)$，其中 $\gamma$ 是衰减参数，$\nu$ 是平滑度参数。一旦指定了 $\beta_0$、$\beta_1$、$\tau^2$、$\sigma_0^2$ 和 $\phi_0$ 的先验，贝叶斯分层模型的定义就完成了。
在 式 (2) 和 式 (3) 下，可以对 $\boldsymbol{\beta_0}$ 进行积分，从而得到模型参数的边缘似然：
$$
\begin{aligned}
L (\beta_0, \beta_1, \tau^2, \sigma_0^2, \phi_0 ; \mathbf{y}) = &x | \sigma_0^2 H_0 (\phi_0) + \tau^2 I |^{-\frac{1}{2}} \
&\times \exp { -\frac{1}{2} (\mathbf{y}-\beta_0 \mathbf{1}-\beta_1 \mathbf{x})^T (\sigma_0^2 H_0(\phi_0)+\tau^2 I)^{-1} \times(\mathbf{y}-\beta_0 \mathbf{1}-\beta_1 \mathbf{x})} \tag{4}
\end{aligned}
$$
其中 $\mathbf{x}=(x(\mathbf{s}_1), \ldots, x(\mathbf{s}_n))^T$。
我们可以关注下其与普通高斯随机效应模型的对比：
普通模型为 $Y_{i j}=\beta_0+\beta_1 x_{i j}+\alpha_i+e_{i j}$，其中 $\alpha_i$ 为 $N(0, \sigma_\alpha^2)$ 的独立同分布样本，$\epsilon_{i j}$ 是 $N(0, \sigma_\epsilon^2)$ 的独立同分布样本。此模型通常无法识别（分离）方差分量 $\sigma_\alpha^2$ 和 $\sigma_\epsilon^2$。但在新模型中，由于 $\beta_0 (\mathbf{s}_i)$ 之间存在依赖关系，而噪声分量之间相互独立，因此能够识别不同的方差分量，如 式 (4) 所示。此外，如果 $U(\mathbf{s})$ 表示回归模型中的总误差，则 $U(\mathbf{s})$ 可以被分解为 “截距过程误差” 和 “纯噪声误差”。
如果我们采用随机模拟的方法用 式（4） 中的边缘似然来拟合贝叶斯模型，则必须获取来自于后验分布 $f(\beta_0, \beta_1, \tau^2, \sigma_0^2, \phi \mid \mathbf{y})$ 的样本。但布 $\boldsymbol{\beta}_0 \mid \mathbf{y}$ 的样本仍然可以一对一获得，因为：
$$
\begin{aligned}
f(\boldsymbol{\beta}_0 \mid \mathbf{y})=\int f(\boldsymbol{\beta}_0 \mid \beta_0, \beta_1, \tau^2,.&.\sigma_0^2, \phi_0, \mathbf{y}) \
& \times f(\beta_0, \beta_1, \tau^2, \sigma_0^2, \phi_0 \mid \mathbf{y})
\end{aligned} \tag{5}
$$
其中
$$
\begin{aligned}
f(\beta_0 \mid \beta_0, \beta_1, \tau^2, \sigma_0^2, \phi_0, \mathbf{y}) \
&=\mathcal{N} ( ( \frac{1}{\tau^2} I + \frac{1}{\sigma_0^2} H_0^{-1} ( \phi_0 ))^{-1} \frac{1}{\tau^2} (\mathbf{y} - \beta_0 \mathbf{1}- \beta_1 \mathbf{x}),
(\frac{1}{\tau^2} I + \frac{1}{\sigma_0^2} H_0^{-1} ( \phi_0 ))^{-1} )
\end{aligned}
$$
我们也可以从 $\beta_0(\mathbf{s})$ 过程的后验分布中获取新位置处的样本，例如对于新位置 $\mathbf{s}_{\text {new}}$，我们可以得到 $\beta_0(\mathbf{s})$ 表面的如下插值预测：
$$
f(\beta_0(\mathbf{s}{\text {new }}) \mid \mathbf{y})=\int f(\beta_0(\mathbf{s}{\text {new }}) \mid \boldsymbol{\beta}_0, \sigma_0^2, \phi_0) f(\beta_0, \sigma_0^2, \phi_0 \mid \mathbf{y}) \tag{6}
$$
积分下的第一个密度是一个单变量高斯分布，可以直接从 $\beta_0(\mathbf{s})$ 的定义给出。对于给定 $\mathbf{y}$ 时的预测 $y(\mathbf{s}_{\text {new }})$，我们要求：
$$
\begin{aligned}
f(y(\mathbf{s}{\text {new }}) \mid \mathbf{y})=\int & f(y(\mathbf{s}{n e w}) \mid \beta_0, \beta_1, \beta_0(\mathbf{s}{n e w}), \tau^2) \
& \times f(\beta_0(\mathbf{s}{\text {new }}) \mid \beta_0, \sigma_0^2, \phi_0) \
& \times f(\boldsymbol{\beta}_0, \beta_0, \beta_1, \tau^2, \sigma_0^2, \phi_0 \mid \mathbf{y})
\end{aligned} \tag{7}
$$
积分符号下的第一项是正态密度。同样，从该预测分布中获取样本也很简单。
2.3 斜率的空间随机过程
根据上述改进，我们很快就会想到一个空间变系数模型的设计思想。假设我们令自变量的系数 $\beta_1$ 也出现随机效应：
$$
Y(\mathbf{s})=\beta_0+\beta_1 x(\mathbf{s})+\beta_1(\mathbf{s}) x(\mathbf{s})+\epsilon(\mathbf{s}) \tag{8}
$$
与传统模型的推导很类似，只是在 式（8）中，$\beta_1(\mathbf{s})$ 是一个二阶平稳的零均值高斯过程，其方差为 $\sigma_1^2$，相关函数为 $\rho_1(\cdot ; \phi_1 )$。如果令 $\tilde{\boldsymbol{\beta}}_1(\mathbf{s})=\beta_1+\beta_1(\mathbf{s})$，则现在 $\beta_1(\mathbf{s})$ 可以被解释为在位置 $\mathbf{s}$ 处，对整体斜率 $\beta_1$ 的随机空间调整。等效地，$\tilde{\boldsymbol{\beta}}_1(\mathbf{s})$ 可以被视为关于斜率的一个随机过程。实际上，我们正在使用一个无限维函数来解释 $x(\mathbf{s})$ 和 $Y(\mathbf{s})$ 之间的关系。
式 (8) 对 式 (2) 和 式 (3) 做了明显修改。特别是，由此产生的边缘似然变成了：
$$
\begin{aligned}
&L(\beta_0, \beta_1, \tau^2, \sigma_1^2, \phi_1 ; \mathbf{y}) \
&=|\sigma_1^2 D_x H_1(\phi_1) D_x+\tau^2 I|^{-\frac{1}{2}} \
&\quad \times \exp {-\frac{1}{2}(\mathbf{y}-\beta_0 \mathbf{1}-\beta_1 \mathbf{x})^T(\sigma_1^2 D_x H_1(\phi_1) D_x+\tau^2 I)^{-1}. \
&.\quad \times(\mathbf{y}-\beta_0 \mathbf{1}-\beta_1 \mathbf{x})}
\end{aligned} \tag{9}
$$
其中 $D_x$ 是以 $(D_x)_{i i}=x(\mathbf{s}_i)$ 为对角线的对角矩阵。此外，由于 $\boldsymbol{\beta}_1=(\beta_1(\mathbf{s}_1), \ldots, \beta_1(\mathbf{s}_n))^T$，使用 类似 式 (5) 和 式 (6) 的方式，我们可以得到 $f(\boldsymbol{\beta}1 \mid \mathbf{y})$ 和新位置处 $f(\beta_1(\mathbf{s}{\text {new }}) \mid \mathbf{y})$ 的样本 。
请注意，式 (8) 为数据提供了一个异构的、非平稳的过程，而与过程 $\beta_1(\mathbf{s})$ 的协方差函数选择无关。这里，$\operatorname{Var}(Y(\mathbf{s}) \mid \beta_0, \beta_1, \tau^2, \sigma_1^2, \phi_1)= x^2(\mathbf{s}) \sigma_1^2+\tau^2$， $\operatorname{Cov}(Y(\mathbf{s}), Y(\mathbf{s}^{\prime} ) \mid \beta_0,\beta_1, \tau^2, \sigma_1^2, \phi_1)=$ $\sigma_1^2 x(\mathbf{s}) x(\mathbf{s} ^{\prime}) \rho_1(\mathbf{s}-\mathbf{s}^{\prime} ; \phi_1)$。
我们在实际工作中观察到，只有当 $x(\mathbf{s})>0$ 时，式 (8) 才是合理的。事实上，在其他模型中倡导的中心化和缩放在这里并不合适。使用中心化后的 $x(\mathbf{s})$ ，我们会发现一些站不住脚的现象，方差 $\operatorname{Var}(Y(\mathbf{s}))$ 减少了，而 $x(\mathbf {s})$ 却增加了。更糟糕的是，对于本质上是中心化的 $x(\mathbf{s})$，我们会发现 $Y(\mathbf{s})$ 独立于任意 $Y(\mathbf{s}^{\prime} )$ 。此外，对 $x(\mathbf{s})$ 的缩放也基本无济于事。$\beta_1(\mathbf{s})$ 将被反向重新缩放，因为模型仅识别 $\beta_1(\mathbf{s}) x(\mathbf{s})$。
这导致了对向量 $\mathbf{x}$ 与向量 $\mathbf{1}$ 可能近似共线性的担忧。式 (9) 表明，如果 $\mathbf{x} \approx c \mathbf{1}$，将出现结果较差的似然。幸运的是，我们可以将 式 (8) 重参数化为 $Y(\mathbf{s})=$ $\beta_0^{\prime}+\beta_1^{\prime} \tilde{x}(\mathbf{s})+\beta_1(\mathbf{s}) x(\mathbf{s})+\epsilon(\mathbf{s})$，其中 $\tilde{x}(\mathbf{s})$ 根据 $\beta_0^{\prime}$ 和 $\beta_1^{\prime}$ 的定义中心化和缩放。现在 $\tilde{\boldsymbol{\beta}}_1(\mathbf{s})=$ $\beta_1^{\prime} / s_x+\beta_1(\mathbf{s})$，其中 $s_x$ 是 $x(\mathbf{s})$ 的样本标准差。
就像 式 (4) ，我们可以与普通的纵向数据线性模型进行类比，其中 $Y_{i j}=\beta_0+$ $\beta_1 x_{i j}+\boldsymbol{\beta}{1 i} x{i j}+\epsilon_{ i j}$，即每个样本都对应有随机的斜率。此外， 式 (8) 中回归模型的总误差 $U(\mathbf{s})$ 现在被划分为 “斜率过程误差” 和 “纯噪声误差”。
2.4 单自变量的通用模型
包含 (1) 和 (7) 的更通用的模型定义将是：
$$
Y(\mathbf{s})=\beta_0+\beta_1 x(\mathbf{s})+\beta_0(\mathbf{s})+\beta_1(\mathbf{s}) x(\mathbf{s})+\epsilon(\mathbf{s}) \tag{10}
$$
式 (10) 引入了平行于普通线性模型的 截距过程 和 斜率过程 。新模型需要一个双变量空间过程来定义来确定 $\boldsymbol{\beta_0}$ 和 $\boldsymbol{\beta_1}$ 的联合分布。我们将在 第 3 节 再次探讨此问题。但是，在过程独立性假设下，式(10) 可以很容易地在 $\boldsymbol{\beta_0}$ 和 $\boldsymbol{\beta_1}$ 上被边缘化，生成：
$$
\begin{aligned}
L(\beta_0, \beta_1, \tau^2, \sigma_0^2, \sigma_1^2, \phi_0, \phi_1 ; \mathbf{y}) \
=&|\sigma_0^2 H_0(\phi_0)+\sigma_1^2 D_x H_1(\phi_1) D_x+\tau^2 I|^{-\frac{1}{2}} \
& \times \exp {-\frac{1}{2}(\mathbf{y}-\beta_0 \mathbf{1}-\beta_1 \mathbf{x})^T.\
& \times(\sigma_0^2 H_0(\phi_0)+\sigma_1^2 D_x H_1(\phi_1) D_x+\tau^2 I)^{-1} \
&.\times(\mathbf{y}-\beta_0 \mathbf{1}-\beta_1 \mathbf{x})}
\end{aligned} \tag{11}
$$
同样，通过模拟方法抽取来自 $f(\boldsymbol{\beta}0 \mid \mathbf{y})$ 和 $f(\boldsymbol{\beta}1 \mid \mathbf{y})$ 的样本，预测 $f(\boldsymbol{\beta}0(\mathbf {s}{\text {new }}) \mid \mathbf{y})$、 $f(\beta_1(\mathbf{s}{\text {new }}) \mid \mathbf{y})$, 和 $f(y(\mathbf{s}{\text {new }}) \mid \mathbf{y})$ 也很直接。
在任意位置预测空间表面的能力，为贝叶斯推断方法提供了一个令人信服的案例。经典克里金法无法解决这个问题，因为 $\beta_0(\mathbf{s}_i)$ 或 $\beta_1(\mathbf{s}_i)$ 无法被观测到。此外，在 式（10） 中，总误差被划分为三个独立组分，且具有明确的解释。在任意 $\mathbf{s}_i$ 处，误差分量的相对大小可以通过三个部分来研究：
- 一是 $E(\beta_0(\mathbf{s}_i) \mid \mathbf{y})$
- 二是 $E(\beta_1(\mathbf{s}_i) x(\mathbf{s}_i) \mid \mathbf{y})$
- 三是 $E(\epsilon(\mathbf{s}_i) \mid \mathbf{y})$。
三者之间的相对贡献，可以通过计算期望值 $E(\tau^2 \mid \mathbf{y})$, $E(\sigma_0^2 \mid \mathbf{y})$ 和 $\bar{x}^2 E(\sigma_1^2 \mid \mathbf{y})$ 来衡量，其中 $\bar{x}=\sum x(\mathbf{s} _i) / n$。在 $\rho_0$ 和 $\rho_1$ 都是各向同性情况下，范围的后验推断结果可以被用来对截距过程和斜率过程的变程作比较。两个系数的后验均值曲面 $E(\beta_0(\mathbf{s}) \mid \mathbf{y})$ 和 $E(\beta_1(\mathbf{s}) \mid \mathbf{y })$ 可以在任意位置的网格处获得，并使用常规软件以图形方式显示。
使用 Agarwal 和 Gelfand (2001) 讨论的切片 Gibbs 采样很容易拟合使用 式 (4)、 式 (9) 或 式 (11) 中边缘似然的贝叶斯模型（也参见 Neal 2002）。
特别是，在 式 (11) 下，我们可以添加一个辅助变量 $U \sim \operatorname{unif}(0, L(\beta_0, \beta_1, \tau^2, \sigma_0^2, \sigma_1^2, \phi_0, \phi_1 ; \mathbf{y}))$，则后验 $f(\beta_0, \beta_1, \tau^2, \sigma_0^2, \sigma_1^2, \phi_0, \phi_1, U \mid \mathbf{y}) \propto I( U<L) \pi(\beta_0),\beta_1, \tau^2, \sigma_0^2, \sigma_1^2, \phi_0, \phi_1)$，其中 $\pi$ 是参数的先验。切片 Gibbs 采样使用均匀抽取来更新 $U$，并使用受指标限制的先验抽取来更新其他参数。该算法是 “现成的”，不需要调整。它比 Metropolis 替代方案收敛得更快，并且可以避免这些替代方案经常出现的自相关问题。
3 多自变量的空间变系数模型
本节我们转向在位置 $\mathbf{s}$ 处存在 $p$ 个协变量的更复杂场景。$\mathbf{X}(\mathbf{s})$ 代表点 $\mathbf{s}$ 处的 $p \times 1$ 协变量向量，为了向量化表示方便起见，我们在 $\mathbf{X}(\mathbf{s})$ 中再增加一个 $1$ 作为截距项。则 式 (10) 可以泛化为：
$$
Y(\mathbf{s})=\mathbf{X}^T(\mathbf{s}) \tilde{\boldsymbol{\beta}}(\mathbf{s})+\epsilon(\mathbf{s}) \tag{12}
$$
其中 $\tilde{\boldsymbol{\beta}}(\mathbf{s})$ 被假设为服从 $p$ 元变量的空间过程模型。
对于观测到的位置 $\mathbf{s}_1, \mathbf{s}_2, \ldots, \mathbf{s}_n$，令 $X^T$ 为 $n \times np$ 的块对角矩阵，其中块为第 $i$ 行的 $\mathbf{X}^T (\mathbf{s}_i)$。我们可以写成 $\mathbf{Y}=X^T \tilde{\boldsymbol{\beta}}+\epsilon$，其中 $\tilde{\boldsymbol{\beta}}$ 为 $n p \times 1$ 的向量，是 $\tilde {\boldsymbol{\beta}}(\mathbf{s})$ 串接后形成的向量，而 $\epsilon \sim \mathcal{N}(0, \tau^2 I)$ 。
在实际工作中，假设 $\tilde{\boldsymbol{\beta}}(\mathbf{s})$ 中的组分过程之间相互独立可能是不合适的。典型的例子是，在简单线性回归的简单情况下，经常会看到斜率和截距之间的负相关而非独立。在 第 6 节 的示例表明，在结合依赖性时模型性能得到显著改进。为了给 $\tilde{\boldsymbol{\beta}}(\mathbf{s})$ 指定多元高斯过程，我们需要均值和交叉协方差函数。对于前者，依据 第 2 节 ，我们将其设为 $\mu_{\boldsymbol{\beta}}=(\beta_1, \ldots, \boldsymbol{\beta}_p)^T$。对于后者，我们需要一个有效选择 $p$ 个变量。在后续工作中，我们将使用 Mardia 和 Goodall (1993) 的可分离设计（在这方面另见 Banerjee 和 Gelfand 2002）。
更准确地说，令 $\boldsymbol{C}(\mathbf{s}, \mathbf{s}^{\prime})$ 为 $p\times p$ 矩阵，其中第 $(l, m)$ 个元素为 $\operatorname{Cov}(\tilde{\boldsymbol{\beta}}_l(\mathbf{s}), \tilde{\boldsymbol{\beta}_m}(\mathbf{s}^{ \prime}))$ ，并令
$$
C(\mathbf{s}, \mathbf{s}^{\prime}){l m}=\rho(\mathbf{s}-\mathbf{s}^{\prime} ; \phi) \tau{l m} \tag{13}
$$
其中 $\rho$ 是一个有效的二维标量值的相关函数，$T_{p \times p}$ 使得 $(T){l m}=\tau{l m}$ 是正定对称的。换句话说，$T$ 是与任何空间位置的观测向量相关联的协方差矩阵，$\rho$ 捕获了跨空间关联的衰减。和上一节一样，如果我们将 $\rho(\mathbf{s}-\mathbf{s}^{\prime}; \phi)$ 的集合收集到一个 $n \times n$ 的矩阵 $H(\phi) $ 中，那么，当 $\otimes$ 表示 Kronecker 积时，$\boldsymbol{\beta}$ 的分布是
$$
\tilde{\boldsymbol{\beta}} \sim \mathcal{N}(\mathbf{1}{n \times 1} \otimes \boldsymbol{\mu}{\boldsymbol{\beta}}, H(\phi) \otimes T) \tag{14}
$$
与上一节一样，如果 $\tilde{\boldsymbol{\beta}} = \boldsymbol{\beta} + \mathbf{1}{n \times 1} \otimes \boldsymbol{\mu}{\boldsymbol{\beta}}$，那么我们可以将 式 (12) 写为
$$
Y(\mathbf{s})=\mathbf{X}^T(\mathbf{s}) \boldsymbol{\mu}_{\boldsymbol{\beta}}+\mathbf{X}^T(\mathbf{s}) \boldsymbol{\beta}(\mathbf{s})+\epsilon(\mathbf{s}) \tag{15}
$$
在式 (15) 中，回归模型中的总误差被划分为 $p+1$ 个组分，每个组分都有一个明显的解释。根据 第 2 节 ，使用 式 (12) 和 式 (14)，我们可以对 $\boldsymbol{\beta}$ 进行积分以获得
$$
\begin{aligned}
&L(\boldsymbol{\mu}{\boldsymbol{\beta}}, \tau^2, T, \phi ; \mathbf{y}) \
&=|X(H(\phi) \otimes T) X^T+\tau^2 I|^{-\frac{1}{2}} \
&\quad \times \exp {-\frac{1}{2}(\mathbf{y}-X(\mathbf{1} \otimes \mu{\boldsymbol{\beta}}))^T(X(H(\phi) \otimes T) X^T+\tau^2 I)^{-1} \
&\qquad \times(\mathbf{y}-X(\mathbf{1} \otimes \boldsymbol{\mu}_{\boldsymbol{\beta}}))} \tag{16}
\end{aligned}
$$
可能令人生畏的 式 (16) 仍然只涉及 $n \times n$ 的矩阵。
完成贝叶斯模型，需要首先指定 $f(\mu_{\boldsymbol{\beta}}, \tau^2, T, \phi)$ 的先验，我们假设其为乘积的形式 $f(\boldsymbol{\mu}_{\boldsymbol{\beta}}) f(\tau^2) f(T) f(\phi)$。这些分量分别服从高斯、逆gamma、逆 Wishart-gamma、gamma 分布，当使用 Matèrn 相关函数时，$\phi=(\gamma,\nu)$ 。使用 式 (16) 和上述先验，模型很容易通过切片 Gibbs 采样进行基于模拟的拟合。
关于类似于 式 (5) 的预测，$f(\tilde{\boldsymbol{\beta}} \mid \mathbf{y})$ 可以利用来自 $f( \boldsymbol{\mu}{\boldsymbol{\beta}}, \tau^2, T, \phi \mid \mathbf{y})$ 的样本使用 $f(\tilde{\boldsymbol{\beta}} \mid \boldsymbol{\mu} {\boldsymbol{\beta}}, \tau^2, T, \phi, y)$ 一对一获得。即 $N(B \mathbf{b}, B)$，其中 $B=(X^T X / \tau^2 + H^{-1}(\phi) \otimes T^{-1})^{-1}$ 和 $\mathbf{b}=X^T \mathbf{y} / \tau^2 + (H^{-1}(\phi) \otimes T^{-1})(\mathbf{1} \otimes \boldsymbol{\mu}{\boldsymbol{\beta}})$ 。 $B$ 是一个 $np \times np$ 的大型矩阵，但对于 $\tilde{\boldsymbol{\beta}}$ 的采样，只需要 $B$ 的 Cholesky 分解并且仅用于保留的后验样本。新位置 $\mathbf{s}{\text {new }}$ 处的预测，类似于 式 (6)，需要 $f(\tilde{\boldsymbol{\beta}}(\mathbf{s} {\text {new }}) | \tilde{\boldsymbol{\beta}}, \boldsymbol{\mu}{\boldsymbol{\beta}}, \tau^2, T, \phi)$。 将 $\mathbf{h}{\text {new }}(\phi)$ 定义为 $n \times 1$ 向量，其中第 $i$ 行元素 $\rho(\mathbf{s}i- \mathbf{s}{\text {new }} ; \phi)$，这个分布是均值为 $\boldsymbol{\mu}{\boldsymbol{\beta}}+(\mathbf{h}{\text {新 }}^T(\phi) \otimes T)(H^{-1}(\phi) \otimes T^{-1})(\tilde{\boldsymbol{\beta}}-\ mathbf{1}{n x 1} \otimes \boldsymbol{\mu}{\boldsymbol{\beta}})=$ $\boldsymbol{\mu}{\boldsymbol{\beta}}+(\mathbf{h}{\text {new } }^T(\phi) H^{-1}(\phi) \otimes I)(\tilde{\boldsymbol{\beta}}-\mathbf{1}{n x 1} \otimes \boldsymbol{\mu }{\boldsymbol{\beta}})$ 和协方差矩阵为 $T-(\mathbf{h}{\text {new }}^T(\phi) \otimes T)(H^{-1 }(\phi) \otimes T^{-1})(\mathbf{h}{\text {new }}(\phi) \otimes T)=$ $(I- \mathbf{h}{\text {new }}^T(\phi) H^{-1}(\phi) \mathbf{h}{\text {new }}(\phi)) T$ 的高斯分布。最后，$Y(\mathbf{s}{\text {new }})$ 的预测分布 $f(Y(\mathbf{s}_{\text {new }}\ right) \mid \mathbf{y})$ 可以用类似于 式 (7) 的方法采样。
当我们在位置 $\mathbf{s}$ 处重复测量时，注意 式 (12) 的扩展。也就是说，假设我们有
$$
Y(\mathbf{s}, l)=\mathbf{X}^T(\mathbf{s}, l) \boldsymbol{\beta}(\mathbf{s})+\epsilon(\mathbf{s}, l) \tag{17}
$$
其中 $l=1, \ldots, L_{\mathbf{s}}$ ，$L_{\mathbf{s}}$ 是在 $\mathbf{s}$ 处测量的次数， $\epsilon(\mathbf{ s}, l)$ 仍然是白噪声。例如，在我们之前提到的房地产估值场景中，$\mathbf{s}$ 可能表示一个公寓楼的位置，$l$ 索引该楼中已售出的公寓，其中第 $l$ 套公寓具有 $\mathbf{X}(\mathbf{s}, l)$。进一步假设 $\mathbf{Z}(\mathbf{s})$ 为表示测点处各种等级特征的 $r \times 1$ 向量。对于一个公寓楼，这些特征可能包括 “提供的便利设施” 、“与中央商务区的距离” 等。然后 式 (17) 可以扩展到 Goldstein (1995) 或 Raudenbush 和 Bryk (2002) 中所定义的 多级模型（Multilevel Model）。特别是，我们可以写
$$
\boldsymbol{\beta}(\mathbf{s})=(\begin{array}{c}
\mathbf{Z}^T(\mathbf{s}) \boldsymbol{\gamma}_1 \
\vdots \
\mathbf{Z}^T(\mathbf{s}) \boldsymbol{\gamma}_p
\end{array})+\mathbf{W}(\mathbf{s}) \tag{18}
$$
在 式 (18) 中，$\boldsymbol{\gamma}_j$ 是与 $\tilde{\boldsymbol{\beta}}_j(\mathbf{s})$ 和 $\mathbf{W}(\mathbf{s})$ 相关的 $r\times 1$ 向量，服从均值为 $\mathbf{0}$ 的多元高斯空间过程。在 式 (18) 中，如果 $\mathbf{W}(\mathbf{s})$ 是独立的，那么我们将有一个普通的多级模型。在 $\mathbf{Z}(\mathbf{s})$ 只是一个捕获截距的标量时，我们就回到了本节的初始模型。
4 时空数据的空间变系数模型
前面部分建模的自然扩展是数据在空间位置上随时间相关的情况。此类数据经常出现在生态、环境和气象环境中。如果假设时间在一定尺度上被离散化为一组有限的等距点，那么我们可以概念化一个时间序列的空间过程，这些空间过程仅在空间位置 $\mathbf{s}_1, \ldots, \mathbf {s}_n$ 处进行观测。
采用与 式 (10) 相似的通用符号，令
$$
Y(\mathbf{s}, t)=\mathbf{X}^T(\mathbf{s}, t) \tilde{\boldsymbol{\beta}}(\mathbf{s}, t)+\epsilon(\mathbf{s}, t), \quad t=1,2, \ldots, M \tag{19}
$$
也就是说，我们引入了时空变化的截距和斜率。如果我们写成 $\boldsymbol{\beta}(\mathbf{s}, t)=\boldsymbol{\beta}(\mathbf{s}, t)+\mu_{\boldsymbol{\beta}}$，那么总误差将被划分为 $p+1$ 个包含 $\epsilon(\mathbf{s}, t)$ 的空间-时间截距片段，每个片段都有明显的解释。我们继续假设 $\epsilon(\mathbf{s}, t)$ 独立同分布且服从 $\mathcal{N}(0, \tau^2)$，但需要为 $\tilde{\boldsymbol{\beta}}(\mathbf{s}, t)$ 指定一个模型。无论如何，式 (19) 定义了一个非平稳过程，该过程的数学期望为 $\mathbb{E}[Y(\mathbf{s}, t)] = \mathbf{X}^T(\mathbf{s}, t) \tilde{\boldsymbol{\beta } }(s, t)$，方差为 $\operatorname{Var}(Y(\mathbf{s}, t)) = \mathbf{X}^T(\mathbf{s}, t) \Sigma_{\tilde{\boldsymbol{\beta}}(\mathbf{s}, t)} \mathbf{X}(\mathbf{s},t) + \tau^2$, 而协方差为 $\operatorname{Cov}(Y(\mathbf{s } , t), Y(\mathbf{s}^{\prime}, t^{\prime})) = \mathbf{X}^T(\mathbf{s}, t) \Sigma_{\tilde{\boldsymbol{\beta}}(\mathbf{s}, t), \tilde{\boldsymbol{\beta}}(\mathbf{s}^{\prime}, t^{\prime})} \mathbf{X}(\mathbf{s}^{\prime}, t^{\prime})$。
我们为 $\boldsymbol{\beta}(\mathbf{s}, t)$ 提出了四个模型。
（1）模型一
当时间序列通常很短时，平行于惯用的纵向数据建模假设，我们可以设置
$$
\boldsymbol{\beta}(\mathbf{s}, t)=\boldsymbol{\beta}(\mathbf{s})
$$
其中 $\boldsymbol{\beta}(\mathbf{s})$ 的建模与前述小节相同。模型 1 可以被看作是一个局部线性模型。
（2） 模型二
接下来，我们有
$$
\boldsymbol{\beta}(\mathbf{s}, t)=\boldsymbol{\beta}(\mathbf{s})+\alpha(t)
$$
其中 $\boldsymbol{\beta}(\mathbf{s})$ 再次与 模型 1 中一样。在建模 $\boldsymbol{\alpha}(t)$ 时，我们考察了两种似然。第一种似然将 $\alpha_k(t)$ 视为时间虚拟变量，将这组 $pM$ 个变量视为先验独立且同分布的。第二种似然将 $\alpha(t)$ 建模为随机游走或自回归过程，我们可以假设组件在 $k$ 上是独立的，但为了更通用，我们将其视为依赖的，使用 式 (13) 的类比，将 $\mathbf{s}$ 替换为 $t$ 和 $\rho$ ，现在是一个有效的一维相关函数。例如，Gelfand、Ecker、Knight 和 Sirmans (2003) 在错误结构的通常建模情况下讨论了这种空间和时间的可加性。
（3）模型三
在 模型 3 中，我们考虑了 Waller、Carlin、Xia 和 Gelfand (1997) 中面元嵌套效应（另见 Gelfand 等，2002 年）。特别是有：
$$
\boldsymbol{\beta}(\mathbf{s}, t)=\boldsymbol{\beta}^{(t)}(\mathbf{s})
$$
在这里，我们有嵌套在时间中的空间变化系数过程。假设这些过程在$t$（本质上是时间虚拟过程）中是独立的，并允许系数过程的时间演化。在第 3 节之后，过程 $\boldsymbol{\beta}^{(t)}(\mathbf{s})$ 将是均值 0 ，在时间 $t$, $C^{(t )}(\mathbf{s}, \mathbf{s}^{\prime})$, 其中 $(C^{(t)}(\mathbf{s}, \mathbf{ s}^{\prime})){l m}=\rho(\mathbf{s}-\mathbf{s}^{\prime} ; \phi^{(t)} ) \tau{l m}^{(t)}$。我们已经指定模型 3 具有跨时间的通用 $\mu_{\boldsymbol{\beta}}$，这使得与我们提出的其他模型具有一定的可比性。但是，我们可以通过将 $\mu_{\boldsymbol{\beta}}$ 替换为 $\mu_{\boldsymbol{\beta}}^{(t)}$ 来增加灵活性。
（4）模型四
最后，模型 4 提出了空间和时间上的可分离协方差规范，扩展了 Gelfand、Zhu 和 Carlin (2001) 的工作：
$$
\boldsymbol{\beta}(\mathbf{s}, t)
$$
使得 $\Sigma_{[\boldsymbol{\beta}(\mathbf{s}, t), \boldsymbol{\beta}(\mathbf{s}^{\ prime}, t^{\prime})]}=\rho^{(1)}(\mathbf{s}-\mathbf{s}^{\prime} ; \phi) \times$ $\rho^{(2)}(t-t^{\prime} ; \gamma) T
$
其中$\rho^{(1)}$ 是有效的二维相关函数，$\rho^{(2)}$ 是有效的一维选择，$T$ 是正定对称的。这里 $\rho^{(1)}$ 获得了与前面部分一样的空间关联，它随着时间的推移被 $\rho^{(2)}$ 衰减。得到的全向量 $\boldsymbol{\beta}$ 的协方差矩阵，被站点和站点内的时间阻塞，具有方便的形式 $H_2(\gamma) \otimes H_1(\phi) \otimes T$。
在上述每个模型中，我们都可以像前面部分中所做的那样边缘化 $\boldsymbol{\beta}(\mathbf{s}, t)$。根据模型，按站点或按时间对数据分块在计算上可能更方便。我们省略了细节，只注意到对于 $n$ 个站点和 $T$ 个时间点，得到的似然将涉及 $nT \times nT$ 矩阵的行列式和逆矩阵。
请注意，所有上述建模都可以应用于横截面数据的情况，其中观测到的位置集随 $t$ 而变化。例如，房地产数据就是这种情况。我们仅在交易时观测售价。在 $t$ 年有 $n_t$ 个位置，除 模型 3 之外的所有似然都将涉及 $\sum n_t \times \sum n_t$ 矩阵。
5 模型比较
在纯空间情况或时空情况下，我们可能会寻求比较模型。例如，在 第 2 节 和 第 3 节 中的空间情况下，我们可以考虑只有系数的子集在空间上变化的模型，其中系数过程是独立的，或者假设多元依赖版本。在时空情况下，我们可以考虑模型 1-4（在模型 2 和 3 的情况下可能包含子模型）。
我们使用 Gelfand 和 Ghosh (1998) 的后验预测损失方法。在时空上下文中说明，对于 式 (19) 中的每个 $(\mathbf{s}, t)$，令 $Y_{\text {new }}(\mathbf{s}, t)$ 是在该位置获得的新观测值和给定 $\mathbf{X}(\mathbf{s}, t)$ 的时间。如果在特定模型下，$\mu(\mathbf{s}, t)$ 和 $\sigma^2(\mathbf{s}, t)$ 表示给定 $Y_{obs}$（所有观测到的 $Y$ 的集合）时 $Y_{\text {new }}(\mathbf{s}, t)$ 的预测分布的均值和方差，则准则变为
$$
D_k=\sum_{(\mathbf{s}, t)}(Y_{o b s}(\mathbf{s}, t)-\mu(\mathbf{s}, t))^2+k \sum_{(\mathbf{s}, t)} \sigma^2(\mathbf{s}, t) \tag{20}
$$
在 式 (20) 中，第一项是拟合优度分量 $(G)$，第二项是模型复杂度分量 $(P)$ 的惩罚。常数 $k$ 对这些组件进行加权，通常设置为 $1$ 。选择使 式 (20) 产生最小值的模型。
6 应用案例
（略）
7 广义线性模型设置
我们简要考虑 式（12） 的广义线性模型版本，将高斯第一阶段替换为：
$$
f(y(\mathbf{s}_i) \mid \theta(\mathbf{s}_i))=h(y(\mathbf{s}_i)) \exp (\theta(\mathbf{s}_i) y(\mathbf{s}_i)-b(\theta(\mathbf{s}_i))) \tag{21}
$$
其中，使用规范链接，$\theta(\mathbf{s}_i)=\mathbf{X}^T(\mathbf{s}_i) \tilde{\boldsymbol{\beta}} (\mathbf{s}_i)$.该规范生成 Diggle 等人的模型。 （1998 年）。在 (19) 中，我们可以包含一个几乎没有额外复杂性的离散度参数。
由此产生的第一阶段似然变为
$$
L (\tilde{\boldsymbol{\beta}}: \mathbf{y})=\exp {\sum y(\mathbf{s}_i) \mathbf{X}^T(\mathbf{s}_i) \tilde{\boldsymbol{\beta}}(\mathbf{s}_i)-b(\mathbf{X}^T(\mathbf{s}_i) \tilde{\boldsymbol{\beta}}(\mathbf{s}_i))} \tag{22}
$$
选择 式（14） 中 $\tilde{\boldsymbol{\beta}}$ 上的先验, 贝叶斯模型可以在 $\phi$、$T$ 、$\mu_{\boldsymbol{\beta}}$ 的先验基础上实现。
可以使用 Gibbs 采样器形式的概念上简单的马尔可夫链蒙特卡罗算法来拟合该模型，该算法将使用自适应拒绝采样（Gilks 和狂野 1992）。在$T$ 上使用逆Wishart，得到的$T$ 的完整条件再次是逆Wishart。更新 $\phi$ 通常是很麻烦的，因为它在 (14) 中以 Kronecker 形式输入。 Metropolis 更新很难设计，但可能提供了最好的似然。同样有问题的是 $\tilde{\boldsymbol{\beta}}$ 的重复组件更新。这种分层居中的参数化（Gelfand, Sahu, and Carlin 1995, 1996）比使用 $\mu_{\boldsymbol{\beta}}$ 和 $\boldsymbol{\beta}$ 更可取，但该算法在自相关问题上仍然运行得很慢。
组件更新的替代方法是通过 Langevin 扩散引入阻塞的 Metropolis 更新（参见，例如，Roberts、Gelman 和 Gilks 1997；Christensen、Möller 和 Waagepetersen 2000）。更新整个 $n p \times 1$ 向量 $\tilde{\boldsymbol{\beta}}$ 似乎不切实际。通过组件进程阻止 $\tilde{\boldsymbol{\beta}}$ 在计算上本质上是难以处理的，并且通过站点阻止 $\tilde{\boldsymbol{\beta}}$ 会产生强自相关问题。显然，需要做更多的工作来为这些模型提供有效的模型拟合策略。
8 大 N 问题
可参见高斯过程的大 n 解决方案。
9 软件包
（1）spBayes
Finley 等（2015 年）专注于制定计算高效且灵活的 MCMC 算法，以估计一系列时空高斯过程 (GP) 模型。但尽管所提出的采样算法非常通用，但在 spBayes 软件包中只实现了很有限的一组模型。具体来说，用户只能在模型的截距项参数上设定单变量或多变量的高斯过程，而无法将其扩展到所有系数（即空间变系数模型）。现在，通过将 spSVC 函数添加到 spBayes 软件包，使用户可以在任何模型回归系数集合上设置高斯过程，从而为用户提供更多选项。
（2）相关文献
文献情况：
- 通用变系数模型方面：
- Hastie, T., Tibshirani, R., 1993. Varying-coefficient models. J. R. Stat. Soc. Ser. B 55 (4), 757–796.
- Fan, J., Zhang, W., 2008. Statistical methods with varying coefficient models. Stat. Interface 1 (1), 179–195.
- 空间变系数模型方面：
- 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), 387–396.
（3）类似软件
此类功能并非独一无二，已经有几个 R 包能够拟合空间变系数模型。其中有一些专门用于处理空间或时空数据，而另一些则面向更为通用的可变化变量集，只是这些变量可以是坐标系中的索引。这些包大多采用了一些基于样条或核的回归方法，以允许预测变量的不同影响。
- 地理加权回归类：
- spgwr 包：该软件包（Bivand 和 Yu，2017）实现了地理加权回归，如 Fotheringham 等（2002 年）提出的那样。
- 样条回归类：
- mgcv 包（Wood，2017）
- svcm 包（Heim，2007）
- np 包（Hayfield 和 Racine，2008）
- mboost 包（Hothorn 等，2018）。
- 基于树的方法：
- vcrpart 包：Bürgin 和 Ritschard (2017) 最近开发了一种基于树的变系数模型 (TVCM) 算法。
- 时空贝叶斯方法：
- walker 包 (Helske, 2019; Vihola et al., 2017)
- spTimer 包（Bakar and Sahu， 2018)
- spTDyn 包（Bakar et al., 2017, Bakar et al., 2015a, 2015b; 2015b)
- 其他使用通用贝叶斯软件的方法
- INLA (Rue et al., 2009; Lindgren and Rue, 2015; Bakka et al., 2018)
- Stan (Carpenter et al., 2017; Stan Development Team, 2018)
