空间过程的贝叶斯建模分析方法综述

【阅读建议】本文重点介绍点参考空间数据的贝叶斯建模和分析方法，尤其是贝叶斯分层建模框架。点参考数据（也被称为地统计数据）主要指在固定空间位置观测到的随机变量数据。过去二十年中，此类数据在空间和时间上的收集量已经大大增加，随之而来的是分析此类数据的大量方法。本文尝试对其中的贝叶斯方法进行回顾。此类分析方法的好处是能够进行全面而准确的推断，并对不确定性进行适当评估。地统计建模的测站数据虽然比较复杂，涉及单变量和多变量、连续型和类别型、静态和动态以及大量长时间观测结果等，但在贝叶斯分层模型框架内，可以统一进行描述和阐释。本文另一亮点在于对大规模观测数据的建模问题做了综述，介绍了降秩方法（高斯预测过程模型）和近邻方法（近邻高斯过程模型）两类主要的处理策略。

【引文信息】 A. E. Gelfand and S. Banerjee, “Bayesian Modeling and Analysis of Geostatistical Data,” Annu Rev Stat Appl, vol. 4, pp. 245–266, 2017, doi: 10.1146/annurev-statistics-060116-054155.

1. 引言

在过去的二十年里，空间分析以惊人的速度增长，主要是由于廉价、高速计算的可用性增加（就像许多其他领域的情况一样）。这使得在许多领域都能够收集到大型空间和时空数据集，复杂地理信息系统软件得到广泛使用，以及探究（拟合和推断）具有挑战性的、更适合、更现实模型的能力。作为结果，空间统计方法已经从一个似乎有点特殊的领域，变成了一个更受模型驱动的领域。

空间数据通常被视为以下三种类型之一：点参考数据、面元数据和随机空间点模式数据。本文只考虑点参考数据的情况，亦即地统计数据模型。本文大量文献内容将聚焦在分层建模上，因为我们认为，在地统计数据分析方面，分层模型是当前也是将来的最主要的前沿技术。

通常，地统计建模可以被视为一种分层定义，因此自然会导致采用贝叶斯框架进行推断，并采用合适的吉布斯采样/马尔可夫链蒙特卡罗 (MCMC) 等方法进行模型拟合。这并不奇怪，因为更普遍地，分层建模已经成为 21 世纪科学工作的建模范式（Clark & Gelfand 2006 ^[10]），尤其是空间统计。本综述将探讨如下地统计建模问题： - （1）空间线性地统计模型（含克里金法） - （2）空间广义线性地统计模型 - （3）多变量地统计模型 - （4）多源空间数据融合方法 - （5）时空数据的建模 - （6）时空大数据问题 ## 2. 单变量建模基本地统计模型将位置 $\mathbf{s}$ 处的响应 $Y(\mathbf{s})$ 指定为： $$ Y(\mathbf{s}) = \mathbf{x}^T(\mathbf{s}) \boldsymbol{\beta} + w(\mathbf{s}) + \epsilon(\mathbf{s}) \tag{1} $$ 其中： - $\mathbf{x}(\mathbf{s})$ 是空间参考协变量/预测变量的 $p\times 1$ 向量，回归系数 $\boldsymbol{\beta}$ 也是一个 $p\times 1$ 向量。例如，如果 $Y(\mathbf{s})$ 是位置 $\mathbf{s}$ 的温度，则 $\mathbf{x}(\mathbf{s})$ 可能是包括海拔等预测变量的向量。 - 残差分为两个部分： - 空间相关部分： $w(\mathbf{s})$ 是地统计的主体，通常被指定为方差为 ${\sigma}^2$ 的平稳高斯过程 (GP)； - 非空间相关部分：$\epsilon(\mathbf{s})$ 通常是一个方差为 ${\tau}^2$ 的加性白噪声，即所谓的块金效应（参见 Cressie & Wikle 2011 ^[12] ，Banerjee 等，2014 ^[2] ）。 - ${\tau}^2$ 和 ${\sigma}^2$ 是方差的两个分量，并且白噪声 $\epsilon(\mathbf{s})$ 的空间无关性会导致 $Y(\mathbf{s})$ 表面处处不连续。 ### 2.1 分层建模法下面以地统计学中常使用的高斯过程为例，来说明传统贝叶斯建模与分层建模之间的区别。 **（1）传统的贝叶斯建模方法** 如果我们有数据 $Y (\mathbf{s}_i ), i = 1, \ldots , n$，并令向量 $\mathbf{Y} = (Y (\mathbf{s}_1), \ldots , Y (\mathbf{s}_n))^T$。假设 $w(\mathbf{s})$ 对应于高斯过程，其具有二阶平稳的协方差函数 $\operatorname{cov}(w(\mathbf{s}), w(\mathbf{s}^{\prime})) = {\sigma}^2{\rho}(\mathbf{s} − \mathbf{s}^{\prime}; \phi)$，式中 ${\sigma}^2$ 是前面提到的空间相关的方差，${\rho}$ 是一个被 $\phi$ 参数化了的相关性函数，同时令独立同分布的白噪声为 $\epsilon(\mathbf{s}) \sim \mathcal{N}(0, {\tau}^2)$ 则根据高斯过程原理，我们可以得到点 $\mathbf{s}_i$ 处的边缘化分布 ${w(\mathbf{s}_i)}$，其协方差矩阵为 $\Sigma = {\sigma}^2H(\phi) + {\tau}^2I$。式中的 $H(\phi)_{ij}= \operatorname{corr}(w(\mathbf{s}_i),w(\mathbf{s}_j))$，即 $H(\phi)$ 的第 $(i, j)$ 个元素对应于 $w( \mathbf{s}_i)$ 和 $w(\mathbf{s}_j)$ 的相关系数。可以看出，协方差矩阵 $\Sigma$ 是对角占优的，因此在计算似然所需的求逆和求行列式方面表现良好。根据 `第 2.1 节` 的基础形式，结合高斯过程，可以得到模型的似然为： $$ \mathbf{Y} \mid \boldsymbol{\theta} \sim \mathcal{N}(X \boldsymbol{\beta}, \sigma^2H(\phi) + \tau^2I) \tag{2} $$ 其中 $X$ 是 $n\times p$，将 $\mathbf{x}(\mathbf{s}_i)$ 采集为行向量，$\boldsymbol{\theta}$ 是模型中所有参数的集合。在贝叶斯设置中，通常进一步假设不同参数之间相互独立，进而有： $p(\boldsymbol{\theta}) = p(\boldsymbol{\beta}) p({\sigma}^2) p({\tau}^2) p(\phi)$。举例来说，我们可以假设相关函数 ${\rho}(·)$ 为各向同性的，是距离 $\|\mathbf{s} - \mathbf{s}^{\prime} \|$ 的函数，衰减参数为 $\phi$，则可以作出如下参数先验假设： - $\boldsymbol{\beta}$：先验习惯于选择多元高斯分布（甚至平坦分布） - ${\sigma}^2$ 和 ${\tau}^2$：选择逆伽玛分布作为先验 - $\phi$：先验取决于对 ${\rho}$ 的选择，常用设置为均匀分布或伽马分布。在实际工作中，这里的 ${\sigma}^2$、${\tau}^2$ 和 $\phi$ 不能有不恰当的先验或后验结果（Berger 等， 2001 ^[5] ）。 **（2）分层建模思想** 如果将上述似然重写为以空间随机效应为条件的形式，对于建模和科学解释是非常有帮助的。也就是说，我们可以将上述模型重塑为等价的分层方式： - 第一阶段（数据建模）：$\mathbf{Y} \mid \boldsymbol{\theta}, \mathbf{w} \sim \mathcal{N}(X \boldsymbol{\beta} + \mathbf{w}, {\tau}^2 I)$。此阶段重点考虑的问题是，如果空间过程 $w(\mathbf{s}_i)$ 已知，那么观测数据是怎么产生的？通常我们假设观测 $Y(\mathbf{s}_i)$ 是独立同分布的。 - 第二阶段（过程建模）：$\mathbf{w} \mid {\sigma}^2, \phi \sim \mathcal{N}(\mathbf{0}, {\sigma}^2H(\phi))$。此阶段重点考虑的是空间随机过程（或随机效应模型），可以有多种空间过程模型，此处沿用高斯过程作为第二阶段的过程模型。 - 第三阶段（参数建模）：$(\boldsymbol{\beta}, {\tau}^2, {\sigma}^2, \phi)$ 。此阶段重点考虑的是前两个阶段的模型参数该如何设置先验？基础概率分布通常可以参数化指定，那么这些参数该如何设置？在贝叶斯方法中，有关参数分布的参数被称为超参数，因此第三阶段通常也被称为超先验指定阶段。上述分层模型和 `式（2）` 的边缘模型在推断上是等价的。传统上，人们可能更习惯于直接拟合 `式（2）` 的边缘模型（参见 Banerjee 等，2014 ^[2] ）。但在本文中，`分层模型（Hierarchical Model）` 形式作为 `多层次过程模型（Multilevel Process Model）` 的一种特例，将会贯穿全文并体现出更为丰富的科学内涵（参见 Berliner 2000 ^[6]）。 $$ [\text{data} \mid \text{process}, \text{parameters}][\text{process} | \text{parameters}] [\text{parameters}] $$ 在过去二十年中，这种 `多层次建模（Multilevel Modelling）` 的科学案例比比皆是。 Banerjee 等 (2014 ^[2]) 展示了一系列带有相关参考资料的演示案例，其中包括使用房屋和社区特征对房价进行建模、使用环境回归模拟物种的存在/不存在或丰度、环境污染物联合建模及气象变量解释等上述所有示例都显示出了强烈的空间依赖性，并且空间方差 ${\sigma}^2$ 支配着纯误差 ${\tau}^2$ 的贡献。

一种解释是 $\sigma^2$ 体现了中尺度的不确定性，而 $\tau^2$ 体现了微尺度和测量误差带来的不确定性。

### 2.2 经典克里金预测克里金问题（ Krige 1951 ^[33] ）是一种最优空间预测：给定一个随机场的观测值 $\mathbf{Y} = (Y(\mathbf{s}_1), \ldots, Y(\mathbf{s}_n))^T$，如何最好地预测出测站 $\mathbf{s}_0$ 处未被观测到 $Y$ 值？由此可见，普通克里金法只关注 $Y$ ，并不关注其与解释变量 $X$ 之间的关系，或者认为该关系已知。普通克里金法基于线性形式（ $\Sigma \ell_i Y(\mathbf{s}_i) + δ_0$ ）的无偏性原则和最小均方误差原则来获得 $Y(\mathbf{s}_0)$ 的最优无偏线性预测。该方法基于本征平稳性假设，假设偏差的均值为零（ $\mathbb{E}[Y(\mathbf{s} + h) - Y(\mathbf{s})] = 0$），并且偏差的方差只与步长（也称间距、滞后距离等）有关。也就是说，该空间随机场存在仅与步长有关的变异函数 ${\boldsymbol{\gamma}}(h) ≡ E[Y(\mathbf{s} + h) - Y(\mathbf{s})]^2$。基于此假设，克里金法可以明确地将最优线性预测表示成 ${\boldsymbol{\gamma}}(h)$ 的函数。普通克里金方法不再需要其他假设，只是在进行最优估计前，不许先利用观测数据估计（半）变异函数 ${\boldsymbol{\gamma}}(h)$。在 `式（2）` 的边缘分布模型中，普通克里金法寻求能够最小化均方预测误差 $\mathbb{E}[(Y(\mathbf{s}_0) - h(\mathbf{Y}))^2|\mathbf{Y}]$ 的函数 $h(\mathbf{Y})$，从而有无偏估计 $h (\mathbf{Y}) = \mathbb{E}[Y(\mathbf{s}_0)|\mathbf{Y}]$。因此 $Y(\mathbf{s}_0)$ 的后验均值和方差具有如下明确的形式： $$ \mathbb{E}\left[Y\left(\mathbf{s}_0\right) \mid \mathbf{Y}\right]=\mathbf{x}\left(\mathbf{s}_0\right)^T \boldsymbol{\beta}+\boldsymbol{\gamma}^T \Sigma^{-1}(\mathbf{Y}-X \boldsymbol{\beta}) \tag{3} $$ $$ \operatorname{var}\left[Y\left(\mathbf{s}_0\right) \mid \mathbf{Y}\right]=\sigma^2+\tau^2-\boldsymbol{\gamma}^T \Sigma^{-1} \boldsymbol{\gamma} \tag{4} $$ 在各向同性假设下，$\boldsymbol{\gamma}^T=\left(\sigma^2 \rho\left(d_{01} ; \phi\right), \ldots, \sigma^2 \rho\left(d_{0 n} ; \phi\right)\right)$ ，其中 $d_{0 j}=\left\|\mathbf{s}_0-\mathbf{s}_j\right\|$。普通克里金法中的后验均值使用了一个线性预测器，并假定该线性预测器的参数 $\boldsymbol{\beta}$ 已知。如果该线性预测器自身的参数也需要被估计时，情况就不再如此了，需要采用该完全贝叶斯克里金法。 ### 2.3 贝叶斯克里金预测 **（1）完全贝叶斯克里金模型** 贝叶斯克里金法的建模对象，是在给定 $\mathbf{s}_0$ 处的预测变量 $\mathbf{x}(\mathbf{s}_0)$ 以及观测数据 $\mathbf{Y}$ 和 $\mathbf{X}$ 的情况下，推断目标变量 $Y(\mathbf{s}_0)$ 的预测性分布： $$ p(Y(\mathbf{s}_0)|\mathbf{x}(\mathbf{s}_0),\mathbf{Y}, X) = \int p(Y(\mathbf{s}_0)|\mathbf{x}(\mathbf{s}_0),\mathbf{Y}, \boldsymbol{\theta}) p(\boldsymbol{\theta} \mid \mathbf{Y},\mathbf{X}) d \boldsymbol{\theta} \tag{5} $$ 其中积分项 $p(Y(\mathbf{s}_0)|\mathbf{x}(\mathbf{s}_0),\mathbf{Y}, \boldsymbol{\theta})$ 是由 $Y(\mathbf{s}_0)$ 和原始数据 $\mathbf{Y}$ 的联合多元高斯分布产生的高斯分布。 **（2）预测分布的计算** 可以很容易地通过模拟方法（或称采样方法）获得 `式 5` 的近似估计。假设我们从模型参数的后验分布 $p(\boldsymbol{\theta} |\mathbf{Y}, X)$ 中抽取后验样本 $\boldsymbol{\theta}^{(1)}, \boldsymbol{\theta}^{(2)}, \ldots, \boldsymbol{\theta}^{(B)}$。那么 `式 5` 可以近似为以下形式的蒙特卡洛求和方式： $$ \widehat{p}\left(Y\left(\mathbf{s}_0\right) \mid \mathbf{x}\left(\mathbf{s}_0\right), \mathbf{Y}, X\right)=\frac{1}{B} \sum_{b=1}^B p\left(Y\left(\mathbf{s}_0\right) \mid \mathbf{x}\left(\mathbf{s}_0\right), \mathbf{Y}, \boldsymbol{\theta}^{(b)}\right) \tag{6} $$ 在实际工作中，我们使用组合采样方式来抽取 $Y\left(\mathbf{s}_0\right)$ 的预测样本，每个 $\boldsymbol{\theta}^{(b)}$ 对应一个样本点，即 $Y\left(\mathbf{s}_0\right)^{(b)}$ 是 $p\left(Y\left(\mathbf{s}_0\right) \mid \mathbf{x}\left(\mathbf{s}_0\right), \mathbf{Y}, \boldsymbol{\theta}^{(b)}\right)$ 的一个样本。这样产生的集合 $\{Y(\mathbf{s}_0)^{(1)}, Y(\mathbf{s}_0)^{(2)}, \ldots, Y(\mathbf{s}_0)^{(B)},\}$ 作为后验预测分布的随机样本集，这使直方图或核密度估计能够得到近似，避免了 `式 6` 的复杂蒙特卡洛求和。

