spBayes--贝叶斯空间变系数模型的 R 软件包

【摘 要】 本文描述并说明了在 spBayes（版本 0.4-2）R 包中拟合空间变系数模型的新功能。新的 spSVC 函数使用计算效率高的马尔可夫链蒙特卡罗算法，并扩展了当前仅适用于空间变化截距回归模型的 spBayes 函数，以适用于回归设计矩阵中任何一组列的独立或多元高斯过程随机效应。讨论和说明了新添加的用于 spSVC 的 OpenMP 并行化选项，以及用于联合和逐点预测和模型拟合诊断的辅助函数。使用中欧的 PM10 分析说明了所提出模型的效用。

【原 文】 A. O. Finley and S. Banerjee, “Bayesian spatially Varying coefficient models in the spBayes R package,” Environmental Modelling & Software, vol. 125, p. 104608, 2020, doi: 10.1016/j.envsoft.2019.104608.

1 简介

在本文中，我们描述并说明了 spBayes（Finley、Banerjee 和 Gelfand 2015）R 包中核心函数最近重新设计和重写的扩展功能。 spBayes 包为空间索引的高斯和非高斯响应提供了一套单变量和多变量回归模型。

到目前为止，有许多提供类似功能的 R 包。最近的 “空间数据分析” CRAN 任务视图（Bivand 2019）列出了 $46$ 个用于地统计分析的包，但这也并不是可用于此类分析的完整列表。Finley 等 (2015) 专注于制定计算高效且灵活的 MCMC 算法，以估计一系列时空高斯过程 (GP) 模型，不过当时虽然提出的一种通用的采样算法，但在 spBayes 中只实现了一小部分模型。具体来说，用户只能在模型的截距项上指定单变量（或多变量）高斯过程。现在，我们将 spSVC 函数添加到了 spBayes 中，可以将单变量（或多变量）高斯过程放置在任意回归系数（集合）上，给用户提供更多的选项。

此功能并不是唯一的，有几个 R 包也能够拟合空间变系数 (SVC) 模型。其中有些是专门为空间或时空数据而设计的，而另一些则提供了更多灵活性和通用性，允许系数随一些通用变量集变化（当然也包括对坐标索引的支持）。这些包中的大多数都采用了某种基于样条或核的回归方法，以允许预测变量产生不同的效应。 在理论方法层面， Hastie 和 Tibshirani (1993) 以及 Fan 和 Zhang (2008) 提供了通用变系数模型的发展，Gelfand、Kim、Sirmans 和 Banerjee (2003a) 提供了特定于空间设置的处理方法。

关于 R 中的实现：

• 地理加权回归：spgwr 包（Bivand 和 Yu 2017）实现了地理加权回归，如 Fotheringham、Brunsdon 和 Charlton（2002 年）中最初详述的那样。
• 样条方法：基于样条的包选项包括 mgcv (Wood 2017)、svcm (Heim 2007)、np (Hayfield and Racine 2008) 和 mboost (Hothorn、Buehlmann、Kneib、Schmid 和 Hofner 2018)。
• 基于树的方法：Bürgin 和 Ritschard（2017 年）最近开发了一种基于树的变系数模型 (TVCM) 算法和相关的 vcrpart 包。
• 时空数据：walker（Helske 2019；Vihola、Helske 和 Franks 2017）、spTimer （ Bakar 和 Sahu（2018）和 spTDyn（Bakar、Kokic 和 Jin 2017、2015a、b）提供贝叶斯时间和时空 SVC 模型。
• 基于通用软件包的实现：
  • INLA（Rue、Martino 和 Chopin 2009；Lindgren 和 Rue 2015；Bakka、Rue、Fuglstad、Riebler、Bolin、Illian、Krainski、Simpson 和 Lindgren 2018）
  • Stan（ Carpenter、Gelman、Hoffman、Lee、Goodrich、Betancourt、Brubaker、Guo、Li 和 Riddell 2017 年；Stan 开发团队 2018 年），可以从各自的 R 包中调用。

原文作者在此处没有提 基于特征向量的空间滤波方法 及其典型软件包 spmoran，有关空间滤波方面的文献请参见：

• Griffith, D.A. (2003) Spatial autocorrelation and spatial filtering: gaining understanding through theory and scientific visualization. Springer Science & Business Media.
• Murakami, D. and Griffith, D.A. (2015) Random effects specifications in eigenvector spatial filtering: a simulation study. Journal of Geographical Systems, 17 (4), 311-331.
• Murakami, D. and Griffith, D.A. (2019a) Eigenvector spatial filtering for large data sets: fixed and random effects approaches. Geographical Analysis, 51 (1), 23-49.
• Murakami, D. and Griffith, D.A. (2019b) Spatially varying coefficient modeling for large datasets: Eliminating N from spatial regressions. Spatial Statistics, 30, 39-64.

软件包见：

• D. Murakami, “Spatial regression modeling using the spmoran package: Boston housing price data examples,” 2017. arXiv:1703.04467v11

spBayes 包中的 spSVC 函数使用了高效的采样算法，提供基于马尔可夫链蒙特卡罗 (MCMC) 的 SVC 推断。该算法的效率源于仅对协方差参数的更新（即回归系数和随机效应被积分出来）、计算并行化以及调谐和/或多线程矩阵代数库的使用。

2 模型和软件

设 $y(\mathbf{s})$ 为位置 $\mathbf{s}$ 处的因变量（响应或结果），现在考虑空间变化的回归模型，

$$y(\mathbf{s}) = (\beta_1 + \delta_1 w_1 (\mathbf{s})) + \sum^p_{j=2} x_j(\mathbf{s}) \{\beta_j + \delta_j w_j(\mathbf{s})\} + (\mathbf{s}) \tag{1}$$

$x_j(\mathbf{s})$ 是位置 $\mathbf{s}$ 处预测变量的已知值（对于每个 $j = 2，\ldots, p, p ≥ 1$ ）, $\beta_j$ 是对应于 $x_j(\mathbf{s})$ 的回归系数, $\beta_1$ 是截距, $\epsilon(\mathbf{s})$ 是对所有 $\mathbf{s}$ 独立分布的高斯测量误差过程。 $w_1 (\mathbf{s})$ 和 $w_j(\mathbf{s})$ 是对应于截距和预测变量的空间过程，从而产生空间变化的回归模型。

我们进一步考虑到并非所有预测变量都会对结果产生空间变化影响。因此，式 (1) 中的 $\delta$ 是二值的指示器，如果相关预测变量具有空间变化的回归系数，则 $\delta = 1$，否则为$\delta = 0$。为了后面的方便，当相应的 $\delta = 1$ 时，我们定义 $\tilde{\beta}_1(\mathbf{s}) = \beta_1 + \delta_1 w_1 (\mathbf{s})$ 和 $\tilde{\beta}_j(\mathbf{s}) = \beta_j + \delta_j w_j(\mathbf{s})$ 作为空间变化回归系数。让 \mathcal{S} = \{s_{1}, s_{2}, \ldots , s_n} 是观测到 $y(\mathbf{s}_i)$ 和预测变量的空间位置集。设 $\mathbf{w}$ 是通过堆叠 $\mathbf{w}(\mathbf{s}_i)$ 获得的 $nr × 1$ 向量，其中每个 $\mathbf{w}(\mathbf{s}_i)$ 是一个 $r × 1$ 向量，具有第 $j$ 个条目 $w_j(\mathbf{s}_i)$，$j = 1, 2, \ldots , r$ 并且 $i = 1, 2,\ldots, n$。我们将 $\mathbf{w}(\mathbf{s})$ 视为多元高斯过程（例如，参见 Banerjee、Carlin 和 Gelfand 2014），因此矩阵 $\mathbf{K}_{\theta}$ 是一个 $nr × nr$ 空间协方差矩阵，构造为具有第 $(i, j)$ 个块的块矩阵从指定多元空间过程 $\mathbf{w}(\mathbf{s})$ 的 $r × r$ 交叉协方差矩阵 $\mathbf{K}_{\theta}(\mathbf{s}_i, \mathbf{s}_j)$ 获得。另外，$\boldsymbol{\beta}$ 是 $\mathbf{X}$ 对应的 $p×1$ 回归系数，$\boldsymbol{\theta}$ 和 $\boldsymbol{\tau}$ 分别是 $\mathbf{K}_{\boldsymbol{\theta}}$ 和 $\mathbf{D}_{\boldsymbol{\tau}}$ 中的参数。

考虑从 式 (1) 构建的贝叶斯分层模型：

$$p({\boldsymbol{\tau}} , {\boldsymbol{\theta}}) × N ({\boldsymbol{\beta}} | {\boldsymbol{\mu}}_{\boldsymbol{\beta}}, \boldsymbol{\Sigma}_{\boldsymbol{\beta}}) × N (\mathbf{w} | \mathbf{0}, \mathbf{text}_{\boldsymbol{\theta}}) × N (\mathbf{y} | \mathbf{X}{\boldsymbol{\beta}} + \mathbf{Zw}, \mathbf{D}_{\boldsymbol{\tau}} ) \tag{2}$$

其中 $\mathbf{y}$ 是 $n × 1$，第 $i$ 个元素 $y(\mathbf{s}_i)$，$\mathbf{X}$ 是 $n × p$ 矩阵，第一列为 $\mathbf{1}$，其余 $p − 1$ 列对应于 式 (1) 中的预测变量 $x_j(\mathbf{s}_i)$。矩阵 $\mathbf{Z}$ 是 $n × nr$ 矩阵，其中 $r ≤ p$，恰好是 $\mathbf{X}$ 中 $\delta_j = 1$ 的那些列。

我们为用户提供了矩阵 $\mathbf{X}$ 的缩放和中心化功能。请注意，$\mathbf{Z}$ 是从 $\mathbf{X}$ 构建的，因此，对于缩放和中心化的 $\mathbf{X}$，$\mathbf{Z}$ 中使用的预测变量也将缩放和中心化。缩放和中心化通常可以提高数值稳定性，并为空间变化的回归模型提供更稳健的估计（例如，参见 Gelfand、Kim、Sirmans 和 Banerjee 2003b）。

一些进一步的规范是有序的。在 spSVC 中，我们将固定 $D_{\tau} = {\tau}^2\mathbf{I}_n$，因此 {\boldsymbol{\tau}} = \{\tau}^2\} 是表示测量误差方差或地统计学中的 “块金” 的标量。互协方差矩阵 $\mathbf{K}_{\boldsymbol{\theta}}(\mathbf{s}, \mathbf{t})$，其中 $\mathbf{t}$ 是一个普通的位置，通常使用共区域化线性模型 (LMC) 建模。在这里，我们将建模 $K_{\boldsymbol{\theta}}(\mathbf{s,t}) = \mathbf{A} \boldsymbol{\Gamma}(\mathbf{s,t})\mathbf{A}^T$，其中 $\mathbf{A}$ 是 $r × r$ 下三角矩阵，$\boldsymbol{\Gamma}(\mathbf{s,t})$ 是对角矩阵，其中 $\rho_j(\mathbf{s,t})$ 是第 $j$ 个对角线元素，其中 $\rho_j(\mathbf{s,t})$ 是一个空间相关函数，其参数特定于 $w_j(\mathbf{s})$。这里，式 (2) 中的 ${\boldsymbol{\theta}}$ 对应于 $\{\mathbf{A}, \{\boldsymbol{\phi}_j\}^j_{r=1}\}$，其中每个 $\boldsymbol{\phi}_j$ 是空间相关函数中参数的集合。例如，对于 Matérn 协方差函数，每个 $\boldsymbol{\phi}_j$ 都包含一个空间衰减参数和一个平滑度参数。

任何位置 $\mathbf{s}$ 内 $\mathbf{w}(\mathbf{s})$ 的协方差结构由 $\mathbf{AA}^T$ 捕获，它与 $\operatorname{Var}{\mathbf{w}(\mathbf{s})}$ 的 Cholesky 分解一致。通常，我们将先验指定为

$$p({\tau}^2, {\boldsymbol{\theta}}) = \text{IG}({\tau}^2 | a_{\tau} , b_{\tau} ) × \text{IW}(\mathbf{AA}^T | r_a, \mathbf{S}_a) × \prod^r_{j=1} p(\boldsymbol{\phi}_j) \tag{3}$$

其中 $\text{IG}$ 是 逆 Gamma 分布， $\text{IW}$ 是 逆 Wishart分布，每个 $p(\boldsymbol{\phi}_j)$ 可以是 spBayes 提供的几种分布之一。 spSVC 提供的另一个特定选择指定 $\mathbf{A} = \text{Diag}(σ_1,\ldots, σ_r)$，因此 $\mathbf{K}_{\boldsymbol{\theta}}(\mathbf{s,t})$ 是对角线元素为 $\sigma^2_j \rho_j(\mathbf{s,t})$ 成对角矩阵（ $j = 1, 2, \ldots,r$ ），在这种情况下我们假设 $\text{IG}(\sigma^2_j | a_σ, b_σ)$ , $\rho_j$ 的选择包括 spBayes 提供的任何标准相关函数。

2.1 参数估计和计算注意事项

(1) 的贝叶斯推断涉及从参数 ${\boldsymbol{\theta}}$、${\boldsymbol{\beta}}$ 和 $\mathbf{w}$ 的边缘后验分布中对其进行采样。这种采样算法需要对密集矩阵进行昂贵的运算，例如分解和乘法。因此需要注意使用有效的数值算法（例如 Cholesky 分解）、使用三角系统并避免冗余操作。