关于参数推断

预测分布的计算基于对协方差函数（或半变异函数）参数的准确推断，而现代的贝叶斯工具箱提供了大量计算效率高的手段，如： MCMC 采样方法、变分推断方法等详情参见相关文献，在此不赘述。

此处采用贝叶斯方法实现推断，和传统地统计学中的泛克里金法不同。总体来说，传统克里金方法虽然基于高斯过程模型，但预测器的拟合采用频率派的最优化方法实现，与本文的贝叶斯主旨有所不同。

### 2.4 一个说明性的例子我们使用来自美国农业部林务局的森林调查数据集进行说明。该数据集是 `R` 中 `spBayes` 包的一部分，包含 1991 年和 2002 年 $437$ 个地块的森林调查数据。我们将总森林生物量的对数变换作为结果变量，将其回归到五个预测变量：坡度、海拔和缨帽变换的亮度、绿度和湿度分量。我们使用 `R` 包 `spBayes` 中的 `spLM` 函数执行完全贝叶斯克里金法，如 `第 2.3 节` 所述。该函数可从 [CRAN](https://cran.r-project.org/web/packages/spBayes/index.html) 获得。有关数据集和代码的更多详细信息，请参见 Banerjee 等（2014 ^[2]）的 `第 6.3 节`，完整的数据集可在 spBayes 中获得。我们拟合了 `第 2.1 节` 中描述的那个高斯过程的贝叶斯分层模型： - 数据模型：回归系数 $\boldsymbol{\beta}$ 采用平坦先验 - 过程模型：相关性函数采用指数函数 ${\rho}(\mathbf{s} − s^{\prime}; \phi) = \exp(−\phi \|s − s^{\prime} \|)$ - 参数模型：${\sigma}^2$ 和 ${\tau}^2$ 分别采用居中逆伽玛先验 $IG(2, 0.08)$ 和 $IG(2, 0.02)$， $\phi$ 采用均匀先验 $U(3/1500, 3/50)$，建议基于 $50$ 到 $1500$ 米的近似范围研究领域的维度。实验使用标准 MCMC 算法获得模型参数的后验样本后，使用 `式 5` 和 `式 6` 进行贝叶斯克里金计算。`图 1` 显示了输出的插值后预测表面，即研究域上的对数(公吨)生物量。 ![单变量展示](https://xishansnowblog.oss-cn-beijing.aliyuncs.com/images/images/stats-20221110173750-02ac.webp) > 图 1：从 Bartlett 实验森林数据中观测到的对数(公吨)生物量的插值后表面。 ### 2.5 避免不恰当的渐近一维时域的推断与二维空间域的推断截然不同：在前者中具有全序，但在二维（或更多）中只有偏序。这暗示着：时间序列的大样本分析通常会设想一个渐进增加的时间域，而且时间域趋向于 $\infty$。相比之下，地统计数据的大样本分析通常会设想一个固定区域，并且在该区域中填充越来越多的点，即所谓的渐近内填充（ infill asymptotics），参见 Stein 1999a ^[47],b ^[48]；Loh 2005 ^[35] ；Zhang & Zimmerman 2005 ^[60] 。当针对时域使用递增渐进结果时，通常假设随着收集到更多数据，我们可以了解到随时间间隔增加的时态关联特性；但当我们应用渐近内填充结果时，由于在研究域内最大距离是固定的，我们无法了解更远距离的关联，因此只能将指定位置的空间预测做得越来越好。无论如何，我们对过程参数的推断将是有界的；Fisher 信息不会无穷大，Cramér-Rao 下界和渐近方差也不会达到零。上述想法鼓励了在贝叶斯框架中对模型的精确推断，并尽量规避任何用于推断的渐近理论。当然，“信息有界” 这一事实意味着，数据永远无法像贝叶斯假设中那样压倒先验。没有免费的午餐，可能需要先验的敏感性分析。 ### 2.7 非高斯变量建模有时候，数据集本身可能不支撑第一阶段采用高斯建模，因为高斯分布要求 $Y(\mathbf{s})$ 是连续的，但实际工作中，存在很多离散型的随机变量。例如，$Y(\mathbf{s})$ 可以是观测位置处的二进制或计数变量（定序型变量或计数型变量），也可能是提供存在/不存在或丰度的物种信息（名义型变量）。此时，我们一般会将第一阶段数据模型中的高斯分布替换为某种指数族分布（有可能带有离散度参数）；而第二阶段的过程模型（通常是高斯过程）依然保持不变（参见 Diggle 等， 1998 ^[18]）。此类模型被成为空间广义线性模型（Spatial Generalized Linear Models，SGLM）。现在有： - **第一阶段**：给定 $\boldsymbol{\beta}$ 和 $w\left(\mathbf{s}_i\right)$ 时，$Y\left(\mathbf{s}_i\right)$ 是条件独立的（独立观测假设） - $f\left(y\left(\mathbf{s}_i\right) \mid \boldsymbol{\beta}, w\left(\mathbf{s}_i\right), \gamma\right)=h\left(y\left(\mathbf{s}_i\right), \gamma\right) \exp \left \{\gamma\left[y\left(\mathbf{s}_i\right) \eta\left(\mathbf{s}_i\right)-\Psi\left(\eta\left(\mathbf{s}_i\right)\right)\right]\right\}$ - $f(\cdot)$ 是一个链接函数，$f\left( \mathbb{E}\left[Y \left(\mathbf{s}_i\right)\right]\right)=\eta\left(\mathbf{s}_i\right)=\mathbf{x}^T\left(\mathbf{s}_i \right) \boldsymbol{\beta} + w\left(\mathbf{s}_i\right)$ - $\gamma$ 是一个离散度参数。 - **第二阶段**：$w(\mathbf{s})$ 为一个高斯过程: $\mathbf{w} \sim \mathcal{N}\left(\mathbf{0}, \sigma^2 H(\phi)\right)$。 - **第三阶段**：提供模型参数的超先验。

第一阶段模型没看懂，本身不复杂，为何表示的如此复杂？

现在，我们拥有的不再是随机过程 $Y(\mathbf{s})$ （相对于`第 2.1 节`的三阶段模型而言）,只有一个有效的有限维联合分布。但这对模型推断来说没有任何困难，因为我们只需要为有限数量的点拟合模型，并且只在有限数量的点处进行插值。将纯误差项 $\epsilon(\mathbf{s})$ 添加到链接变换的均值定义是不明智的，它们会随着第一阶段的条件独立而变得多余。此外，第二阶段建模的空间随机效应促使邻近位置的空间变量的均值彼此相似。对于位置 $\mathbf{s}$ 和 $\mathbf{s}^\prime$，在 $Y(\mathbf{s})$ 和 $Y(\mathbf{s}^\prime)$ 之间引入了边缘空间依赖性。 ## 3. 多变量地统计建模多变量点参考空间数据的一个常见示例来自环境监测站，这些站可能会同时对臭氧、一氧化二氮、一氧化碳、PM2.5 等进行测量。有了这些数据，我们可以预测特定位置以及跨位置测量值之间的依赖性（跨位置的不同属性随机变量之间有时也可能存在依赖性，例如一个位置的属性 A 受周边其他位置的属性 B、C 等影响的场景）。 ### 3.1 传统建模 **（1）建模思路** 令 $\mathbf{Y}(\mathbf{s})$ 表示 $\mathbf{s}$ 处随机变量的 $p\times 1$ 向量。我们感兴趣的是 $\{\mathbf{Y}(\mathbf{s}): \mathbf{s} \in D\}$。我们设想表面之间具有依赖性的 $p$ 个空间表面，因此需要一个 $p$ 维随机过程。高斯过程似乎是一个自然的选择，但对于 $n$ 个位置来说，我们需要一个 $np$ 维的多元高斯分布。关键对象是互协方差，$C(\mathbf{s}, \mathbf{s}^\prime) = \operatorname{cov}(\mathbf{Y}(\mathbf{s}),\mathbf{Y}(\mathbf{s}^\prime))$，这是一个不需要对称的 $p\times p$ 矩阵，即任意两个测站处的任意两个变量相互之间的相关性不一定相等， $\operatorname{cov}(Y_j(\mathbf{s}), Y_{j^\prime} (\mathbf{s}^\prime))$ 不必等于 $\operatorname{cov}(Y_{j^\prime} (\mathbf{s}), Y_j(\mathbf{s}^\prime))$。此外，任意两个测点之间的互协方差矩阵 $C(\mathbf{s}, \mathbf{s}^\prime )$ 不一定是正定的（除非在限制意义上）； $C(\mathbf{s}, \mathbf{s})$ 是与 $\mathbf{Y}(\mathbf{s})$ 相关联的自协方差矩阵。如果继续结合 $p$ 个变量的均值，我们将具有所需的 $np$ 维高斯分布。

注解：此处用互协方差矩阵（或交叉协方差矩阵）来对任意两个测点处的多种随机变量之间的相关性建模。当 $\mathbf{s} == \mathbf{s}^\prime$ 时，$C(\mathbf{s},\mathbf{s})$ 代表某个测站内部各随机变量之间的协方差，因此称为自协方差。