2.1.1 采样过程参数

从 式 (2) 中采样采用 MCMC 方法，特别是 Gibbs 采样和随机游走 Metropolis 步骤（例如 Robert 和 Casella 2004）。为了更快收敛，我们从模型中整合出 ${\boldsymbol{\beta}}$ 和 $\mathbf{w}$，并从 $p({\boldsymbol{\theta}} | \mathbf{y}) \propto p({\boldsymbol{\theta}}) × N (\mathbf{y} | \mathbf{X} {\boldsymbol{\mu}}_{\boldsymbol{\beta}}, \boldsymbol{\Sigma}_{\mathbf{y}|{\boldsymbol{\theta}}})$ 中提取第一个样本， 其中 $\boldsymbol{\Sigma}_{\mathbf{y}| {\boldsymbol{\theta}}} = \mathbf{X} \boldsymbol{\Sigma}_{\boldsymbol{\beta}}\mathbf{X}^T + \mathbf{ZK}_{\boldsymbol{\theta}}\mathbf{Z}^T + \mathbf{D}_{\boldsymbol{\tau}}$ 。每次更新 ${\boldsymbol{\theta}}$ 都需要构建这个矩阵。 $\mathbf{D}_{\tau}$ 是对角矩阵，$\mathbf{X} \boldsymbol{\Sigma}_{\boldsymbol{\beta}}\mathbf{X}^T$ 是固定的，所以计算涉及到矩阵 $\mathbf{ZK}_{\boldsymbol{\theta}} \mathbf{Z}^T$ 需要 $rn^2$ 个 flops（浮点运算）。

在转换参数以支持整个实线后，我们采用带有多元正态提议（与 ${\boldsymbol{\theta}}$ 中的参数相同的维度）的随机游走 Metropolis 步骤。这涉及估值

\log p({\boldsymbol{\theta}} | \mathbf{y}) = \text{const} + \log p({\boldsymbol{\theta}}) − \frac{1}{2} \log |\boldsymbol{\Sigma}_{y | {\boldsymbol{\theta}}} | − \frac{1}{2} Q({boldsymbol{\theta}) \tag{4}

其中 $Q(\boldsymbol{\theta}) = (\mathbf{y} − \mathbf{X}{\boldsymbol{\mu}}_{\boldsymbol{\beta}})^T \boldsymbol{\Sigma}^{−1}_{y | \theta}(\mathbf{y} − \mathbf{X}{\boldsymbol{\mu}}_{\beta})$。

通常，我们计算 $L = chol(\boldsymbol{\Sigma}y | {\boldsymbol{\theta}})$，其中 $chol(\boldsymbol{\Sigma}y | {\boldsymbol{\theta}})$ 返回 $\boldsymbol{\Sigma}_{\mathbf{y} |\theta}$ 的下三角 Cholesky 因子 $\mathbf{L}$ 。这涉及 $O(^n3/3)$ 个 flops。接下来，我们得到 $u = trsolve(L, y − X{\boldsymbol{\mu}}{\boldsymbol{\beta}})$，它求解三角方程组 $Lu = y − X{\boldsymbol{\mu}}{\boldsymbol{\beta}}$。这涉及 $O(n^2)$ 次 flops，而 $Q({\boldsymbol{\theta}}) = u^Tu$ 需要另外 $2n$ 次触发器。 式 (4) 中的对数行列式计算为 $2 \Sigma^n_{i=1} \log l_{i,i}$，其中 $l_{i,i}$ 是 $\mathbf{L}$ 中的对角线项。由于 $\mathbf{L}$ 已经获得，对数行列式还需要 $n$ 步。因此，Cholesky 分解占主导地位，计算 式 (4) 在 $O(n^3)$ flops 中实现。如果 ${\boldsymbol{\beta}}$ 是平坦的，即 $\boldsymbol{\Sigma}^{−1}_{\boldsymbol{\beta}} = \mathbf{O}$，则式 (4) 分布的模拟是

$$\log p({\boldsymbol{\theta}} | \mathbf{y}) = \text{const} + \log p({\boldsymbol{\theta}}) − \frac{1}{2} \log |\mathbf{X}^T \boldsymbol{\Sigma}^{−1}_{y | {\beta},{\theta}} \mathbf{X}| − \frac{1}{2} \log |\boldsymbol{\Sigma}_{y | {\beta},{\theta}}| − \frac{1}{2} Q({\boldsymbol{\theta}}) \tag{5}$$

其中\boldsymbol{\Sigma}{y | {\beta},{\theta}} = \mathbf{ZK}{\boldsymbol{\theta}} \mathbf{Z}^T + \mathbf{D}_{\tau} 和 $Q({\boldsymbol{\theta}}) = \mathbf{y}^T \boldsymbol{\Sigma}^{−1}_{y | {\beta},{\theta}} \mathbf{y} − b^T(\mathbf{X}^T \boldsymbol{\Sigma}^{−1}_{y |{\beta},{\theta}} \mathbf{X})^{−1}\mathbf{b}$ 和 $\mathbf{b} = \mathbf{X}^T \boldsymbol{\Sigma}^{−1}_{\mathbf{y} | {\beta},{\theta}} \mathbf{y}$。计算过程与上述类似。我们首先评估 $\mathbf{L} = chol(\boldsymbol{\Sigma}_{y | {\beta},{\theta}})$，然后获得 $[\mathbf{v:U}] = trsolve(\mathbf{L}, [\mathbf{y} :\mathbf{X}])$，因此 $\mathbf{Lv} = \mathbf{y}$ 和 $\mathbf{LU} = \mathbf{X}$。接下来，我们评估 $\mathbf{W} = chol( \mathbf{U}^T \mathbf{U})$, $\mathbf{b} = \mathbf{U}^T \mathbf{v}$ 并求解 $\mathbf{b} = trsolve(\mathbf{W,b})$。最后，式（5） 被评估为

$$\log p({\boldsymbol{\theta}}) − \sum^p_{ i=1 } \log w_{i,i} − \sum^n_{i=1} \log l_{i,i} − \frac{1}{2} (\mathbf{v^Tv} − \tilde{\mathbf{b}}^T \tilde{\mathbf{b}})$$

其中 $w_{i,i}$ 和 $l_{i,i}$ 分别是 $\mathbf{W}$ 和 $\mathbf{L}$ 中的对角线元素。flops 的数量也是 $n$ 的立方阶。

重要的是，我们上面的策略避免了计算逆。我们使用 Cholesky 分解并仅求解三角系统。如果 $n$ 不大，比如说 \sim $10^2$，这个策略是可行的。如 第 2.2 节 所述和 第 3 节 所示，使用高效和并行化的数值线性代数例程可以显着缩短计算时间。

2.1.2 采样斜率和随机效应

一旦我们从 $p(\boldsymbol{\theta} | \mathbf{y})$ 中获得了边缘后验样本 $\boldsymbol{\theta}$，我们就可以使用 组合采样 抽取 $\boldsymbol{\beta}$ 和 $\mathbf{w}$ 的后验样本。假设 $\{\boldsymbol{\theta}^{(1)}, \boldsymbol{\theta}^{(2)}, \ldots, \boldsymbol{\theta}^{(M)}\}$ 是来自 $p({\boldsymbol{\theta}} | \mathbf{y})$ 的 $M$ 个样本。抽取 ${\boldsymbol{\beta}}^{(k)} \sim p({\boldsymbol{\beta}} | {\boldsymbol{\theta}}^{(k)}, \mathbf{y})$ 和 $\mathbf{w}^{(k)} \sim p(\mathbf{w} | {\boldsymbol{\theta}}^{(k)}, \mathbf{y})$，$k = 1, 2, \ldots, M$ 分别从 $p({\boldsymbol{\beta}} | \mathbf{y})$ 和 $p(\mathbf{w} | \mathbf{y})$ 产生 $M$ 个样本。我们只需要保留和存储 MCMC 算法收敛以后（即老化后）获得的 ${\boldsymbol{\theta}}$ 样本。

为了进一步阐明，请注意 ${\boldsymbol{\beta}} | {\boldsymbol{\theta}}, \mathbf{y} \sim N_p(\mathbf{Bb,B})$ 具有均值 $\mathbf{Bb}$ 和方差-协方差矩阵 $B$，其中

$$\begin{align*} \mathbf{b} &= \boldsymbol{\Sigma}^{-1}_{\boldsymbol{\beta}} {\boldsymbol{\mu}}_{\boldsymbol{\beta}} + \mathbf{X}^T \boldsymbol{\Sigma}^{-1}_{\mathbf{y} | {\boldsymbol{\beta}},{\boldsymbol{\theta}}} \mathbf{y}\\ B &= ( \boldsymbol{\Sigma}^{-1}_{\boldsymbol{\beta}} + \mathbf{X}^T \boldsymbol{\Sigma}^{-1}_{y | {\boldsymbol{\beta}},{\boldsymbol{\theta}}} \mathbf{X})^{−1} \tag{6} \end{align*}$$

对于每个 $k = 1, 2, \ldots, M$，我们在当前值 ${\boldsymbol{\theta}}^{(k)}$ 处计算 $\mathbf{B}$ 和 $\mathbf{b}$ 并得出 ${\boldsymbol{\beta}}^{(k)} \sim N_p(\mathbf{Bb,B})$。这是通过计算 $\mathbf{b} = \boldsymbol{\Sigma}^{-1}_{\boldsymbol{\beta}} {\boldsymbol{\mu}}_{\boldsymbol{\beta}} + \mathbf{U}^T \mathbf{v}$ 实现的，其中 $\mathbf{L} = chol(\boldsymbol{\Sigma}_{y | {\boldsymbol{\beta}},{\boldsymbol{\theta}}^{(k)}})$ 和 $[\mathbf{v:U}] = trsolve(\mathbf{L}, [\mathbf{y:X}])$。接下来，我们生成 $p$ 个独立的标准正态变量，将它们收集到 $\mathbf{z}$ 中并设置

$${\boldsymbol{\beta}}^{(k)} = trsolve ( \mathbf{L}^T_{B}, trsolve(\mathbf{L}_{B}, \mathbf{b}) ) + trsolve(\mathbf{L}^T_{B}, \mathbf{z}) \tag{7}$$

其中 $\mathbf{L}_B = chol ( \boldsymbol{\Sigma}^{-1}_{\boldsymbol{\beta}} + \mathbf{U}^T \mathbf{U})$ 。这样就完成了第 $k$ 次迭代。经过 $M$ 次迭代后，我们得到 $\{ \boldsymbol{\beta}^{(1)}, \boldsymbol{\beta}^{(2)},\ldots , \boldsymbol{\beta}^{(M)} \}$，它们是来自 $p(\boldsymbol{\beta} | \mathbf{y})$ 的样本。

通过对空间随机效应的点或区间估计制图，通常有助于识别缺失的回归变量和/或更好地理解模型的充分性。 $\boldsymbol{\Sigma}_{\mathbf{y}| w,{\boldsymbol{\theta}}} = \mathbf{X} \boldsymbol{\Sigma}_{\boldsymbol{\beta}} \mathbf{X}^T + \mathbf{D}_{\boldsymbol{\tau}}$ 并注意 $\mathbf{w} | {\boldsymbol{\theta}}, \mathbf{y} \sim N (\mathbf{Bb,B})$，其中

$$\begin{align*} \mathbf{b} &= \mathbf{Z}^T \boldsymbol{\Sigma}^{-1}_{y|w,{\boldsymbol{\theta}}}(\mathbf{y} − \mathbf{X} {\boldsymbol{\mu}}_{\boldsymbol{\beta}})\\ B &= ( \mathbf{K}^{-1}_{\boldsymbol{\theta}} + \mathbf{Z}^T \boldsymbol{\Sigma}^{-1}_{y | w,{\boldsymbol{\theta}}} \mathbf{Z})^{-1} \tag{8} \end{align*}$$

这里的向量 $b$ 的计算类似于 ${\boldsymbol{\beta}}$。对于每个 $k = 1, 2, \ldots , M$ 我们现在评估 $\mathbf{L} = chol(\boldsymbol{\Sigma}_{y | w,{\boldsymbol{\theta}}^{(k)}}), [\mathbf{v:U}] = trsolve(\mathbf{L}, [\mathbf{y} − \mathbf{X} {\boldsymbol{\mu}}_{\boldsymbol{\beta}} : \mathbf{Z}({\boldsymbol{\theta}}^{(k)})])$ 并设置 $\mathbf{b} = \mathbf{U}({\boldsymbol{\theta}}^{(k)})^T \mathbf{v}$。对于计算 $\mathbf{B}$，可以按照 ${\boldsymbol{\beta}}$ 的方式进行，但这会涉及 $chol(K({\boldsymbol{\theta}}))$，对于某些协方差函数（例如，具有大 $ν$ 的高斯函数或 Matérn），它可能在数值上变得不稳定。对于稳健的软件性能，我们定义 $G^{-1}_{\boldsymbol{\theta}} = \mathbf{Z}^\prime \boldsymbol{\Sigma}^{-1}_{y | w,{\boldsymbol{\theta}}} \mathbf{Z}$ 并利用恒等式（Henderson 和 Searle 1981）

$$( \mathbf{K}^{-1}_{\boldsymbol{\theta}} + \mathbf{G}^{-1}_{\boldsymbol{\theta}})^{-1} = \mathbf{G}_{\boldsymbol{\theta}} − \mathbf{G}_{\boldsymbol{\theta}} (\mathbf{K}_{\boldsymbol{\theta}} + \mathbf{G}_{\boldsymbol{\theta}})^{-1} \mathbf{G}_{\boldsymbol{\theta}}$$

以设计一个数值稳定的算法。对于每个 $k = 1, 2, \ldots, M$ ，我们评估 $L = chol(\mathbf{K}^{(k)}_{\boldsymbol{\theta}} + \mathbf{G}^{(k)}_{\boldsymbol{\theta}})$, $\mathbf{W} = trsolve(\mathbf{L}, \mathbf{G}^{(k)}_{\boldsymbol{\theta}} )$ 和 $\mathbf{L}_B = chol(\mathbf{G}^{(k)}_{\boldsymbol{\theta}} − \mathbf{W}^T \mathbf{W})$。如果 $z$ 是独立标准正态变量的 $r × 1$ 向量，则我们设置 $\mathbf{w}^{(k)} = \mathbf{L}_B \mathbf{L}^T_B \mathbf{b} + \mathbf{L}_B \mathbf{z}$。结果 $\{\mathbf{w}^{(1)}, \mathbf{w}^{(2)}, \ldots , \mathbf{w}^{(M)}\}$ 是来自 $p(\mathbf{w} | \mathbf{y})$ 的样本。

2.1.3 空间预测

为了预测与 $n_0 × p$ 预测变量矩阵 $\mathbf{X}_0$ 关联的随机 $n_0 × 1$ 向量 $\mathbf{y}_0$，我们假设

$$\left[\begin{array}{c} \mathbf{y} \\ \mathbf{y}_0 \end{array}\right] \mid \boldsymbol{\beta}, \boldsymbol{\theta} \sim N_{n_0+n}\left(\left[\begin{array}{c} \mathbf{X} \\ \mathbf{X}_0 \end{array}\right] \boldsymbol{\beta},\left[\begin{array}{cc} \mathbf{C}_{11}(\boldsymbol{\theta}) & \mathbf{C}_{12}(\boldsymbol{\theta}) \\ \mathbf{C}_{12}(\boldsymbol{\theta})^{\top} & \mathbf{C}_{22}(\boldsymbol{\theta}) \end{array}\right]\right) \tag{9}$$

其中 $\mathbf{C}_{11}(\boldsymbol{\theta})=\boldsymbol{\Sigma}_{y \mid \beta, \theta}, \mathbf{C}_{12}(\boldsymbol{\theta})$ 是 $\mathbf{y}$ 和 $\mathbf{y}_0$ 之间的 $n × n_0$ 交叉协方差矩阵， $\mathbf{C}_{22}(\boldsymbol{\theta})$ 是 $\mathbf{y}_0$ 的方差-协方差矩阵。一个有效的联合分布将提供条件分布 $p(\mathbf{y}_0 | \mathbf{y}, {\boldsymbol{\beta}}, {\boldsymbol{\theta}})$，它服从以下均值和方差

$$\begin{align*} \boldsymbol{\mu}_p &=\mathbf{X}_0 \boldsymbol{\beta}+\mathbf{C}_{12}(\boldsymbol{\theta})^{\top} \mathbf{C}_{11}(\boldsymbol{\theta})^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) \\ \boldsymbol{\Sigma}_p &=\mathbf{C}_{22}(\boldsymbol{\theta})-\mathbf{C}_{12}(\boldsymbol{\theta})^{\top} \mathbf{C}_{11}(\boldsymbol{\theta})^{-1} \mathbf{C}_{12}(\boldsymbol{\theta}) \tag{10} \end{align*}$$