**（2）共区域化模型** 有大量为 $C(\mathbf{s}, \mathbf{s}^\prime)$ 创建有效定义的文献（例如，参见 Chiles & Delfiner 1999 ^[8]）。一个简单的建设性（有效）选择是 `共区域化（coregionalization）` 的线性模型，最初由 Wackernagel（1994 ^[54]，2003 ^[55]）建议用于降维，随后由 Gelfand 等 (2004 ^[25]) 提出用于多元空间建模。该模型定义如下：令多元随机过程 $\mathbf{v}(\mathbf{s}) = A \mathbf{w}(\mathbf{s})$，其中 $\mathbf{w}(\mathbf{s}) = (w_1(\mathbf{s}),w_2(\mathbf{s}), \ldots ,w_p(\mathbf{s}))$ ， $w_j(\mathbf{s}),j= 1, 2, \ldots , p$ 指 $p$ 个独立的平稳空间过程，具有各自的相关函数 $\rho_j(\mathbf{s} − \mathbf{s}^\prime)$。$A$ 是任意 $p\times p$ 非奇异矩阵。此多元随机过程中的其互协方差矩阵可以被定义为 $C(\mathbf{s}-\mathbf{s}^\prime)=\sum_{j=1}^{\rho} \rho_j(\mathbf{s}-\mathbf{s}^\prime)\mathbf{a}_j \mathbf{a}_j^T$ （其中 $\mathbf{a}_j$ 表示矩阵 $A$ 的一列）。该线性变换为各空间过程之间引入了依赖性，并且还为每个分量提供了不同的协方差函数。一种特殊的情况是：对于所有 $j$ 都具有共同的相似函数 $\rho (\mathbf{s} − \mathbf{s}^\prime )$ ，此时互协方差矩阵可以写为可分离的形式：$C(\mathbf{s} − \mathbf{s}^\prime) = \rho (\mathbf{s} − \mathbf{s}^\prime)AA^T$。由此，我们得出了一个通用的 `共区域化多元地统计模型`： $$ \mathbf{Y}(\mathbf{s}) = \boldsymbol{\mu}(\mathbf{s}) + \mathbf{v}(\mathbf{s}) + \boldsymbol{\varepsilon }(\mathbf{s}) $$ 其中: - $\boldsymbol{\varepsilon}(\mathbf{s}) \sim \mathcal{N}(\mathbf{0}, \mathbf{D}_{\varepsilon}), (D_{\varepsilon})_{jj}=\tau_j^2$ - $\mathbf{v}(\mathbf{s}) = A \mathbf{w}(\mathbf{s})$，其中 $w_j(\mathbf{s})$ 是均值为零且具有独立相关函数的高斯过程 - $\boldsymbol{\mu}(\mathbf{s})$ 通常以回归形式出现，即 $\mu_j(\mathbf{s})=\mathbf{x}_j^T(\mathbf{s})\boldsymbol{\beta}_j$ 形式。可以看到，上式使我们可以得到一个以多元空间过程模型 $\mathbf{v}(\mathbf{s})$ 为条件的分层模型。 **（3）将空间变系数模型视为多变量建模** 空间变系数模型 (Gelfand 等， 2003 ^[24] ) 是展示层次建模能力的另一个有吸引力的例子。假设将 $\mathbf{s}$ 处的单变量观测值 $Y(\mathbf{s})$ 建模为 $Y(\mathbf{s}) = \mathbf{x}(\mathbf{s})^T\boldsymbol{\beta}(\mathbf{s}) + \varepsilon (\mathbf{s})$。则在空间变系数模型中，多元空间过程是针对系数 $\boldsymbol{\beta}(\mathbf{s})$ 的。此处假设 $p = 2$ 且 $\mathbf{x}(\mathbf{s})$ 中有一列为 $1$，则有： $$ Y(\mathbf{s}) = \beta_0 (\mathbf{s}) + x(\mathbf{s})\beta_1 (\mathbf{s}) + \varepsilon (\mathbf{s}) $$ 式中有一个随空间变化的截距项（如空间随机效应）和一个随空间变化的斜率。由于希望斜率和截距能够和位置相关，因此，我们需要一个双变量（截距和斜率）的空间过程，进而可以获得一类非常丰富的非线性模型（但对于每一个空间点局部是线性的）。值得注意的是，我们可以在仅观测单变量 $Y(\mathbf{s})$ 空间过程的情况下，推断多变量的高斯过程 $\boldsymbol{\beta}(\mathbf{s})$。只是，由于我们只有 $n$ 个位置的观测向量，因此无法详细了解多元空间过程 $\boldsymbol{\beta}(\mathbf{s})$ 中的每一个组分。在实际工作中，我们经常会发现其中截距组分的过程表面是信息性的，而斜率组分的过程表面往往是平坦的。 ### 3.2 核卷积方法另一种创建多元地统计模型的方法是通过移动平均 (Ver Hoef & Barry 1998 ^[53] )，也称为核卷积 (Higdon 1998 ^[28])。事实上，这是创建平稳随机过程的一般性方法的特例（Yaglom 1987 ^[59] ，第 8 章）。假设 $k_l(·)$，$l= 1, \ldots, p$ 是 $\mathbb{R}^2$ 上的一组 $p$ 个平方可积的核函数，且 $k_l(0) = 1$。令 $w(\mathbf{s})$ 是一个均值为零、方差为 $1$ 的高斯过程，其相关函数还是表示为 $\rho$，则我们可以定义 $p$ 元空间过程 $\mathbf{v}(\mathbf{s})$ 的每一个组分空间过程为： $v_l(\mathbf{s}) = \sigma_l∫k_l(\mathbf{s} − \mathbf{t})w(\mathbf{t})d \mathbf{t}, l= 1, \ldots, p, \sigma_l> 0$。因此，$\mathbf{v}(\mathbf{s})$ 是均值为 $\mathbf{0}$、互协方差函数为 $C(\mathbf{s}, \mathbf{s}^\prime )$ 的多元高斯过程，其中互协方差函数定义如下： $$ \left(C\left(\mathbf{s}, \mathbf{s}^{\prime}\right)\right)_{ll^{\prime}}=\sigma_l \sigma_{l^{\prime}} \iint k_l(\mathbf{s}-\mathbf{t}) k_{l^{\prime}}\left(\mathbf{s}^{\prime}-\mathbf{t}^{\prime}\right) \rho\left(\mathbf{t}-\mathbf{t}^{\prime}\right) \mathrm{d} \mathbf{t} \mathrm{d} \mathbf{t}^{\prime} \tag{7} $$ 根据变量变化原理，$(C(\mathbf{s}, \mathbf{s}^\prime))_{ll^\prime}$ 只依赖于 $\mathbf{s} − \mathbf{s}^\prime$ ，也就是说，此时的 $\mathbf{v}(\mathbf{s})$ 是一个弱平稳过程。如果 $k_l$ 仅通过间隔长度 $\|s−\mathbf{s}^\prime\|$ 依赖于 $\mathbf{s}−\mathbf{s}^\prime$ 且 $\rho$ 是各向同性的，则 $C(\mathbf{s} − \mathbf{s}^\prime)$ 也是各向同性的。要使用此过程，通常会使用一组有限的节点位置做离散近似。具体来说，选择一组有限的位置 $\mathbf{s}^*_1,\ldots,\mathbf{s}^*_r$ 并定义每一个分量过程 $v_l(\mathbf{s})=\sigma_l \sum_{j=1}^r k_l(\mathbf{s}-\mathbf{s}^*_j)w(\mathbf{s}^*_j)$。此时 $(C(\mathbf{s}, \mathbf{s}^\prime ))_{ll^\prime}$ 是一个双求和表达式。由于离散近似，其结果过程不再平稳。 ### 3.3 协同克里金预测在经典多元地统计学中，协同克里金指的是在新位置 $\mathbf{s}_0$ 处对 $\mathbf{Y}(\mathbf{s}_0)$ 或分量 $Y_j(\mathbf{s}_0)$ 进行空间预测。也就是说，我们可以在给定所有观测向量的情况下，预测整个向量 $\mathbf{Y}(\mathbf{s}_0)$；我们可以根据观测到的向量以及剩余分量的值来预测新位置处的某个分量（协同克里金法）。读者可以在 Wackernagel (2003 ^[55]) 中找到详细的评论。与 `第 2.3 节` 类似，如果假设一个多元高斯空间过程，我们可以获得单个协同克里金问题的显式解。为简单起见，假设 $\mathbf{Y}(\mathbf{s})$ 将均值对齐为 $\boldsymbol{0}$。互协方差函数提供了互协方差矩阵 $\Sigma_{\mathbf{Y}}$，即数据 $\mathbf{Y} = (\mathbf{Y}(\mathbf{s}_1)^T,\mathbf{Y}(\mathbf{s}_2)^T, \ldots , \mathbf{Y}(\mathbf{s}_n)^T)^T$ 的 $np \times np$ 协方差矩阵。假设现在要预测 $Y_l(\mathbf{s}_0)$，则互协方差函数提供了一个 $np\times 1$ 向量 $\mathbf{C}_0$，多样本堆叠成 $\mathbf{C}_{0j},\quad j= 1, 2, \ldots ,n$ （其第 $l$ 个元素 $c_{0j, l}= \operatorname{cov}(Y_1 (\mathbf{s}_0), Y_l(\mathbf{s}_j))$ ）。通过 $(\mathbf{Y}, Y_1(\mathbf{s}_0))$ 的多元高斯分布，可以得到协同克里金估计： $$ \mathrm{E}\left(Y_1\left(\mathbf{s}_0\right) \mid \boldsymbol{Y}\right)=\mathbf{C}_0^T \Sigma_Y^{-1} \boldsymbol{Y} \tag{8} $$ 相关的方差 $\operatorname{var}(Y_1(\mathbf{s}_0)|\mathbf{Y})$ 也立即可用。类似的，我们可以开发协同位置的克里金形式；不过我们省略了细节。在上面的多元地统计模型中，我们可以实现类似于 `第 2.3 节` 的完全贝叶斯克里金法（详见 Banerjee 等， 2014 ^[2]）。 ### 3.4 地统计回归案例回归本身感兴趣的是条件分布 $Y(\mathbf{s}_0)|X(\mathbf{s}_0)$，而不是 $Y(\mathbf{s}_0)$ 的克里金预测。例如，我们可能会寻找树木直径和土壤湿度之间的局部关系。但在实际工作中，我们有时候不会像 `式 (1)` 那样直接建立一个地统计模型，而是采用多变量建模方法分别建立两个独立的随机过程，这是为什么呢？原因有以下两个： - 首先，我们也可能对逆回归问题也感兴趣，即给定 $Y(\mathbf{s}_0)$ 时的 $X(\mathbf{s}_0)$。因此，有必要对 $X(\mathbf{s})$ 进行随机建模。特别是当 $X(\mathbf{s})$ 为连续型时，我们需要为 $X$ 和 $Y$ 分别引入一个过程模型，从而形成一个多变量过程模型。 - 其次，我们需要考虑数据存在缺失的可能性（等效的，出现错位的情况）。也就是说，有一些 $X$ 没有关联的 $Y$，或有一些 $Y$ 没有关联的 $X$。我们希望使用所有可用数据，而不是只使用完整数据。为了进行基于模型的插补，我们需要一个针对缺失数据的双（多）变量模型（例如，Little & Rubin 2002 ^[34]）。当不存在缺失（错位）时，设 $\mathbf{X} = (X(\mathbf{s}_1), \ldots ,X(\mathbf{s}_n))^T$ and $\mathbf{Y} = (Y(\mathbf{s}_1), \ldots ,Y(\mathbf{s}_n))^T$ 分别是协变量和响应变量的测量值。我们设想（也可能经过适当的转换后）一个双变量高斯空间过程，$\mathbf{W}(\mathbf{s}) = (X(\mathbf{s}), Y(\mathbf{s}))^T$，其均值为 $\boldsymbol{\mu}(\mathbf{s}) = (\mu_X(\mathbf{s}),\mu_Y(\mathbf{s}))^T$ ，为了简单起见，采用前面提到的可分离互协方差函数定义。则 $Y(\mathbf{s})$ 和 $X(\mathbf{s})$ 在位置 $\mathbf{s}$ 的联合分布为： $$ \mathbf{W}(\mathbf{s})=\left(\begin{array}{l} X(\mathbf{s}) \\ Y(\mathbf{s}) \end{array}\right) \sim \mathrm{N}(\boldsymbol{\mu}(\mathbf{s}), T) \tag{9} $$ 当存在缺失（错位）的情况时，令 $\mathbf{X}$ 表示我们观测到的 $X(\mathbf{s})$ 向量，$\mathbf{Y}$ 为我们观测到的 $Y(\mathbf{s})$ 向量；另外令 $\mathbf{X}^\prime$ 表示相对于 $Y$ 缺失了 $X$ 的观测向量，$\mathbf{Y}^\prime$ 表示关于 $X$ 缺失了 $Y$ 的观测向量。根据前面的讨论，我们可以用增广向量 $\mathbf{X}_{\text{aug}}= (\mathbf{X},\mathbf{X}^\prime)$ 和 $\mathbf{Y}_{\text{aug}}= (\mathbf{Y},\mathbf{Y}^\prime)$ 来替换 $\mathbf{X}$ 和 $\mathbf{Y}$。通过对 $X$ 和 $Y$ 方向的排列，可以将其收集到向量 $\mathbf{W}_{\text{aug}}$ 中。在贝叶斯模型中，我们可以将未观测到的 $\mathbf{X}^\prime$ 和 $\mathbf{Y}^\prime$ 视为隐向量。当实现用于模型拟合的 Gibbs 采样时，我们在给定 $\mathbf{X}^\prime$ 和 $\mathbf{Y}^\prime$（即 $\mathbf{W}_{\text{aug}}$）的情况下先更新模型参数，然后在给定 $\mathbf{X}$、$\mathbf{Y}$ 和模型参数的情况下更新 $(\mathbf{X}^\prime,\mathbf{Y}^\prime)$。对于回归问题，我们感兴趣的分布是 $f(E[Y(\mathbf{s}_0)|X(\mathbf{s}_0)] | X(\mathbf{s}_0),\mathbf{Y},\mathbf{X})$。为简单起见，假设在测站处 $\boldsymbol{\mu}(\mathbf{s}) = (\mu_1,\mu_2)^T$ 保持不变。根据 `式 (9)` ，对于 $(X(\mathbf{s}), Y(\mathbf{s}))$ 来说, $p(Y(\mathbf{s})|X(\mathbf{s}), \beta_0 , \beta_1 , \sigma^2)$ 是一个高斯分布 $N(\beta_0 + \beta_1 X(\mathbf{s}), \sigma^2)$。即，$E[Y(\mathbf{s})|X(\mathbf{s})] = \beta_0 + \beta_1 X(\mathbf{s})$，其中 $$ \begin{align*} \beta_0 &=\mu_2-\frac{T_{12}}{T_{11}} \mu_1\\ \beta_1 &=\frac{T_{12}}{T_{11}} \\ \sigma^2&=T_{22}-\frac{T_{12}^2}{T_{11}} \tag{10} \end{align*} $$ 因此，给定来自 $(\mu_1,\mu_2, T,\phi )$ 的联合后验分布的样本，我们直接从 `式 (10)` 中的参数的后验分布中获取样本，进而从 $E[Y(\mathbf{s})|X(\mathbf{s})]$ 的后验分布中获取样本。显然，我们可以颠倒 $X$ 和 $Y$ 的角色来求解逆向问题。此外，多变量情况的扩展很清楚在地统计数据集中，面向二值型响应变量时的扩展可以按照如下方式进行。现在，在每个测站处，根据观测到的是 “失败” 还是 “成功”，$Z(\mathbf{s})$ 等于 $0$ 或 $1$。因此，该过程的实现称为双色图或二值图 (DeOliveira 2000 ^[17])。采用隐变量方法进行概率建模（例如，DeOliveira 2000 ^[17] ），让 $Y(\mathbf{s})$ 成为与位置相关的隐空间过程，让 $X(\mathbf{s})$ 成为标量的协变量（同样，很容易扩展到多变量情况）。当且仅当 $Y(\mathbf{s}) > 0$ 时，让 $Z(\mathbf{s}) = 1$。我们设想一个分布如 `式 (9)` 的双变量空间过程 $W(\mathbf{s}) = (X(\mathbf{s}), Y(\mathbf{s}))^T$，但现在 $\mu( s) = (\mu_1, \mu_2 + \boldsymbol{\alpha}^T \boldsymbol{U}(\mathbf{s}))^T$，其中 $\boldsymbol{U}(\mathbf{s})$ 是一个 $p\times 1$ 的固定协变量向量。给定 $X(\mathbf{s})$ 时 $Y(\mathbf{s})$ 的条件均值表明，相对于回归系数，给定 $X(\mathbf{s})$ 时 $Y(\mathbf{s})$ 的条件方差不可识别。因此，不失一般性，我们设置 $T_{22} = 1$，使得 $T$ 矩阵只有两个参数。最后，probit 空间回归模型为： $$ \begin{aligned} &P\left(Z(\mathbf{s})=1 \mid x(\mathbf{s}), \mathbf{U}(\mathbf{s}), \boldsymbol{\alpha}, \mu_1, \mu_2, T_{11}, T_{12}\right) \\ &=\Phi\left(\left[\beta_0+\beta_1 X(\mathbf{s})+\boldsymbol{\alpha}^T \mathbf{U}(\mathbf{s})\right] / \sqrt{1-\frac{T_{12}^2}{T_{11}}}\right) \tag{11} \end{aligned} $$ 其中，如 `式 (10)` 所示，$\beta_0 = \mu_2 − (T_{12}/T_{11})\mu_1$，并且 $\beta_1 = T_{12}/T_{11}$。 ## 4. 数据融合空间和时空数据的 `融合（Fusion）` 或 `同化（Assimilation）`（在本文中两个属于可以互换）可以以各种形式出现。有大量关于气候建模数据同化的文献（Kalnay 2003）^[30]以及越来越多的关于环境暴露评估数据融合的文献（本文重点）。从某种意义上说，所有这些工作都包含在计算机模型输出的数据同化的更大范畴之中。数据融合的总体目标是综合关于同一地区、同一时间段内的多个空间和时空响应数据源。例如：环境应用中对污染物（如 PM2.5）的建模，我们可能同时有站点监测数据、卫星数据和计算机模型结果，期望能够通过数据融合改进克里金法，并改善短期预报效果。数据层在空间和时间上几乎总是存在错位的情况，例如，一个数据是点参考的，而另一个数据则是面元单元的，或者一个是小时数据，而另一个日数据。为了方便说明，我们下面仅局限在两个数据层的情况：其中一个数据层来自监测站数据（即地统计数据源），而另一个数据层来自于按照某种空间分辨率的网格单元作出的环境模拟数据（不是实际数据）。 ### 4.1 积分方法 **（1）随机积分** 在 Poole & Raftery (2000 ^[40]) 的早期贝叶斯融合工作基础上，Fuentes & Raftery (2005 ^[22]) 提出了一种融合方法。它概念化了一个真实的监测数据表面，并查看监测站数据以及模型输出数据随真实表面以适当方式发生变化时的情况。特别地，网格单元 $A$ 中的平均值（由 $Z(A)$ 表示）来自于对一个被称为 “支撑问题变化” 的空间过程实现的随机积分 (Gelfand 等， 2001 ^[26]) $$ Z(A) = \frac{1}{|A|} \int_A Z(\mathbf{s}) d \mathbf{s} \tag{12} $$ 其中 $|A|$ 表示网格单元 $A$ 的面积，$Z(A)$ 称为块平均。作为某个随机过程实现的一个积分，$Z(A)$ 是一个随机变量且永远无法显式使用。Fuentes 和 Raftery（2005 ^[22]）提出了一种使用 `式 (12)` 的融合建模方案： **（2）测站数据的建模（点参考数据）** 令 $Z(\mathbf{s})$ 和 $Y(\mathbf{s})$ 分别表示测站 $\mathbf{s}$ 处的观测数据和相关联的真实过程值。第一个模型假设是： $$ Z(\mathbf{s}) = Y(\mathbf{s}) + \varepsilon(\mathbf{s}) \tag{13} $$ 其中 $\varepsilon (\mathbf{s}) \sim N(0, \sigma_\varepsilon)$ 表示位置 $\mathbf{s}$ 处的测量误差。假设真实过程为 $$ Y(\mathbf{s})= \mu(\mathbf{s}) + \eta(\mathbf{s}) \tag{14} $$ 其中 $\mu(\mathbf{s})$ 提供了空间均值表面，通常以趋势、高程等为特征；误差项 $\eta(\mathbf{s})$ 为零均值的高斯过程。 **（3）测站数据的面元化** 计算机模型输出的面元尺度数据通常被假设为有偏的和/或未校准的。因此，该模型在概念上的点参考输出（表示为 $Q(\mathbf{s})$ ）被假定为与真实表面有关： $$ Q(\mathbf{s}) = a(\mathbf{s}) + b(\mathbf{s}) Y(\mathbf{s}) + \delta (\mathbf{s}) \tag{15} $$ 其中 $a(\mathbf{s})$ 和 $b(\mathbf{s})$ 分别表示局部加法和乘法偏差（参考 `第 3 节` 中的空间变系数模型）。假设误差项 $\delta (\mathbf{s})$ 是一个由 $\mathcal{N}(0, \sigma_\delta)$ 给出的白噪声过程。同时，由于计算机模型输出数据是在网格 $A_1, \ldots ,A_J$ 上提供的，所以通过对 `式(15)` 的随机积分，可以将点参考过程转换为网格数据，即 $Q(A_j) = \int_{A_j} a(\mathbf{s}) d \mathbf{s} + \int_{A_j} b(\mathbf{s})Y(\mathbf{s}) d \mathbf{s} + \int_{A_j} \delta(\mathbf{s}) d \mathbf{s}$。根据观察，不稳定的模型拟合（由于可识别性问题）主要来自于空间变系数 $b(\mathbf{s})$，因此令其为常数 $b(\mathbf{s}) = b$s。 **（4）融合预测** 这样，新位置 $\mathbf{s}_0$ 处的空间预测可以通过后验预测分布 $f(Y(\mathbf{s}_0)|\mathbf{Z},\mathbf{Q})$ 获得，其中 $\mathbf{Z}$ 表示所有测站数据，$\mathbf{Q}$ 表示所有网格级的计算机模型输出数据 $\mathbf{Q} = \{Q(A_1), \ldots , Q(A_J)\}$。当存在大量网格单元时，此融合策略在计算上变得不可行；每个随机积分都需要近似，并且引入的空间随机变量总数变得非常庞大（参见 `第 6 节` 中对大数据问题的讨论）。在多个时间段的动态实施就更加不可行了。事实上，虽然我们有非常多的网格单元，但监测站数量相对稀少。因此提议应该比例缩小而不是积分变大。 ### 4.2 降尺度方法降尺度器模型由 Berrocal 等（2010 ^[7] ）提出。为了便于解释，我们将其与臭氧应用关联。示例考虑了两个地面臭氧数据源（每日 $8$ 小时最大值，以 $ppb$ 为单元）。其中一个是地统计数据，来自国家空气监测网络；第二个数据源是 CMAQ，这是一个估算每日 $8$ 小时最大臭氧浓度的数值模型，CMAQ 模型的输出是以网格单元（ $12 \times 12$ 公里）为分辨率拼贴形成的表面。在美国东部，有 $861$ 个监测站和 $40,044$ 个网格单元，网格单元比监测站多得多；此时实施 Fuentes & Raftery (2005 ^[22]) 的融合模型需要做大量的块平均，使其在计算上不堪重负。对于降尺度器，我们令 $Z(\mathbf{s})$ 表示在 $\mathbf{s}$ 处观测到的臭氧浓度的平方根（服务于工作习惯）。我们使用 $X(B)$ 表示网格单元 $B$ 上的数值模型输出的平方根。每个点 $\mathbf{s}$ 都与其所在的 CMAQ 网格单元相关联。我们通过以下方式将观测数据与 CMAQ 模型输出相关联。对于网格 $B$ 中的每个 $\mathbf{s}$，我们假设 $$ Z(\mathbf{s}) = \tilde{\beta}_0(\mathbf{s}) + \tilde{\beta}_1(\mathbf{s}) X(B) + \varepsilon(\mathbf{s}), \varepsilon(\mathbf{s}) \sim \mathcal{N}(0,\tau^2) \tag{16} $$ 其中 $\tilde{\beta}_0(\mathbf{s}) = \beta_0 +\beta_0 (\mathbf{s})$，$\tilde{\beta}_1(\mathbf{s}) = \beta_1 +\beta_1 (\mathbf{s})$，$\varepsilon (\mathbf{s})$ 是块方差为 $\tau^2$ 的白噪声过程。 $\beta_0$ 和 $\beta_1$ 代表 CMAQ 模型的整体截距和斜率偏差，$\beta_0 (\mathbf{s})$ 和 $\beta_1 (\mathbf{s})$ 分别是对加法偏差和乘法偏差的局部调整。空间变系数 $\beta_0 (\mathbf{s})$ 和 $\beta_1 (\mathbf{s})$ 依次使用共区域化方法建模为双变量均值为零的高斯空间过程，如 `第 3 节` 所述。总而言之，$\tilde{\beta}_0(\mathbf{s}) + \tilde{\beta}_1(\mathbf{s}) X(B)$ 扮演了在上一小节中真实值 $Y(\mathbf{s})$ 的角色。此模型很简单，但可以提供局部级别的校准，并赋予空间过程 $Z(\mathbf{s})$ 灵活的非平稳协方差结构。它比贝叶斯融合更容易拟合，因为不需要处理随机积分。此外，为了拟合模型，我们只需要处理与监测站相关的响应，这与网格单元的数量相比数量要少得多。 ## 5. 时空地统计模型在过去几年中，由于空间和时间索引的数据集的激增，时空建模受到了越来越多的关注，并且随之而来的是需要了解它们。例如，在空气污染研究中，我们不仅对污染物表面的空间性质感兴趣，而且对这个表面如何随时间变化感兴趣。当时间被离散化为通常的整数间隔时，我们可以通过两种方式查看时空索引的地统计数据 $Y(\mathbf{s}, t)$。写作 $Y(\mathbf{s}, t) = Y_s(t)$，我们有一个空间变化的时间序列模型。写作 $Y(\mathbf{s}, t) = Y_t(\mathbf{s})$，我们有一个随时间变化的空间模型。例如，对于气候数据，我们可能对特定时间的温度或降水的空间模式感兴趣，但也对天气的动态模式感兴趣。对于房地产市场，我们可能对单户住宅销售市场每季度或每年的变化情况感兴趣。在这种情况下，我们不会在每个时间段观测到相同的位置；数据是横截面的，而不是纵向的。在时间序列的情况下，我们可以将数据视为每个位置的多元测量向量，然后像上一节一样尝试使用多元空间数据模型。对于短系列，这可能是合理的；对于更长的序列，我们可能想要介绍通常时间序列建模的各个方面。缺失数据问题在时空环境中变得更加复杂。现在我们可能会遇到缺少空间信息的时间点、缺少特定（可能是未来）时间点信息的位置，或者它们的组合。在这里，贝叶斯分层方法特别有用，因为它不仅有助于组织我们对模型的思考，而且还可以适当地为离观测数据（空间或时间）较远的预测附加不确定性。 ### 5.1 连续时间的地统计建模具有连续时间的地统计数据建模自然扩展了第 2 节的建模。让 $Y(\mathbf{s}, t)$ 表示时间 $t$ 位置 $\mathbf{s}$ 的测量值。对于假设大致呈高斯分布的连续数据，我们可以写出一般形式 $$ Y(\mathbf{s},t) = \mu(\mathbf{s},t) + e(\mathbf{s},t) \tag{17} $$ 其中 $\m\mu(\mathbf{s}, t)$ 表示均值结构，$e(\mathbf{s}, t)$ 表示残差。如果 $\mathbf{x}(\mathbf{s}, t)$ 是与 $Y(\mathbf{s}, t)$ 相关联的协变量向量，那么我们可以设置 $\mu(\mathbf{s}, t) = \mathbf{x}(\mathbf{s}, t)^T \boldsymbol{\beta}(\mathbf{s}, t)$。请注意，这种形式允许时空变化系数（参见 `第 3 节`）。在很多应用中，经常采用 $\boldsymbol{\beta}(\mathbf{s}, t) = \boldsymbol{\beta}$ 的时空固定系数，但设置 $\boldsymbol{\beta}(\mathbf{s}, t) = \boldsymbol{\beta}(\mathbf{s})$ 可以产生空间变化的系数；设置 $\boldsymbol{\beta}(\mathbf{s}, t) = \boldsymbol{\beta}(t)$ 的情况在文献中不存在，除非 $t$ 被离散化为下面 `第 5.3 节` 的情况。最后，$e(\mathbf{s}, t)$ 通常会重写为 $w(\mathbf{s}, t) + \varepsilon (\mathbf{s}, t)$，其中 $\varepsilon (\mathbf{s}, t)$ 是高斯白噪声过程，$w(\mathbf{s}, t)$ 是均值为零的时空过程。我们可以将 `式 (17)` 视为一个在给定 $\{\mu(\mathbf{s}, t)\}$ 和 $\{w(\mathbf{s}, t)\}$ 时，第一阶段具有条件独立性的分层模型。然后，本着 `第 2.5 节` 的精神，我们可以将高斯第一阶段模型替换为另一种第一阶段模型（例如，指数族模型），并写成 $Y(\mathbf{s}, t) \sim f(y(\mathbf{s}, t) | \mu(\mathbf{s}, t), w(\mathbf{s}, t))$, 其中 $$ f(y(\mathbf{s}, t) \mid \mu(\mathbf{s}, t), w(\mathbf{s}, t))=b(y(\mathbf{s}, t)) \exp \{\gamma[\eta(\mathbf{s}, t) y(\mathbf{s}, t)-\chi(\eta(\mathbf{s}, t))]\} \tag{18} $$ 式中 $γ$ 是正值的离散参数，$g(η(\mathbf{s}, t)) = \mu(\mathbf{s}, t) + w(\mathbf{s}, t)$ 为某种链接函数 $g$。因此，挑战在于指定时空过程 $w(\mathbf{s}, t)$。简单版本将指定 $w(\mathbf{s}, t) = w(\mathbf{s}) + α(t)$（加法）或 $w(\mathbf{s}, t) = w(\mathbf{s})α(t)$（乘法），其中 $w(\mathbf{s})$ 和 $α(t)$ 分别是空间和时间上的独立均值零高斯过程。更常见的是在 $\mathbb{R}^2 \times \mathb{R}^{+}$ 上提供时空过程。仅限于高斯过程，我们只需要指定一个有效的时空协方差函数。在这里，意味着对于任何一组位置和任何一组时间点，随机变量的结果集的协方差矩阵是正定的。这里重要的一点是将 $(\mathbf{s}, t)$ 视为三维向量并在 $\mathbf{R}^3$ 上提出有效的相关函数是不明智的。因为无论缩放比例如何，空间距离与时间尺度上的距离无关。因此，假定固定时空协方差规范采用以下形式：$\operatorname{cov}(w(\mathbf{s}, t),w(\mathbf{s}^\prime, t^\prime)) = c(\mathbf{s}−\mathbf{s}^\prime , t−t^\prime)$。 [各向同性形式设置为 $\operatorname{cov}(w(\mathbf{s}, t),w(\mathbf{s}^\prime, t^\prime)) = c(\|s−s ^\prime\|, |t− t^\prime|)$]。一个经常使用的选择是采用可分离形式： $$ \operatorname{cov}(w(\mathbf{s}, t),w(\mathbf{s}^\prime, t^\prime)) = \sigma^2 \rho^{(1)}(\mathbf{s}−\mathbf{s}^\prime;\phi) \rho^{(2)}(t−t^\prime; \psi) \tag{19} $$ 其中 $\rho^{(1)}$ 是有效的二维相关函数，$\rho^{(2)}$ 是有效的一维相关函数（参见 Mardia & Goodall 1993 ^[36] 及其中的参考文献）。 `式 (19)` 表明依赖性在空间和时间上以乘法方式衰减。为什么 `式 (19)` 有效？对于位置 $\mathbf{s}_1, \ldots ,\mathbf{s}_I$ 和时间 $t_1, \ldots ,t_J$，将变量收集到一个向量 $\mathbf{w}_{\mathbf{s}}^T = (\mathbf{w}^T(\mathbf{s}_1), \ldots, \mathbf{w}^T(\mathbf{s}_I))$ 中，其中 $\mathbf{w}(\mathbf{s}_i) = (w(\mathbf{s}_i,t_1), \ldots ,w(\mathbf{s}_i,t)J))^T$ , $\mathbf{w}_s$ 的协方差矩阵为： $$ \sum_{\mathbf{w}_s}\left(\sigma^2, \boldsymbol{\phi}, \boldsymbol{\psi}\right)=\sigma^2 H_s(\boldsymbol{\phi}) \otimes H_t(\boldsymbol{\psi}) \tag{20} $$ 其中 $\otimes$ 再次表示 Kronecker 积。在 `式 (20)` 中，$H_s(\phi)$ 为 $I \times I$ 矩阵，其元素为 $H_s(\psi))_{ii^\prime} = \rho^{(1)}(t_i− t_{i^\prime};\phi)$； $H_t(\psi)$ 是 $J \times J$ 矩阵，其元素为 $H_t(\psi))_{jj^\prime} = \rho^{(2)}(t_j− t_{j^\prime};\psi)$。 `式 (20)` 表明 $\Sigma_{\mathbf{w}_s}$ 是正定的。因此，$\mathbf{w}_s$ 将遵循 $IJ$ 维多元高斯分布，其均值向量为 $\boldsymbol{\mu}_s(\boldsymbol{\beta})$ 和协方差矩阵来自于 `式（20）`。基于相关数据 $\mathbf{Y}_s^T=(\mathbf{Y}^T(\mathbf{s}_1), \ldots,\mathbf{Y}^T(\mathbf{s}_I))$，在给定 $\boldsymbol{\beta}$、$\sigma^2$、$\phi$ 和 $\psi$ 先验时，贝叶斯模型就完全确定了。可以采用类似于静态空间情况中基于仿真的模型拟合方法，不过需要注意 $Y_s$ 引起的对数似然 $$ -\frac{1}{2} \log \left|\sigma^2 H_s(\boldsymbol{\phi}) \otimes H_t(\boldsymbol{\psi})\right|-\frac{1}{2 \sigma^2}\left(\boldsymbol{Y}_s-\boldsymbol{\mu}_s(\boldsymbol{\beta})\right)^T\left(H_s(\boldsymbol{\phi}) \otimes H_t(\boldsymbol{\psi})\right)^{-1}\left(\boldsymbol{Y}_s-\boldsymbol{\mu}_s(\boldsymbol{\beta})\right) $$ 然而，$|\sigma^2 H_s(\phi) \otimes H_t(\psi)| = (\sigma^2)^{IJ} |H_s(\phi)|^J |H_t(\psi)|^I$ 和克罗内克积的性质。换句话说，即使 `式 (20)` 是 $IJ \times IJ$，我们只需要 $I \times I$ 和 $J \times J$ 矩阵的行列式和逆矩阵，从而加快似然评估和吉布斯抽样。预测很简单；参见 Banerjee 等 (2014 ^[2]) 了解详情。 ### 5.2 不可分离的时空模型在许多应用程序中，我们期望在捕获数据中的依赖结构方面进行时空交互。也就是说，我们期望空间依赖性随时间变化，或者，我们期望时间依赖性随空间位置变化。例子包括大气污染物沉积的时空趋势、土壤水分含量空间格局的时间演变、疾病和污染物暴露的时空格局，以及生态学中种群动态的空间特征。在这方面， `式 (19)` 中时空协方差函数的可分离形式虽然便于计算且具有吸引人的解释，但限制了时空相互作用的性质。没有时空依赖性。 Cressie & Huang (1999 ^[11] ) 介绍了一类灵活的不可分平稳协方差函数，允许时空相互作用。然而，它们在频谱域中工作，并且要求可以显式计算 $c(\mathbf{s} − \mathbf{s}^\prime, t− t^\prime)$。也就是说，傅里叶逆变换可以得到封闭形式，这只发生在非常特殊的情况下。 Gneiting (2002 ^[27]) 的工作采用了类似的方法，但获得了非常普遍的有效时空模型类别，这些模型不依赖于封闭形式的逆傅立叶。一个简单的例子是类 $c(\mathbf{s} − \mathbf{s}^\prime, t− t^\prime) = \sigma^2 (|t− t^\prime| + 1)^{−1} \exp(− \| \mathbf{s} − \mathbf{s}^\prime \| (|t− t^\prime| + 1)^{−\boldsymbol{\beta}/2})$。其中，$\boldsymbol{\beta}$ 为时空相互作用参数； $\boldsymbol{\beta}= 0$ 提供了可分离的形式。 Stein (2005 ^[49] ) 也在谱域中工作，提供了一类谱密度，其产生的时空协方差函数与灵活的分析行为不可分离。特别地，谱密度是 $$ \hat{c}(\mathbf{w}, v) \propto\left[c_1\left(\alpha_1^2+\|\mathbf{w}\|^2\right)^{\alpha_1}+c_2\left(\alpha_2+v^2\right)^{\alpha_2}\right]^{-v} $$ 不幸的是，无法显式计算相关的协方差函数；快速傅里叶变换提供了最好的计算前景。与 Gneiting 的类不同，可分离性不会作为特例或极限情况出现。如需进一步讨论不可分离的规范，请读者参阅 Gelfand 等（2010 ^[23]）的第 23 章。 ### 5.3 动态时空模型现在，我们转向离散时间。在这里，动态线性模型（通常称为状态空间模型）提供了一个通用框架来拟合一系列时间结构模型（West & Harrison 1997 ^[56] ）。我们简要回顾一下一般的动态线性建模框架。设 $\mathbf{Y}_t$ 为时间 $t$ 的 $m \times 1$ 可观测向量。 $\mathbf{Y}_t$ 通过测量方程与 $p \times 1$ 向量 $θ_t$ 相关，称为状态向量。通常，$θ_t$ 的元素是不可观测的，而是由一阶马尔可夫过程产生的，从而产生一个转移方程。因此，我们可以将上述框架描述为 $$ \boldsymbol{Y}_t=F_t \boldsymbol{\theta}_t+\boldsymbol{\varepsilon}_t, \quad \varepsilon_t \sim N\left(\mathbf{0}, \sum_t^{\varepsilon}\right) $$ $$ \boldsymbol{\theta}_t=G_t \boldsymbol{\theta}_{t-1}+\boldsymbol{\eta}_t, \boldsymbol{\eta}_t \sim N\left(\mathbf{0}, \sum_t^\eta\right) $$ 其中 $F_t$ 和 $G_t$ 分别是 $m \times p$ 和 $p \times p$ 矩阵。在测量方程中，$\varepsilon_t$ 是一个 $m \times 1$ 向量，由均值为 $0$ 的连续不相关高斯变量和一个 $m \times m$ 协方差矩阵组成，在转移方程中，$η$ 是序列不相关的零中心高斯扰动的 $p\times 1$ 向量和相应的 $p\times p$ 协方差矩阵。关联结构可以跨时间显式计算，例如 $\operatorname{cov}\left(\boldsymbol{\theta}_t, \boldsymbol{\theta}_{t-1}\right)=G_t \operatorname{var}\left(\boldsymbol{\theta}_{t-1}\right)$ and $\operatorname{cov}\left(\boldsymbol{Y}_t, \boldsymbol{Y}_{t-1}\right)=F_t G_t \operatorname{var}\left(\boldsymbol{\theta}_{t-1}\right) F_t^T$。 $F_t$ 和 $G_t$ 称为系统矩阵，允许随时间变化。它们可能涉及未知参数，但给定参数后，时间演化会以预定方式出现。矩阵 $F_t$ 通常由手头问题的设计指定，而 $G_t$ 通过建模假设指定。例如，$G_t= Ip$，$p\times p$ 单元矩阵，为 $θ_t$ 提供随机游走。无论如何，系统是线性的，对于任何时间点 $t$，$\mathbf{Y}_t$ 可以表示为现在的 $\varepsilon_t$ 与现在和过去的 $η_t$ 的线性组合。时空模型的形式化考虑一组测站 $S= \{\mathbf{s}_1, \ldots ,\mathbf{s}_{N_s}\}$ 和时间点 $T= {t_1, \ldots ,t_{N_t}\}$，产生观测值 $Y(\mathbf{s}, t)$ 和协变量向量 $\mathbf{x}(\mathbf{s}, t)$，对于每个 $(\mathbf{s}, t) \in S\times T$。响应 $Y(\mathbf{s}, t)$ 首先通过测量方程建模，该方程包含测量误差 $\varepsilon(\mathbf{s}, t)$，作为序列和空间不相关的零中心高斯扰动。转移方程现在涉及协变量的回归参数（斜率）。斜率向量，比如 $\tilde{\beta}(\mathbf{s}, t)$，被分解为纯时间分量 $\boldsymbol{\beta}_t$ 和时空分量 $\boldsymbol{\beta}(\mathbf{s}, t)$。这两者都是通过转移方程生成的，及时捕获它们的马尔可夫依赖性。纯时间分量的转移方程与通常的状态空间建模一样，而时空分量是由多元高斯空间过程生成的。因此，我们可以将时空建模框架写成 $$ Y(\mathbf{s}, t)=\mu(\mathbf{s}, t)+\varepsilon(\mathbf{s}, t) ; \varepsilon(\mathbf{s}, t) \stackrel{\text { ind }}{\sim} N\left(0, \sigma^2\right), $$ $$ \mu(\mathbf{s}, t)=\mathbf{x}^T(\mathbf{s}, t) \tilde{\boldsymbol{\beta}}(\mathbf{s}, t) $$ $$ \tilde{\boldsymbol{\beta}}(\mathbf{s}, t)=\boldsymbol{\beta}_t+\boldsymbol{\beta}(\mathbf{s}, t) $$ $$ \beta_t=\beta_{t-1}+\boldsymbol{\eta}_t, \quad \boldsymbol{\eta}_t \stackrel{\text { ind }}{\sim} N_p\left(\mathbf{0}, \sum_\eta\right) $$ $$ \boldsymbol{\beta}(\mathbf{s}, t)=\boldsymbol{\beta}(\mathbf{s}, t-1)+\mathbf{w}(\mathbf{s}, t) $$ 我们引入了 $\mathbf{w}(\mathbf{s}, t)$ 的共区域化线性模型（如第 3 节中所述），即 $\mathbf{w}(\mathbf{s}, t) = A \mathbf{v} (\mathbf{s}, t)$，其中 $\mathbf{v}(\mathbf{s}, t) = (v_1(\mathbf{s} , t), \ldots ,v_p(\mathbf{s}, t))^T$，产生 $\Sigma_w= AA^T$。 $v_l(\mathbf{s}, t)$ 是具有单位方差和相关函数 $\rho_l(·; \phi_l)$ 的高斯过程的连续独立复制，此后记为 $\operatorname{GP}(0, \rho_l(·; \phi_l))$，对于 $l= 1， \ldots ,p$ 并且在 $l$ 上独立。允许空间变化的系数 $\boldsymbol{\beta}(\mathbf{s}, t)$ 与 $\mathbf{x}(\mathbf{s}, t)$ 相关联，为 $X_s$ 关于 $Y_s$ 提供了任意丰富的解释关系。我们可以根据上述规范计算 Y 的一般关联结构。例如，$cov (Y(\mathbf{s}, t), Y(\mathbf{s}^\prime, t− 1)) = \mathbf{x}^T(\mathbf{s}, t) \Sigma_{\tilde{\beta}(\mathbf{s},t),\tilde{\beta}(\mathbf{s}^\prime,t−1)} \mathbf{x}(\mathbf{s}, t − 1)$, 其中 $$ \sum_{\tilde{\boldsymbol{\beta}}(\mathbf{s}, t), \tilde{\boldsymbol{\beta}}\left(\mathbf{s}^{\prime}, t-1\right)}=(t-1)\left(\sum_\eta+\sum_{l=1}^p \rho_l\left(\mathbf{s}-\mathbf{s}^{\prime} ; \phi_l\right) \mathbf{a}_l \mathbf{a}_l^T\right) $$ 且 $$ \operatorname{var}(Y(\mathbf{s}, t))=\mathbf{x}^T(\mathbf{s}, t) t\left[\sum_\eta+A A^T\right] \mathbf{x}(\mathbf{s}, t) $$ 贝叶斯层次模型可以用先验指定完成，例如 $$ \beta_0 \sim N\left(\mathbf{m}_0, C_0\right) \text { and } \boldsymbol{\beta}(\cdot, 0) \equiv 0, $$ $$ \sum_\eta \sim I W\left(a_\eta, B_\eta\right), \quad \sum_w \sim I W\left(a_{\mathbf{w}}, B_{\mathbf{w}}\right) \text { and } \sigma_{\varepsilon}^2 \sim I G\left(a_{\varepsilon}, b_{\varepsilon}\right) $$ $$ \mathbf{m}_0 \sim N\left(\mathbf{0}, \sum_0\right) ; \sum_0=10^5 \times I_p $$ 其中 $B_η$ 和 $B_w$ 是逆 Wishart 分布的 $p \times p$ 精度（超参数）矩阵。用于拟合此类层次模型的 Gibbs 采样器相对简单（有关详细信息，请参见 Banerjee 等，2014 ^[2]）。 ## 6. 大数据问题毫不奇怪，随着计算能力的提高，我们正在收集越来越大的空间和时空数据集。我们现在着眼于以非常高的时间分辨率监测大型空间网络的环境污染物。我们现在以相对较高的分辨率在全球范围内研究海面温度。数据集现在达到 $10^6$ 或更高阶。模型拟合变得越来越具有挑战性。特别是，实施 Gibbs 采样或其他 MCMC 模型拟合算法需要重复评估各种完整条件密度函数。在使用高斯过程从随机效应构建的分层模型的情况下，这需要重复评估在高斯过程下出现的似然和/或联合密度和条件密度。对于较大的 $n$，与生成的 $n\times n$ 矩阵相关的计算可能不稳定，并且重复计算（针对基于采样的模型拟合）可能非常缓慢，甚至不可行。矩阵求逆的 Cholesky 分解需要大约 $O(n^3/3)$ flops。这种情况被非正式地称为 “大 $n$ 问题”。对于多变量模型，例如，在某个位置进行 $p$ 个测量会导致 $np\times np$ 的矩阵。对于时空模型（例如，在 $T$ 时间点的空间时间序列）导致 $nT\times nT$ 的矩阵。通过建模策略也许可以简化成 $Tn\times n$ 矩阵或一个 $n\times n$ 和一个 $T\times T$ 矩阵，但如果 $n$ 很大，问题仍然存在。 Banerjee 等（2014 ^[2] ，第 12 章）对大 $n$ 问题进行了讨论。因为我们的地统计建模下的挑战是处理空间过程 $w(\mathbf{s})$ ，所以我们考虑的重点是 “使用计算上更可行的高斯过程替换常规高斯过程” 的策略： - 一种方法是通过 `降维` 来形成预测过程 (Banerjee 等， 2008 ^[3] , Finley 等， 2009 ^[21] )，该方法在 Sang & Huang (2012 ^[43] ) 和 Katzfuss (2016 ^[32] ) 中得到了进一步扩展； - 第二种方法是通过似然来引入稀疏性，并将其扩展到有效的空间过程（Datta 等， 2016a ^[14] ,b ^[16] ）。 ### 6.1 降秩模型 **（1）降秩的基本思想** 处理大型空间数据集的一种常见方法是设计能够实现降维的模型（例如，参见 Wikle 2010 ^[57] 和其中的参考资料，对分层模型中的降维进行了很好的回顾）。其基本想法是将空间过程 $w(\mathbf{s})$ 替换为一个降秩的近似空间过程 $\tilde{w}(\mathbf{s})$。新的降秩（或低秩）模型通常建立在某种基于较小位置集实现的隐过程实现之上。准确地说，令新过程在 $\mathbf{s}$ 点的实现为： $$ \tilde{w}(\mathbf{s})=\sum_{j=1}^m l\left(\mathbf{s}, \mathbf{s}_j^*\right) Z\left(\mathbf{s}_j^*\right) \tag{31} $$ 其中 $Z(\mathbf{s})$ 是一个理想中的完整过程。空间过程的某次实现 $\tilde{w}(\mathbf{s})$ 完全由函数 $l(\cdot, \cdot)$ 和 $m$ 个位置组成的变量集合 $\{Z(\mathbf{s}_j^*)$，$j=1,2, \ldots, m\}$ 决定，$s_j^*$ 的集合被称为 “结（Knots）”。对于原始数据集中的 $n$ 个位置，与其关联的新过程可以表示成向量 $\tilde{\mathbf{w}}=(\tilde{w}(\mathbf{s}_1), \tilde{w}(\mathbf{s}_2), \ldots, \tilde{w}(\mathbf{s}_n))^T$, 则有如下的向量化表示形式： $$ \tilde{\mathbf{w}}=\mathbf{L} \mathbf{z}^* \tag{32} $$ 其中 $\mathbf{L}$ 是 $n \times m$ 的矩阵（ $m m$，$\tilde{\mathbf{w}}$ 的联合分布就是奇异的。但无论如何，我们确实创建了一个有效的随机过程。特别地，协方差函数可以改写为： $$ \operatorname{cov}\left(\tilde{w}(\mathbf{s}), \tilde{w}\left(\mathbf{s}^{\prime}\right)\right)=\mathbf{l}(\mathbf{s})^T \sum_{\mathbf{z}^*} \mathbf{l}\left(\mathbf{s}^{\prime}\right) \tag{33} $$ 其中 $\mathbf{l}(\mathbf{s})$ 是元素为 $l(\mathbf{s}_i, \mathbf{s}_j^*)$ 的 $m \times 1$ 向量。从 `式（33）` 可以看到，即使 $Z(\mathbf{s})$ 是平稳的，引入的协方差函数也可能是非平稳的。但其中存在一个特例：如果 $\mathbf{z}^*$ 是高斯分布，则 $\tilde{w}(s)$ 是一个高斯过程。 **（2）核卷积** 如果 $Z$ 是独立同分布的，且服从均值为 $0$ 、方差为 $\sigma^2$ 的高斯分布，即 $Z(\mathbf{s})$ 是一个白噪声过程，则 `式 33` 的协方差函数可以简化为 $\sigma^2 \mathbf{l}(\mathbf{s})^T \mathbf{l}(\mathbf{s}^{\prime})$。这种形式曾出现在 Barry & Ver Hoef (1996 ^[4]) 和 Higdon (1998 ^[28] ) 中，前者称之为 `移动平均模型`，后者称之为 `核卷积`。特别地，$l$ 的一个自然选择是核函数，如 $K(\mathbf{s}-\mathbf{s}^{\prime})$，它给更靠近 $\mathbf{s}$ 的 $\mathbf{s}^{\prime}$ 赋予更多权重。核具有参数（这导致参数化的协方差函数）并且可能在空间上发生变化（Higdon 2002 ^[29]，Paciorek & Schervish 2006 ^[38] ）。降秩形式可以看作是形式如 $\tilde{w}(\mathbf{s})=\int_{R^2} K(\mathbf{s}-\mathbf {s}^{\prime}) Z(\mathbf{s}^{\prime}) \mathrm{d} \mathbf{s}^{\prime}$ 的空间过程离散化(见 Xia & Gelfand 2006 ^[58]年关于此离散近似的讨论）。高斯核是经常被使用的，尽管其会导致高斯协方差函数，并且产生的空间过程过于平滑以至于在实践中不能令人满意（Stein 1999a ^[47]，Paciorek ＆Schervish 2006 ^[38] ）。核卷积所能得到的随机过程范围具有局限性；例如，广泛使用的指数协方差函数就并非来自核卷积。 **（3）克里金预测过程** 与之相比，Banerjee 等 (2008 ^[3] ) 受克里金思想启发提出了另外一类模型。考虑一组结 $\mathscr{S}^*=\left\{\mathbf{s}_1^*, \ldots, \mathbf{s}_m^*\right\}$，这可能但不一定是完整观测位置集 $\mathscr{S}=\left\{\mathbf{s}_1, \mathbf{s}_2, \ldots, \mathbf{s}_n\right\}$ 的一个子集。假设 $w(\mathbf{s}) \sim \operatorname{GP}(0, C(\cdot ; \boldsymbol{\theta}))$ ， $\mathbf{w}^*$ 是 $w(\mathbf{s})$ 在 $\mathscr{S}^*$ 上的一次实现。也就是说，$\mathbf{w}^*$ 是元素为 $w\left(\mathbf{s}_i^*\right)$ 的 $m \times 1$ 向量，并且服从多元高斯分布 $\mathbf{w}^* \sim \operatorname{MVN}(\mathbf{0}, C^*(\boldsymbol{\theta}))$，其中 $C^*(\boldsymbol{\theta})$ 是 $m \times m$ 协方差矩阵，其元素为 $C(\mathbf{s}_i^*, \mathbf{s}_j^* ; \boldsymbol{\theta})$。则测站 $\mathbf{s}_0$ 的空间插值（克里金法）由下式给出： $$ \tilde{w}\left(\mathbf{s}_0\right)=E\left[w\left(\mathbf{s}_0\right) \mid \mathbf{w}^*\right]=\mathbf{C}^T\left(\mathbf{s}_0 ; \boldsymbol{\theta}\right) C^{*-1}(\boldsymbol{\theta}) \mathbf{w}^* \tag{34} $$ 其中 $\mathbf{C}\left(\mathbf{s}_0 ; \boldsymbol{\theta}\right)$ 是元素为 $C\left(\mathbf{s}_0, \mathbf {s}_j^* ; \boldsymbol{\theta}\right)$ 的 $m \times 1$ 向量。这个单测站插值模型定义了一个具有如下协方差函数的空间过程 $\tilde{w}(\mathbf{s}) \sim \operatorname{GP}(0, \tilde{C}(\cdot))$ $$ \tilde{C}\left(\mathbf{s}, \mathbf{s}^{\prime} ; \boldsymbol{\theta}\right)=\mathbf{C}^T(\mathbf{s} ; \boldsymbol{\theta}) C^{*-1}(\boldsymbol{\theta}) \mathbf{C}\left(\mathbf{s}^{\prime}, \boldsymbol{\theta}\right) \tag{35} $$ 其中 $\mathbf{C}(\mathbf{s} ; \boldsymbol{\theta})=\left[C\left(\mathbf{s}, \mathbf{s}_j^* ; \boldsymbol{\theta} \right)\right]_{j=1}^m$。 $\tilde{w}(\mathbf{s})$ 被称为从父过程 $w(\mathbf{s})$ 派生的预测过程。$\tilde{w}(\mathbf{s})$ 正是基于 $w(\mathbf{s})$ 在 $\mathscr{S}^*$ 上实现的克里金预测。给定父过程的协方差函数和 $\mathscr{S}^*$ ，则该过程是完全指定的。从 `式 35` 可以看出，无论 $w(\mathbf{s})$ 是否平稳，此过程都是非平稳的。与 `式 31` 的联系也很清楚：我们将 $Z(\mathbf{s})$ 视为父进程并且 $\mathbf{l}(\mathbf{s})=\mathbf{C}(\mathbf{s } ; \theta)^T C^{*-1}(\theta)$。因此，多变量的每一个空间过程都可以引入一个预测过程模型（实际上，根据节的数量，可以是任意多个）。后者将前者的过程实现建模至低维子空间，从而减少了计算负担。回顾 `式 1`，我们将 $w(\mathbf{s})$ 替换为 $\tilde{w}(\mathbf{s})$ 来获得预测过程模型： $$ Y(\mathbf{s})=\mathbf{x}^T(\mathbf{s}) \beta + \tilde{w}(\mathbf{s}) + \varepsilon(\mathbf{s}) \tag{36} $$ 因为 $\tilde{w}(\mathbf{s})=\mathbf{C}^T(\mathbf{s}) C^{*-1}(\theta) \mathbf{w}^*$, $\tilde{w}(\mathbf{s})$ 是 $\mathbf{w}^*$ 的一个空间变化线性变换。可以立即看到降维的作用。在拟合 `式 36` 中的模型时，$n$ 个随机效应 $\{w(\mathbf{s}_i), i=1,2, \ldots, n \}$ 被替换为在 $\mathbf{w}^*$ 中的 $m$ 个随机效应；我们可以使用仅涉及 $m \times m$ 矩阵的 $m$ 维联合分布。尽管高斯预测过程是由高斯过程引入的，但它本身就是空间随机效应的一个过程模型。此外，它是 Kullback-Leibler 意义上的最优投影过程。在给定 $\mathscr{S}^*$ 上的 $w(\mathbf{s}^*)$ 时，$\mathscr{S}$ 上的 $w(\mathbf{s})$ 的所有有限维后验中，该预测过程能够最小化逆 Kullback-Leibler 距离 (Csató 2002 ^[13], Banerjee 等， 2008 ^[3] )。不过，该预测过程的结果通常过于平滑。如果 $C(\mathbf{s} ; \theta)$ 在 $\|\mathbf{s}\|>0$ 时无限可微（实践中大多数协方差函数都满足的一个假设），则其实现结果也几乎肯定是无限可微的。 **（4）预测过程模型偏差** 预测过程在 $\mathscr{S}^*$ 上做确定性地插值，即 $\tilde{w}\left(\mathbf{s}_j^*\right)=w\left(\mathbf{s}_j ^*\right)$。在另一个极端，假设 $\mathscr{S}$ 和 $\mathscr{S}^*$ 不相交，则 $\mathbf{w} \mid \mathbf{w}^* \sim \operatorname{MVN}(c^T C^{*-1} \mathbf{w}^*, C-c^T C^{*-1} c)$ 其中 $c$ 是 $m \times n$ 矩阵，其列为 $\mathbf{C}(\mathbf{s}_i)$ ，而 $C$ 是 $\mathbf{w}$ 的一个 $n \times n$ 协方差矩阵。我们可以写成 $\mathbf{w}=\tilde{\mathbf{w}} + (\mathbf{w}-\tilde{\mathbf{w}})$，其中 $\mathbf{w}^*$ 的选择确定了右侧第二项的（条件）可变性，即 $\Sigma_{\mathbf{w}}$ 与 $\Sigma_{\tilde{w}}$ 的接近程度。它还表明，由于 $n$ 个变量由 $m$ 个随机变量确定（$m[21] ）提出了的一种简单的补救措施，即使用经过修改的预测过程： $$ \tilde{w}_{\varepsilon}(\mathbf{s})=\tilde{w}(\mathbf{s})+\tilde{\varepsilon}(\mathbf{s}) \tag{39} $$ 其中 $\tilde{\varepsilon}(\mathbf{s}) \stackrel{\text {iid }}{\sim} N\left(0, \delta^2(\mathbf{s})\right)$ ,并且 $\tilde{\mathcal{E}}(\mathbf{s})$ 独立于 $\tilde{w}(\mathbf{s})$。现在，$\tilde{w}_{\mathcal{E}}(\mathbf{s})$ 的方差与父进程 $w(\mathbf{s})$ 相等该补救措施在计算上是高效的，添加一个独立的空间变化块金不会产生大量的计算费用。 Sang 等 (2011 ^[44] ) 和 Sang & Huang (2012 ^[43] ) 指出，这种补救措施可能无法捕获小规模相关性，即非对角线项，因此他们建议对 $\tilde{\varepsilon}(\mathbf{s})$ 进行逐渐变窄的协方差调整。 Katzfuss (2013 ^[31] ) 在考虑预测过程建模时，扩展了 Sang 等 (2011 ^[44] ) 的逐渐变窄想法。他使用两个项来指定空间误差，结合低秩分量来捕获中长期依赖，并结合锥形残差分量来捕获局部依赖。非平稳 Matérn 协方差函数被用来提供更高的灵活性。与逐步优化残差过程的方法不同，Katzfuss (2016 ^[32]) 使用了一系列嵌套的预测过程，用于开发面向海量空间数据集的多分辨率近似。 ### 6.2 最近邻高斯过程 Datta 等 (2016a ^[14] , b^[16]) 没有采用驻留在 $m$ 维空间中的随机过程方法，而是通过构建稀疏性归纳的空间过程来提供一种替代的高效计算策略，该过程也由给定的高斯过程归纳。特别是，他们展示了（使用了低维条件分布的）似然近似如何生成（具有稀疏精度矩阵的）高斯密度，这些矩阵对应于有限维空间过程实现的分布，因此被称为最近邻高斯过程 (NNGP)。当 $n$ 非常大时，低秩过程中的较小的 $m$ 将趋于过度平滑并且不能期望充分捕获空间过程的表面。在这方面，NNGP 明显优于低秩模型。此外，它具有大规模可扩展性，可以方便地嵌入到分层模型中，用于分析包含超过 $10^5$ 个空间（或时空）点数据。令 $w(\mathbf{s}) \sim \operatorname{GP}(0,C(·,·|θ))$ 表示具有协方差函数 $C(\mathbf{s}, \mathbf{s}^\prime )$ 的以零为中心的高斯过程。设 $\mathscr{R} = {\mathbf{s}_1, \mathbf{s}_2, \ldots ,\mathbf{s}_r}$ 是 $\mathscr(D)$ 中位置的固定有限集合，我们称之为`参考集(Reference Set)`。如果 $w_\mathscr{R} = (w(\mathbf{s}_1),w(\mathbf{s}_2), \ldots ,w(\mathbf{s}_r ))^T$，则 $w_\mathscr{R} \sim N(0,C_\mathscr{R}(θ))$，其中 $C_\mathscr{R}(θ)$ 是一个 $r\times r$ 的正定矩阵，以 $C(\mathbf{s}_i,\mathbf{s}_j)$ 作为其第 $(i, j)$ 个元素。 NNGP 是一种分两步指定大规模可扩展高斯过程的模型。首先，我们指定了一个计算效率非常高的多元高斯分布，它近似于参考集上的过程实现。随后，我们通过基于 $\mathscr{R}$ 中位置的克里金方程将其扩展到域上的过程。这两个步骤如下所述。过程实现 $w_\mathscr{R}$ 的概率密度是 $N(0,C_\mathscr{R}(θ))$。第一步寻求该密度的近似值。根据 Vecchia (1988 ^[52]) 和 Stein 等 (2004 ^[51] ) ，我们将 $\mathbf{w}$ 的联合密度写为$p(\mathbf{w})= \prod_{i=1}^r p(\mathbf{w}(\mathbf{s}_i)) \mid \mathbf{w}_{H(\\mathbf{s}_i)}$ ，其中 $H(\mathbf{s}_1)$ 是空集，$H(\mathbf{s}_i) = \{\mathbf{s}_1, \mathbf{s}_2, \ldots ,\mathbf{s}_{i−1}\}$ 和 $w_{H(\mathbf{s}_i)} = \{ \mathbf{w}(\mathbf{s}_j) \mid \mathbf{s}^j \in H(\mathbf{s}_i)\}$ 对于 $2 \leq i \leq r$。我们用 $w_{N(\mathbf{s}_i)}$ 替换条件集 $w_H(\mathbf{s}_i)$，其中 $N(\mathbf{s}_i)$ 是 $\mathbf{s}$ 从 $H(\mathbf{s}_i)$ 中点的一组最近邻居（根据欧几里德距离定义），只要 $H(\mathbf{s}_i)$ 有多于 $m$ 个点。Datta 等 (2016a ^[14]) 表明所得密度是一个高斯密度。因此， $$ N\left(\mathbf{w}_{\mathscr{R}} \mid \mathbf{0}, C_{\mathscr{R}}(\boldsymbol{\theta})\right)=\prod_{i=1}^r p\left(w\left(\mathbf{s}_i\right) \mid \mathbf{w}_{H\left(\mathbf{s}_i\right)}\right) \approx \prod_{i=1}^r p\left(w\left(\mathbf{s}_i\right) \mid \mathbf{w}_{N\left(\mathbf{s}_i\right)}\right)=N\left(\mathbf{w}_{\mathscr{R}} \mid \mathbf{0}, K(\boldsymbol{\theta})\right) \tag{40} $$ `式（40）` 中近似值背后的基本思想是将条件集从 $H(\mathbf{s}_i)$ 压缩到 $N(\mathbf{s}_i)$。这类似于构造一个系数高斯概率图模型或稀疏贝叶斯网络（参见如，Murphy 2012 ^[37]，第 10 章）。对于参考集中的每个 $\mathbf{s}_i$，它假定条件独立性 $p(w(\mathbf{s}_i) | \mathbf{w}_{H(\mathbf{s}_i)}) = p(w(\mathbf{s}_i) | \mathbf{w}_{N(\mathbf{s}_i)})$。这导致具有稀疏精度矩阵 $K^{−1}(θ)$ 的多元高斯分布。Datta 等 (2016a ^[14]) 以封闭形式导出了 $K^{−1}(θ)$ 的确切结构。它最多有 $rm(m+1)/2$ 个非零元素。较小的 $m$（邻居数量）值会产生更稀疏的 $K^{-1}(θ)$。重要的是，$\operatorname{det}(K(θ))$ 是 $r$ 个最多 $m\times m$ 阶的决定因素的乘积。同样，对于 $m \leq r$，在存储和浮点运算方面的计算节省是显而易见的。第二步将这个近似扩展到一个过程。为此，我们首先为任意测站 $\mathbf{s}$ 定义 $N(\mathbf{s})$ 为 $\mathscr{R}$ 中测站中 $\mathbf{s}$ 的最近邻集合。然后，对于任何 $\mathbf{s} \notin \mathscr{R}$，我们定义 $$ w(\mathbf{s})=\sum_{j=1}^r a_j(\mathbf{s}) w\left(\mathbf{s}_j\right)+\eta(\mathbf{s}), \quad \eta\left(\mathbf{s}_i\right) \stackrel{i n d}{\sim} N(0, f(\mathbf{s})) \tag{41} $$ 其中 $\mathbf{a}_j (\mathbf{s}) = 0$，其中 $\mathbf{s}_j \notin N(\mathbf{s})$。这意味着非零的 $\mathbf{a}_j (\mathbf{s})$ 是克里金 $w(\mathbf{s})$ 在其参考集邻居上的系数，$f(\mathbf{s})$ 是相应的克里金方差（即克里金误差的方差）。 NNGP 在最近邻克里金法方面的发展是在 Datta 等 (2016c ^[15]) 中发展起来的。对应于 NNGP 的协方差函数很容易推导，并在 Datta 等 (2016a ^[14]) 中明确提供。尽管结果 $K(θ)$ 对于我们定义的序列 $H(\mathbf{s}_1) \subseteq H(\mathbf{s}_2) \subseteq \ldots \subseteq H(\mathbf{s}_r)$ 并非是不变的，但 Vecchia (1988 ^[52]) 和 Stein 等 (2004 ^[51]) 断言 `式 (40)` 中的近似值对此排序不敏感。 Datta 等 (2016a ^[14]) 在文献补充资料的模拟实验也证实了这一点。 Datta 等 (2016a ^[14],b ^[16]) 的模拟实验似乎也表明，与 $r$ 相比，$m$ 通常可以被认为非常小,当 $m$ 大于 $15$ 时，似乎并没有表现出任何推断上的优势，即使对于包含超过 $10^5$ 个空间位置的数据集也是如此。为了使 `式 (40)` 中的近似能够有效并确保具有充分代表性的邻居集，参考集的大小 $r$ 需要足够大以表征空间域。然而，这并不妨碍涉及 NNGP 模型的计算，因为浮点运算的存储和数量在 $r$ 中始终是线性的。原则上，参考集 $\mathscr{R}$ 可以是研究域中任何有限的位置集。一个特别方便的选择是简单地将 $\mathscr{R}$ 视为数据集中观测到的位置的集合。Datta 等 (2016a ^[14]) 通过广泛的模拟实验和实际应用证明，这个简单选择似乎非常有效，并且 $\mathscr{R}$ 的替代选择似乎没有明显的好处。因为 NNGP 是一个合适的高斯过程，我们可以将其用于任何分层模型式中空间随机效应的先验。例如，如果 $\mathscr{R}$ 是一组观测到的位置，我们可以指定 $$ p(\boldsymbol{\theta}, \boldsymbol{\tau}, \beta) \times N\left(\mathbf{w}_{\mathscr{R}} \mid \mathbf{0}, K(\boldsymbol{\theta})\right) \times \prod_{i=1}^n N\left(y\left(\mathbf{s}_i\right) \mid \mathbf{x}^T\left(\mathbf{s}_i\right) \boldsymbol{\beta}+w\left(\mathbf{s}_i\right), \tau^2\right), $$ 其中 $y(\mathbf{s}_i)$ 是在位置 $\mathbf{s}_i$ 处测量的结果，$\mathbf{x} (\mathbf{s}_i)$ 是预测变量的向量，$\tau^2$ 是对应于 $\mathbf{y}(\mathbf{s}_i)$ 的块金方差。使用 Datta 等 (2016a ^[14]) 描述的含 Gibbs 步骤和随机游走 Metropolis 步骤的高效 MCMC 采样器，现在已经可以用于更新 $\mathbf{w}(\mathbf{s})$、$\boldsymbol{\beta}$、协方差参数 $\boldsymbol{\theta}$ 和块金 $\tau^2$。通过从后验预测分布中采样来执行任意位置 $\mathbf{s} \notin \mathscr{R}$ 的预测 $$ \int N\left(y(\mathbf{s}) \mid \mathbf{x}^T(\mathbf{s}) \boldsymbol{\beta}+w(\mathbf{s}), \tau^2\right) \times N\left(w(\mathbf{s}) \mid \mathbf{a}^T(\mathbf{s}) \mathbf{w}_{N(\mathbf{s})}, f(\mathbf{s})\right) \times p\left(\mathbf{w}_{\mathscr{R}}, \boldsymbol{\beta}, \tau^2, \boldsymbol{\theta} \mid \mathbf{y}_{\mathscr{R}}\right) $$ 其中 $x^T(\mathbf{s})$ 是位置 $\mathbf{s}$ 处的相应预测变量矩阵，$\mathbf{a}(\mathbf{s})$ 是 (41) 中非零的 $a_j (\mathbf{s})$ 组成的向量。建议读者参考 Datta 等 (2016a ^[14],b ^[16]) 了解详情。在本节中，我们对完全高斯过程（Full GP）、高斯预测过程（Gaussian Predictive Process）和近邻高斯过程（NNGP）进行了简单比较，使用了超过 $2,500$ 个空间位置的模拟高斯随机场。详细信息由 Datta 等 (2016a ^[14]) 提供。 `图 2(a-c)` 显示了来自完全高斯过程、$m = 10$ 的 NNGP 模型、$m=64$ 改进高斯预测过程模型的空间随机效应的后验中值估计。该图显示 $m= 10$ 的 NNGP 模型非常接近完全高斯过程模型，而 $m=64$ 的降秩改进高斯预测过程模型平滑了小尺度中的模式。最后一项观测突出了对降秩模型的主要批评之一（Stein 2014 ^[50]），并说明了为什么当真实表面具有精细空间分辨率细节时，这些模型通常会提供折衷的预测性能。我们将读者归纳至 Datta 等 (2016a ^[14]) 了解更多详情和说明。 ## 7. 其他挑战 ### 7.1 高效且用户友好的计算 MCMC 算法的计算发展（参见，例如，Robert & Casella 2004 ^[42]）极大地促进了贝叶斯分层模型在广泛学科中的普及，空间建模也不例外。然而，地统计模型拟合和推断的自动化实施面临巨大的挑战。 - 首先，需要昂贵的矩阵计算，这对于大型数据集来说可能会变得令人望而却步。 - 其次，非边缘化模型的拟合程序不太适合使用 BUGS 范式中的 Gibbs 采样器直接更新（Finley 等，2007 ^[19]），并会导致链的收敛速度较慢。 - 第三，调查人员经常遇到具有多个空间相关结果的多元空间数据集，其分析涉及矩阵计算的多元空间模型，而这些模型在 BUGS 中实施效果不佳。不过，随着通过 CRAN (http://cran.r-project.org) 统计计算环境中软件包的提供，这些问题开始减弱。从 CRAN 可以轻松获得几个用于点参考数据的贝叶斯方法自动化和 MCMC 算法收敛性诊断的软件包。这其中包括 `geoR` (Ribeiro & Diggle 2001 ^[41])、`geoRglm` (Christensen & Ribeiro 2002 ^[9])、`spTimer` (Bakar & Sahu 2014 ^[1])、`spBayes` (Finley 等， 2015 ^[20])、`spate` (Sigrist 等， 2015 ^[45]) 和 `ramps`（Smith 等，2008 ^[46]）。在本文介绍的分层地统计模型方面，`spBayes` 为用户提供了一套用于 “高斯” 和 “非高斯” 、“单变量” 和 “多变量”的空间数据贝叶斯分层模型以及动态贝叶斯时空模型。它侧重于完全贝叶斯推断的性能问题、采样器收敛速度、使用折叠吉布斯采样器的效率。它通过避免昂贵的矩阵计算来减少采样器的运行时间，通过实现 `预测过程模型` 来增加对大型数据集的可扩展性（第 6 节）。除了对现有模型的这些一般计算改进之外，它还采用一类动态时空模型及其预测过程对应物，来分析在空间和时间上索引的数据，适用于空间连续而时间离散的应用场景。 ### 7.2 多变量解释，多变量响应联合物种分布模型是当前最缺乏有效解决方案的地统计问题（参见如，Pollock 等，2014 ^[39]及其中的参考文献）。说明性数据可能来自大量空间测站（如 $10^3$ 规模）上的大量物种（如 $10^3$ 规模）的存在/不存在（二值数据类型）或丰度（计数型或覆盖百分比）。每个测站中都有一套环境预测器可用，并且在任何测站都会观测到大量 $0$ 值。人们已经达成共识，物种之间在共现或共丰度方面存在着依赖性，因此，在各测站处的响应变量也将存在空间依赖性。目前，尚需要非常严格的假设来解决这个问题。 ## 参考文献