贝叶斯预测通过从后验预测分布 $p(\mathbf{y}_0 | \mathbf{y}) = \int p(\mathbf{y}_0 | \mathbf{y}, \boldsymbol{\beta}, \boldsymbol{\theta}) p(\boldsymbol{\beta}, \boldsymbol{\theta} | \mathbf{y}) d \boldsymbol{\beta} d \boldsymbol{\theta}$ 中采样来进行。对于 $\{ \boldsymbol{\beta}, \boldsymbol{\theta} \}$ 的每个后验样本，我们抽取一个相应的 $\mathbf{y}_0 \sim N (\boldsymbol{\mu}_p, \boldsymbol{\Sigma}_p)$。这会从后验预测分布中产生样本。

后验预测计算仅涉及收敛后保留的 MCMC 样本。对于任何后验样本 $\{ \boldsymbol{\beta}^{(k)}, \boldsymbol{\theta}^{(k)} \}$，我们求解 $[\mathbf{u:V}] = trsolve(\mathbf{L}, [\mathbf{y} − \mathbf{X} \boldsymbol{\beta}^{(k)} : \mathbf{C}_{12}({\boldsymbol{\theta}}^{(k)})])$，其中 $\mathbf{L} = chol (\mathbf{C}_{11}({\boldsymbol{\theta}}^{(k)}))$。接下来，我们设置 ${\boldsymbol{\mu}}^{(k)}_p = \mathbf{X}_0 {\boldsymbol{\beta}}^{(k)} + \mathbf{V}^T \mathbf{u}$ 和 $\boldsymbol{\Sigma}^{(k)}_p= \mathbf{C}_{22}({\boldsymbol{\theta}}^{(k)}) − \mathbf{V}^T \mathbf{V}$ 并抽取 $y^{(k)}_0 \sim N ({\boldsymbol{\mu}}^{(k)}_p , \boldsymbol{\Sigma}^{(k)}_p)$ 。

更新 $\mathbf{y}^{(k)}_0$ 需要对 $\boldsymbol{\Sigma}_p$ 进行 Cholesky 分解，它是 $n_0 × n_0$ 的矩阵，如果 $n_0$ 很大，可能会很昂贵。在大多数实际设置中，取 $n_0 = 1$ 并执行逐点预测就足够了。

2.2 spBayes 软件特色

spSVC 函数包含 spBayes 中的 spLM 函数，并提供额外的用户选项以简化分析和推断。下面的列表重点介绍了其中一些新选项。

1. 任何一组预测变量都可以接收独立的单变量 GP 或多变量 GP。

2. 可以通过从联合或逐点（边缘）后验预测分布中采样来进行预测。

3. openMP（Dagum 和 Menon 1998）支持可通过参数估计、成分采样、模型拟合诊断和预测函数的 n.omp.threads 参数获得

4. 矩阵运算并行化可通过基本线性代数子程序（BLAS；www.netlib.org/blas）和线性代数包（LAPACK；www.netlib.org/lapack）的多线程实现来实现。

5. 用于索引观测和预测位置的坐标系统可以是任意维度——用户以前仅限于使用二维系统。

6. 单变量和多变量随机效应样本和空间变系数作为列表返回，其中元素名称对应于给定的预测变量。

3 案例展示

我们考虑两个分析来说明 spSVC 的关键特性以及支持功能。第一个分析是针对模拟数据集，第二个分析是之前在 Hamm、Finley、Schaap 和 Stein (2015) 以及 Datta、Banerjee、Finley、Hamm 和 Schaap (2016) 中分析过的空气污染数据集。

3.1 模拟数据分析

模拟数据 mvSVCData 在 spBayes 中可用，包含分布在二维单位正方形空间域内的 $n=500$ 个观测值。在某一位置 $\mathbf{s}$ 处，结果生成如下

$$y(\mathbf{s}) = \beta_0 + w_0 (\mathbf{s}) + a(\mathbf{s}) \{\beta_a + w_a(\mathbf{s}) \} + b(\mathbf{s}) \{\beta_b + w_b(\mathbf{s}) \} + \epsilon(\mathbf{s}) \tag{11}$$

预测变量 $a(\mathbf{s})$ 和 $b(\mathbf{s})$ 取自均值为零、方差为一的独立正态分布。回归系数 $\beta_0$、$\beta_a$ 和 $\beta_b$ 分别等于 $1$、$10$ 和 $-10$。与截距和预测变量相关的 $r=3$ 空间随机效应是从不可分离的多变量 GP 中生成的。用于构造多元 GP 的 $nr × nr$ 协方差矩阵中第 $(i, j)$ 个 $r × r$ 块的互协方差函数，即 $\mathbf{K}_{\boldsymbol{\theta}}(\mathbf{s}_i, \mathbf{s}_j) = \mathbf{A} \boldsymbol{\Gamma} (\mathbf{s}_i, \mathbf{s}_j) \mathbf{A}^T$，是

$$\left(\begin{array}{ccc} 1 & 0 & 0 \\ -1 & 1 & 0 \\ 0 & 1 & 0.1 \end{array}\right)\left(\begin{array}{ccc} \exp \left(-4 d_{i, j}\right) & 0 & 0 \\ 0 & \exp \left(-6 d_{i, j}\right) & 0 \\ 0 & 0 & \exp \left(-6 d_{i, j}\right) \end{array}\right)\left(\begin{array}{ccc} 1 & -1 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & 0.1 \end{array}\right)$$

其中 $d_{i,j}$ 是位置 $\mathbf{s}_i$ 和 $\mathbf{s}_j$ 之间的欧氏距离，$\boldsymbol{\Gamma}(\mathbf{s}_i, \mathbf{s}_j)$ 的对角线元素是 $k = 1, 2, \ldots,r$ 的指数相关函数 $\exp(−\phi_k d_{i,j})$。图 2(a)-(c) 显示了 $\mathbf{w}_0$ 、$\mathbf{w}_a$ 和 $\mathbf{w}_b$ 的实现。残差项 $\epsilon(\mathbf{s})$ 是根据均值为零且方差 ${\boldsymbol{\tau}}^2 = 0.1$ 的正态分布模拟的。

下面的代码指定了模型协方差参数的先验分布，以及 MCMC 采样器初始和 Metropolis 提案方差值。在这里，我们对每个支持从 $1$ 到 $10$ 的空间衰减参数使用均匀先验。交叉协方差矩阵的先验是具有自由度 $r$ 和恒等尺度矩阵的 $IW$。测量误差（或块金方差）的先验遵循形状为 $2$ 尺度为 $1$ 的 $IG$。

参数先验、初始值和 Metropolis 采样提议方差分别通过 priors、 starting 和 tuning 参数传递给 spSVC。提议模型是通过 formula 参数指定的，使用的语法类似于 R 中的 lm 中的语法，并添加了 svc.cols 参数，该参数接受整数索引或字符名称向量以指示空间变化的设计矩阵列（即 $\delta_j = 1$ 的 $\mathbf{X}$ 那些列)。例如，在下面对 spSVC 的调用中，传递给 svc.cols 的向量表示我们希望截距和标记为 a 和 b 的列遵循多元 GP（或等效地，可以使用参数值 c(1,2,3) )。

data(SVCMvData.dat)

r <- 3 
n.ltr <- r*(r+1)/2 
priors <- list("phi.Unif"=list(rep(1,r), rep(10,r)), 
               "K.IW"=list(r, diag(rep(1,r))), 
               "tau.sq.IG"=c(2, 1)) 
starting <- list("phi"=rep(3/0.5,r), "A"=c(1,0,0,1,0,1), "tau.sq"=1) 
tuning <- list("phi"=rep(0.1,r), "A"=rep(0.01, n.ltr), "tau.sq"=0.01) 
sim.m <- `spSVC`(y~a+b, coords=c("x.coords","y.coords"), data=SVCMvData.dat, 
               starting=starting, svc.cols=c("(Intercept)","a","b"), 
               tuning=tuning, priors=priors, cov.model="exponential", 
               n.samples=10000, n.report=5000, n.omp.threads=4)

模型输出：

---------------------------------------
General model description 
---------------------------------------
Model fit with 200 observations. 
Number of covariates 3. 
Number of space varying covariates 3. 
Using the exponential spatial correlation model. 
Number of MCMC samples 10000. 

Priors and hyperpriors: 
beta flat. 
K IW hyperpriors: 
df: 3.00000 
S: 
1.000 0.000 0.000 
0.000 1.000 0.000 
0.000 0.000 1.000 

phi Unif lower bound hyperpriors: 1.000 1.000 1.000 
phi Unif upper bound hyperpriors: 10.000 10.000 10.000 

tau.sq IG hyperpriors shape=2.00000 and scale=1.00000 

Source compiled with OpenMP, posterior sampling is using 4 thread(s). 

------------------------------------------------
Sampling
------------------------------------------------
Sampled: 5000 of 10000, 50.00% 
Report interval Metrop. 
Acceptance rate: 34.64% 
Overall Metrop. Acceptance rate: 34.64% 
------------------------------------------------
Sampled: 10000 of 10000, 100.00% 
Report interval Metrop. 
Acceptance rate: 34.24% 
Overall Metrop. Acceptance rate: 34.44% 
-------------------------------------------------

如 第 2 节 所述，spSVC 仅为模型协方差参数计算并返回 MCMC 样本。如果 verbose=TRUE，则将基本模型规格写入终端，然后更新采样器的进度和 Metropolis 算法接受率。使用 n.report 参数控制采样器进度报告间隔。应该调整 Metropolis 采样器提案方差以实现 $~30-50\%$ 之间的接受率（参见 Gelman、Carlin、Stern、Dunson、Vehtari 和 Rubin 2013，了解模型拟合最佳实践）。如果证明难以保持可接受的接受率，则可以添加 amcmc 参数以调用自适应 MCMC 算法（Roberts 和 Rosenthal 2009）自动调整调整以实现目标接受率（有关更多详细信息，请参见手册页).