[1] Bakar KS, Sahu SK. spTimer: spatio-temporal Bayesian modeling using R. J Stat Softw. 2014; 63:132.
[2] Banerjee, S., Carlin, BP., Gelfand, AE. Hierarchical Modeling and Analysis for Spatial Data. 2. Boca Raton, FL: Chapman and Hall/CRC Press; 2014.
[3] Banerjee S, Gelfand AE, Finley AO, Sang H. Gaussian predictive process models for large spatial datasets. J R Stat Soc B. 2008; 70:825–48.
[4] Barry RP, Ver Hoef JM. Blackbox kriging: spatial prediction without specifying variogram models. J Agric Biol Environ Stat. 1996; 1:297–322.
[5] Berger JO, DeOliveira V, Sansó B. Objective Bayesian analysis of spatially correlated data. J Am Stat Assoc. 2001; 96:1361–74.
[6] Berliner LM. Hierarchical Bayesian modeling in the environmental sciences. AStA Adv Stat Anal. 2000; 84:141–53.
[7] Berrocal VJ, Gelfand AE, Holland DM. A spatio-temporal downscaler for outputs from numerical models. J Agric Biol Environ Stat. 2010; 15:176–97.
[8] Chiles, JP., Delfiner, P. Geostatistics: Modeling Spatial Uncertainty. New York: Wiley; 1999.
[9] Christensen OF, Ribeiro PJ Jr. geoRglm: A package for generalised linear spatial models. RNews. 2002; 2:26–28.
[10] Clark JS, Gelfand AE. A future for models and data in ecology. Trends Ecol Evol. 2006; 21:375–80.
[11] Cressie N, Huang H-C. Classes of nonseparable spatio-temporal stationary covariance functions. J Am Stat Assoc. 1999; 94:1330–40.
[12] Cressie, N., Wikle, CK. Statistics for Spatiotemporal Data. New York: Wiley; 2011.
[13] Csató, L. PhD thesis. Aston Univ; Birmingham, UK: 2002. Gaussian processes—iterative sparse approximations.
[14] Datta A, Banerjee S, Finley AO, Gelfand AE. Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets. J Am Stat Assoc. 2016a; 111:800–12.
[15] Datta A, Banerjee S, Finley AO, Gelfand AE. On nearest-neighbor Gaussian process models for massive spatial data. WIREs Comput Stat. 2016c; (8):162–71.
[16] Datta A, Banerjee S, Finley AO, Hamm NAS, Schaap M. Nonseparable dynamic nearest neighbor Gaussian process models for large spatio-temporal data with an application to particulate matter analysis. Ann Appl Stat. 2016b; 10:1286–316.
[17] DeOliveira V. Bayesian prediction of clipped Gaussian random fields. Comput Stat Data Anal. 2000; 34:299–314.
[18] Diggle PJ, Tawn JA, Moyeed RA. Model-based geostatistics (with discussion). J R Stat Soc C. 1998; 47:299–350.
[19] Finley AO, Banerjee S, Carlin BP. spBayes: an R package for univariate and multivariate hierarchical point-referenced spatial models. J Stat Softw. 2007; 19:4.
[20] Finley AO, Banerjee S, Gelfand AE. spBayes for large univariate and multivariate point-referenced spatiotemporal data models. J Stat Softw. 2015; 63:1–28.
[21] Finley AO, Sang H, Banerjee S, Gelfand AE. Improving the performance of predictive process modeling for large datasets. Comput Stat Data Anal. 2009; 53:2873–84.
[22] Fuentes M, Raftery AE. Model evaluation and spatial interpolation by Bayesian combination of observations with outputs from numerical models. Biometrics. 2005; 61:36–45.
[23] Gelfand, AE., Diggle, P., Fuentes, M., Guttorp, P. Handbook of Spatial Statistics. Boca Raton, FL: CRC; 2010.
[24] Gelfand AE, Kim H-J, Sirmans CF, Banerjee S. Spatial modeling with spatially varying coefficient processes. J Am Stat Assoc. 2003; 98:387–96.
[25] Gelfand AE, Schmidt AM, Banerjee S, Sirmans CF. Nonstationary multivariate process modeling through spatially varying coregionalization (with discussion). Test. 2004; 13:263–312.
[26] Gelfand AE, Zhu L, Carlin BP. On the change of support problem for spatio-temporal data. Biostatistics. 2001; 2:31–45.
[27] Gneiting T. Nonseparable, stationary covariance functions for space-time data. J Am Stat Assoc. 2002; 97:590–600.
[28] Higdon DM. A process-convolution approach to modeling temperatures in the north Atlantic Ocean. J Environ Ecol Stat. 1998; 5:173–90.
[29] Higdon, DM. Space and space-time modeling using process convolutions. In: Anderson, C.Barnett, V.Chatwin, PC., El-Shaarawi, AH., editors. Quantitative Methods for Current Environmental Issues. London: Springer-Verlag; 2002. p. 37-56.
[30] Kalnay, E. Atmospheric Modeling, Data Assimilation and Predictability. Cambridge, UK: Cambridge Univ. Press; 2003.
[31] Katzfuss M. Bayesian nonstationary spatial modeling for very large datasets. Environmetrics. 2013; 24:189–200.
[32] Katzfuss, M. A multi-resolution approximation for massive spatial datasets. J Am Stat Assoc. 2016. http://dx.doi.org/10.1080/01621459.2015.1123632
[33] Krige DG. A statistical approach to some basic mine valuation problems on the Witwatersrand. J Chem Metall Min Soc S Afr. 1951; 52:119–39.
[34] Little, RJA., Rubin, DB. Statistical Analysis with Missing Data. New York: Wiley-Interscience; 2002.
[35] Loh W-H. Fixed-domain asymptotics for a subclass of Matérn-type Gaussian random fields. Ann Stat. 2005; 33:2344–94.
[36] Mardia, KV., Goodall, C. Spatio-temporal analyses of multivariate environmental monitoring data. In: Patil, GP., Rao, CR., editors. Multivariate Environmental Statistics. Amsterdam: Elsevier; 1993. p. 347-86.
[37] Murphy, KP. Machine Learning: A Probabilistic Perspective. Cambridge, MA: MIT Press; 2012.
[38] Paciorek CJ, Schervish M. Spatial modelling using a new class of nonstationary covariance functions. Environmetrics. 2006; 17:483–506.
[39] Pollock LJ, Tingley R, Morris WK, Golding N, O’Hara RB, et al. Understanding co-occurrence by modelling species simultaneously with a Joint Species Distribution Model ( JSDM). Methods Ecol Evol. 2014; 5:397–406.
[40] Poole D, Raftery AE. Inference for deterministic simulation models: the Bayesian melding approach. J Am Stat Assoc. 2000; 95:1244–55.
[41] Ribeiro PJ Jr, Diggle PJ. geoR: A package for geostatistical analysis. R News. 2001; 1/2:15–18.
[42] Robert, CP., Casella, G. Monte Carlo Statistical Methods. 2. New York: Springer-Verlag; 2004.
[43] Sang H, Huang JZ. A full scale approximation of covariance functions for large spatial data sets. J R Stat Soc B. 2012; 74:111–32.
[44] Sang H, Jun M, Huang JZ. Covariance approximation for large multivariate spatial data sets with an application to multiple climate model errors. Ann Appl Stat. 2011; 4:2519–48.
[45] Sigrist F, Kuensch HR, Stahel WA. spate: An R package for spatio-temporal modeling with a stochastic advection-diffusion process. J Stat Softw. 2015; 63:1–23.
[46] Smith BJ, Yan J, Cowles MK. Unified geostatistical modeling for data fusion and spatial heteroskedasticity with R package ramps. J Stat Softw. 2008; 25:1–21.
[47] Stein, ML. Interpolation of Spatial Data: Some Theory for Kriging. New York: Springer-Verlag; 1999a.
[48] Stein ML. Predicting random fields with increasingly dense observations. Ann Appl Probability. 1999b; 9:242–73.
[49] Stein ML. Space-time covariance functions. J Am Stat Assoc. 2005; 100:310–21.
[50] Stein ML. Limitations on low rank approximations for covariance matrices of spatial data. Spatial Stat. 2014; 8:1–19.
[51] Stein ML, Chi Z, Welty LJ. Approximating likelihoods for large spatial datasets. J R Stat Soc B. 2004; 66:275–96.
[52] Vecchia AV. Estimation and model identification for continuous spatial processes. J R Stat Soc B. 1988; 50:297–312.
[53] Ver Hoef JM, Barry RP. Constructing and fitting models for cokriging and multivariable spatial prediction. J Stat Plann Inference. 1998; 69:275–94.
[54] Wackernagel H. Cokriging versus kriging in regionalized multivariate data analysis. Geoderma. 1994; 62:83–92.
[55] Wackernagel, H. Multivariate Geostatistics: An Introduction with Applications. 3. New York: Springer-Verlag; 2003.
[56] West, M., Harrison, PJ. Bayesian Forecasting and Dynamic Models. 2. New York: Springer-Verlag; 1997.
[57] Wikle, CK. Low-rank representations for spatial processes. In: Gelfand, AE.Diggle, P.Fuentes, M., Guttorp, P., editors. Handbook of Spatial Statistics. Boca Raton, FL: Chapman and Hall/CRC; 2010. p. 107-18.
[58] Xia, G., Gelfand, AE. Tech Rep. ISDS, Duke Univ; 2006. Stationary process approximation for the analysis of large spatial datasets.
[59] Yaglom, AM. Correlation Theory of Stationary and Related Random Functions: Volume I: Basic Results. London: Springer; 1987.
[60] Zhang H, Zimmerman DL. Towards reconciling two asymptotic frameworks in spatial statistics. Biometrika. 2005; 92:921–36.