上面 spSVC 调用中的 n.omp.threads 参数要求给定 MCMC 迭代中的循环键通过 openMP 使用 4 个线程（Dagum 和 Menon 1998）。如果用户的 R 设置为使用并行版本的 BLAS，则 n.omp.threads 还将控制某些 LAPACK 矩阵运算中的线程数。这种并行化可以大大减少采样器的运行时间。

用于进行此分析的计算机具有 Intel® Core™ i7-8550U CPU @ 1.80GHz 芯片，具有 4 个内核和使用 openMP 编译的 R，正如在调用 spSVC 后打印的 “General model description” 中所确认的那样，其中指出使用 OpenMP 编译的源代码，后验采样使用 4 个线程。如果 R 未使用 openMP 支持编译并且 n.omp.threads 设置为大于 1 的值，spSVC 将发出警告。除了 openMP 支持之外，R 的当前实现使用 openBLAS (Zhang 2016)，它是BLAS 能够利用多个处理器。图 1 显示了在可用 CPU 数量上完成 $10000$ 次 MCMC 迭代所需的运行时间

图 1：最佳运行时间的模拟数据分析评估。

在执行 spSVC 之后，sim.m 对象保存协方差参数 (p.theta.samples) 的 MCMC 样本以及数据和模型拟合细节。 spRecover 函数可能使用老化和变薄后的 p.theta.samples 进行组合采样，以从回归系数 $\boldsymbol{\beta}$ ( p.beta.recover.samples)、空间随机效应 $\mathbf{w}$ ( p.w.recover.samples) 和模型拟合值 ( p.y.samples )。 spRecover 还返回组合采样中使用的 p.theta.samples ( p.theta.recover.samples) 的子集。此外，为方便起见，spRecover 返回空间变化回归系数 $\tilde{\beta}(\mathbf{s})$ 的样本。 spRecover 将这些不同的组合采样输出附加到 spSVC 输入对象，即下面的 spRecover 返回的 sim.m 对象与 spSVC 返回的 sim.m 对象相同，只是添加了组合采样结果。除了为所有模型参数提供后验样本外，还需要分别通过 spPredict 和 spDiag 调用 spRecover 进行后续预测和模型拟合诊断。与 spSVC 一样，spRecover 在可用时通过 openMP 利用多个 CPU。

sim.m <- spRecover(sim.m, start=5000, thin=2, n.omp.threads=4, verbose=FALSE)

spSVC 和 spRecover 将样本作为 coda 对象返回，以简化后验摘要。下面的输出提供了老化后和细化后的中位数，以及 $\beta$ 的 $95\%$ 可信上下限、用于构造 $\mathbf{K}_{\boldsymbol{\theta}}$ 的交叉协方差矩阵、空间衰减 $\phi$ 和 ${\boldsymbol{\tau}}^2$。这些摘要表明模型很好地捕获了所使用的参数值来模拟数据。

round(summary(sim.m$p.beta.recover.samples)$quantiles[,c(3,1,5)],2) 

                50%    2.5%    97.5% 
(Intercept)    0.20   -1.05     1.16 
a             10.53    9.37    11.92 
b            -10.20  -11.08    -9.31

请注意，按照 第 2 节 中的符号，用于模拟数据的交叉协方差矩阵是

$$\mathbf{K}_\theta=\mathbf{A A}^{\top}\left(\begin{array}{ccc} 1 & 0 & 0 \\ -1 & 1 & 0 \\ 0 & 1 & 0.1 \end{array}\right)\left(\begin{array}{ccc} 1 & -1 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & 0.1 \end{array}\right)=\left(\begin{array}{ccc} 1 & -1 & 0 \\ -1 & 2 & 1 \\ 0 & 1 & 1.01 \end{array}\right)$$

协方差参数的后验总结如下。图 2(d)-(f) 给出了观测到的与估计的随机效应。

round(summary(sim.m$p.theta.recover.samples)$quantiles[,c(3,1,5)],2) 

                             50%    2.5%     97.5% 
K[1,1]                      1.11    0.58      2.78 
K[2,1]                     -1.00   -2.49     -0.46 
K[3,1]                     -0.06   -0.52      0.35 
K[2,2]                      1.89    1.17      3.60 
K[3,2]                      0.90    0.35      1.76 
K[3,3]                      1.14    0.63      2.03 
tau.sq                      0.18    0.10      0.33 
phi.(Intercept)             4.16    1.57      8.65 
phi.a                       4.98    2.09      9.81 
phi.b                       6.13    1.55      9.89

图 2：模拟数据分析、观测（真实）和模型估计（拟合）随机效应值。图 (d)-(f) 后中位数和 95% 可信区间分别显示为点和条。

3.2 空气污染数据分析

暂略。

4 总结

新的 spSVC 函数更全面地实现了 Finley 等人详述的计算高效的 MCMC 算法。 (2015) 并提供了一个灵活的软件工具来拟合空间变系数模型。虽然其他软件（其中一些在第 1 节中提到）提供了类似的空间自适应回归，但很少有软件同时提供单变量和多变量 GP 规范以及建议的采样算法和 OpenMP 并行化的使用以及对多变量的可选调用提供的计算效率。较低级别的 BLAS 和 LAPACK 多线程矩阵代数库。 此功能可适应非高斯和多变量结果，以及位置数量无法使用全秩空间 GP 的设置。

参考文献

Bakar KS, Kokic P, Jin H (2015a). “Hierarchical spatially varying coefficient and temporal dynamic process models using spTDyn.” Journal of Statistical Computation and Simulation. URL 10.1080/00949655.2015.1038267.

Bakar KS, Kokic P, Jin H (2015b). “A spatio-dynamic model for assessing frost risk in southeastern Australia.” Journal of the Royal Statistical Society, Series C. URL 10.1111/rssc. 12103.

Bakar KS, Kokic P, Jin H (2017). Spatially varying and spatio-temporal dynamic linear models. R package version 2.0.

Bakar KS, Sahu SK (2018). Spatio-Temporal Bayesian Modeling. R package version 3.3.

Bakka H, Rue H, Fuglstad GA, Riebler AI, Bolin D, Illian J, Krainski E, Simpson DP, Lindgren FK (2018). “Spatial modelling with INLA: A review.” ArXiv e-prints. 1802. 06350.

Banerjee S, Carlin BP, Gelfand AE (2014). Hierarchical Modeling and Analysis for Spatial Data. Second edition. Chapman & Hall/CRC, Boca Raton, FL.

Bivand R (2019). CRAN Task View: Analysis of Spatial Data. 2019-02-25, URL https: //cran.r-project.org/web/views/Spatial.html.

Bivand R, Yu D (2017). spgwr: Geographically Weighted Regression. R package version 0.6-32, URL https://CRAN.R-project.org/package=spgwr.

Brunekreef B, Holgate ST (2002). “Air Pollution and Health.” The Lancet, 360(9341), 12331242.

Bürgin R, Ritschard G (2017). “Coefficient-Wise Tree-Based Varying Coefficient Regression with vcrpart.” Journal of Statistical Software, Articles, 80(6), 1–33. ISSN 1548-7660. doi:10.18637/jss.v080.i06.

Candiani G, Carnevale C, Finzi G, Pisoni E, Volta M (2013). “A Comparison of Reanalysis Techniques: Applying Optimal Interpolation and Ensemble Kalman Filtering to Improve Air Quality Monitoring at Mesoscale.” Science of the Total Environment, 458-460(0), 7–14.

Carpenter B, Gelman A, Hoffman M, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, Riddell A (2017). “Stan: A Probabilistic Programming Language.” Journal of Statistical Software, Articles, 76(1), 1–32. ISSN 1548-7660. doi:10.18637/jss.v076.i01. URL https://www.jstatsoft.org/v076/i01.

Dagum L, Menon R (1998). “OpenMP: an industry standard API for shared-memory programming.” Computational Science & Engineering, IEEE, 5(1), 46–55.

Datta A, Banerjee S, Finley A, Hamm N, Schaap M (2016). “Nonseparable dynamic nearest neighbor Gaussian process models for large spatio-temporal data with an application to particulate matter analysis.” Annals of Applied Statistics, 10(3), 1286–1316. ISSN 19326157.

Denby B, Schaap M, Segers A, Builtjes P, Horalek J (2008). “Comparison of Two Data Assimilation Methods for Assessing PM10 Exceedances on the European Scale.” Atmospheric Environment, 42(30), 7122–7134. European Commission (2015). “European Union Air Quality Standards.” http://ec.europa.eu/environment/air/quality/standards.htm.

Fan J, Zhang W (2008). “Statistical Methods with Varying Coefficient Models.” Statistics and its interface, 1 1, 179–195.

Finley A, Banerjee S, Gelfand A (2015). “spBayes for Large Univariate and Multivariate Point-Referenced Spatio-Temporal Data Models.” Journal of Statistical Software, Articles, 63(13), 1–28. ISSN 1548-7660.

Fotheringham A, Brunsdon C, Charlton M (2002). Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Wiley. ISBN 9780471496168.

Gelfand AE, Ghosh SK (1998). “Model choice: A minimum posterior predictive loss approach.” Biometrika, 85(1), 1–11.

Gelfand AE, Kim HJ, Sirmans CF, Banerjee S (2003b). “Spatial Modeling With Spatially Varying Coefficient Processes.” Journal of the American Statistical Association, 98(462), 387–396.

Gelman A, Carlin J, Stern H, Dunson D, Vehtari A, Rubin D (2013). Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science. Taylor & Francis. ISBN 9781439840955.

Hamm N, Finley A, Schaap M, Stein A (2015). “A Spatially Varying Coefficient Model for Mapping PM10 Air Quality at the European scale.” Atmospheric Environment, 102, 393–405.

Hastie T, Tibshirani R (1993). “Varying-Coefficient Models.” Journal of the Royal Statistical Society. Series B (Methodological), 55(4), 757–796.

Hayfield T, Racine JS (2008). “Nonparametric Econometrics: The np Package.” Journal of Statistical Software, 27(5). URL http://www.jstatsoft.org/v27/i05/.

Heim S (2007). svcm: 2d and 3d space-varying coefficient models in R. R package version 0.1.2.

Helske J (2019). walker: Bayesian Regression with Time-Varying Coefficients. R package version 0.2.4-1, URL http://github.com/helske/walker.

Henderson HV, Searle SR (1981). “On deriving the inverse of a sum of matrices.” SIAM Review, 23(1), 53–60.

Hoek G, Krishnan RM, Beelen R, Peters A, Ostro B, Brunekreef B, Kaufman JD (2013). “Long-term air pollution exposure and cardio- respiratory mortality: a review.” Environmental Health, 12, 43.

Hothorn T, Buehlmann P, Kneib T, Schmid M, Hofner B (2018). mboost: Model-Based Boosting. R package version 2.9-1, URL https://CRAN.R-project.org/package=mboost.

Lindgren F, Rue H (2015). “Bayesian Spatial Modelling with R-INLA.” Journal of Statistical Software, 63(19), 1–25. URL http://www.jstatsoft.org/v63/i19/.

Loomis D, Grosse Y, Lauby-Secretan B, El Ghissassi F, Bouvard V, Benbrahim-Tallaa L, Guha N, Baan R, Mattock H, Straif S (2013). “The Carcinogenicity of Outdoor Air Pollution.” The Lancet Oncology, 14(13), 1262–1263.

R Core Team (2018). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Robert C, Casella G (2004). Monte Carlo Statistical Methods. Springer Texts in Statistics, second edition. Springer-Verlag.

Roberts GO, Rosenthal JS (2009). “Examples of Adaptive MCMC.” Journal of Computational and Graphical Statistics, 18(2), 349–367.

Rue H, Martino S, Chopin N (2009). “Approximate Bayesian Inference for Latent Gaussian Models Using Integrated Nested Laplace Approximations (with discussion).” Journal of the Royal Statistical Society B, 71, 319–392.

Schaap M, Timmermans RMA, Roemer M, Boersen GAC, Builtjes P, Sauter F, Velders G, Beck J (2008). “The LOTOS-EUROS Model: Description, Validation and Latest Developments.” International Journal of Environment and Pollution, 32(2), 270–290.

Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A (2001). “Bayesian Measures of Model Complexity and Fit.”

Stan Development Team (2018). “RStan: the R interface to Stan.” R package version 2.18.2, URL http://mc-stan.org/.

van de Kassteele J, Stein A (2006). “A Model for External Drift Kriging with Uncertain Covariates applied to Air Quality Measurements and Dispersion Model Output.” Environmetrics, 17(4), 309–322.

Vihola M, Helske J, Franks J (2017). “Importance Sampling Type Estimators Based on Approximate Marginal MCMC.” ArXiv e-prints. 1609.02541.

Wood S (2017). Generalized Additive Models: An Introduction with R. 2 edition. Chapman and Hall/CRC.

Zhang X (2016). “An Optimized BLAS Library Based on GotoBLAS2.” https://github. com/xianyi/OpenBLAS/. Accessed 2015-06-01.