基于似然的的方法

【摘 要】在空间统计学中，点参考数据模型通常采用高斯过程（场）建模，而超参数的推断则主要有基于经验的矩量估计法和基于似然的统计推断方法。本文主要介绍基于似然的统计推断方法。文中涉及最大似然估计、受限最大似然估计、组合似然近似估计、渐进特性分析等内容。本文内容摘自 Gelfand 的 《空间统计手册》第四章。

【原 文】 Gelfand, A.E. et al. (2010)， Handbook of spatial statistics (chapter 4). CRC press.

• 关于似然方法的基础资料见：
  • 参见 普度大学机器人视觉实验室的 《最大似然、最大后验与贝叶斯方法的区别》
  • 参见 Reid 等 2010 年的 《似然函数与基于似然的推断》
  • 参见 Cousineau 等 2016 年 《似然概念的扫盲帖》
• 对非似然方法感兴趣的同学，可以阅读：
  • 参见 《近似贝叶斯计算（ABC）索引贴》

4.1 概述

上一章考虑了结合使用矩量法和最小二乘法来估计地质统计模型的参数（参见 《点参考数据及克里金法》 ）。这些方法统称为 “经典地质统计学” ，通常不明确要求任何分布假设，但在任何已知意义上都不是最优的。

在本章中，我们将介绍通过基于似然的方法估计模型参数的地质统计模型。这些方法遵循似然原理，并且在适当规律性条件下，可能预期具有某些最优性能。

在这里，我们暂时仅考虑单空间变量 $Y$ 的情形，定义在（原则上）有界的研究区域 $D \subset \mathbb{R}^d$ 中的所有点上，其中 $d = 2$ 或 $3$ 。我们将其建模为一个随机过程 $Y(·) \equiv \{ Y(\mathbf{s}):\mathbf{s} \in D \}$ 。

假设在 $D$ 的 $n$ 个不同点处都观测到了各自的 $Y(·)$ （可能有误差），那么我们的目标是从这些观测中推断出此随机过程。

基于似然进行推断的方法已经在高斯随机场中得到全面发展。高斯随机场是任意有限维分布均为多元高斯分布的一类随机过程。因此，本章大部分内容（ 第 4.2 ~ 4.7 节 ）仅考虑此类过程的推断问题。在第 4.8 节中，我们简要考虑了非高斯过程的基于似然推断方法。本章全文都强调频率派方法，尽管也提到了一些贝叶斯程序。如果希望从贝叶斯角度更全面地理解和处理此问题，请参阅 第 7 章 。

4.2 最大似然估计

4.2.1 符号与似然的定义

考虑 式（3.1） 的地统计模型，假设均值函数是线性函数，误差过程服从高斯分布。也就是说，假设

$$Y(\mathbf{s}) = X(\mathbf{s})^{\top} \boldsymbol{\beta} + e(\mathbf{s})$$

其中 $X(\mathbf{s})$ 是 $\mathbf{s}$ 处的 $p$ 维可观测协变量向量，$\boldsymbol{\beta}$ 是 $p$ 维未知参数向量，$e(· ) \equiv \{ e(\mathbf{s}):\mathbf{s} \in D \}$ 为零均值高斯过程。

进一步假设 $e(·)$ 具有协方差函数

$$\operatorname{Cov}[e(\mathbf{s}),e(\mathbf{t})] \equiv C(\mathbf{s,t}; \boldsymbol{\theta})$$

其中 $\boldsymbol{\theta}$ 是给定参数空间中的 $m$ 维未知超参数向量，且 $\Theta \subset \mathbb{R}^m$。根据模型定义，协方差函数必须是正定的（即其对应的协方差矩阵必须是非负定的）。$\boldsymbol{\beta}$ 和 $\boldsymbol{\theta}$ 构成的联合参数空间通常被认为是 $\mathbb{R}^p$ 和 $\Theta$ 的笛卡尔积。请注意，此公式并未要求 $e(·)$ 在任何意义上都是平稳的。这与经典地统计学形成了对比，后者需要某种形式的平稳性假设（二阶或本征）。

采用如下符号定义：

• 将在位置 $\mathbf{s}_1,\ldots, \mathbf{s}_n$ 处观测到的过程 $Y(·)$ 的值表示为 $n$ 维向量 $\mathbf{Y} = [Y(\mathbf{s}_1),Y(\mathbf{s}_2) ,\ldots ,Y(\mathbf{s}_n)]^{\top}$
• 将观测到的协变量 $X(·)$ 的值表示为 $n× p$ 矩阵 $\mathbf{X} = [\mathbf{X}(\mathbf{s}_1), \mathbf{X}(\mathbf{s}_2),\ldots , \mathbf{X}(\mathbf{s}_n)]^{\top}$ ，并假设 $n > p$ 且 $\operatorname{rank}(\mathbf{X}) = p$
• 令 $\boldsymbol{\Sigma} = \boldsymbol{\Sigma}(\boldsymbol{\theta})$ 表示 $n × n$ 协方差矩阵 $\operatorname{Var}(\mathbf{Y})$ ，其第 $(i, j)$ 个元素为 $C(\mathbf{s}_i , \mathbf{s}_j ; \boldsymbol{\theta})$

因为 $Y(·)$ 是高斯的，所以 $\mathbf{Y}$ 的联合分布是多元高斯，其均值向量为 $\mathbf{X}\boldsymbol{\beta}$，其协方差矩阵为 $\boldsymbol{\Sigma} (\boldsymbol{\theta})$ 。其似然函数为：

$$l(\boldsymbol{\beta, \theta}; \mathbf{Y}) = (2 \pi)^{-n/2}|\boldsymbol{\Sigma (\theta)}|^{-1/2} \exp\{ -(\mathbf{Y - X}\boldsymbol{\beta})^{\top} \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta})( \mathbf{Y - X}\boldsymbol{\beta})/2 \}$$

注：$\mathbf{Y}$ 仅仅被视为随机过程的一次实现，这与频率派多样本实验的方法论完全不同。这里之所以被称为似然，正是因为它只是建模了一个样本的可能性。虽然此似然在形式上与频率派的边缘似然非常相似，但概念内涵截然不同

$\boldsymbol{\beta}$ 和 $\boldsymbol{\theta}$ 的最大似然 (Maximum Likelihood, ML) 估计被定义为能够使 $l(\boldsymbol{\beta, \theta}; \mathbf{Y})$ 最大化的任意参数值（分别在 $\mathbb{R}^p$ 和 $\Theta$ 中）。在实际工作中，通常会采用对数形式，上述参数值也会同时最大化如下的对数似然函数。

$$\log l(\boldsymbol{\beta}, \boldsymbol{\theta}; \mathbf{Y}) =- \frac{n}{2} \log(2\pi ) - \frac{1}{2} \log |\boldsymbol{\Sigma}(\boldsymbol{\theta})| - \frac{1}{2}(\mathbf{Y - X}\boldsymbol{\beta})^{\top} \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta})(\mathbf{Y - X}\boldsymbol{\beta})$$

4.2.2 最大似然估计的实施

剖面似然（Profile likelihood）： 当进行最大似然估计时，有时为了计算可行，会假设部分参数为固定的已知数，然后只针对剩余部分参数做最大似然估计，此时的似然函数被称为剖面似然。例如，假设 $\boldsymbol{\theta}= \boldsymbol{\theta}_0$ 已知，则上式的剖面似然为 $l(\boldsymbol{\beta}, \boldsymbol{\theta}_0; \mathbf{Y})$，其中 $\boldsymbol{\beta}$ 为待优化的参数。

（1）基于剖面似然的估计策略

如果假设 $\boldsymbol{\theta} \in \Theta$ 固定（ 如 $\boldsymbol{\theta}_0$，即假设 $\boldsymbol{\theta}$ 已知 ），则最大化 $\log l(\boldsymbol{\beta}, \boldsymbol{\theta}_0; \mathbf{Y})$ 的参数 $\boldsymbol{\beta}$ 值由下式给出（推导略）：

$$\hat{\boldsymbol{\beta}} = (\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta}_0) \mathbf{X})^{-1}\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}( \boldsymbol{\theta}_0)\mathbf{Y}$$

如果假设协方差 $\operatorname{Var}(\mathbf{Y})$ 就是 $\boldsymbol{\Sigma} (\boldsymbol{\theta}_0)$ ，则上式得到的估计 $\hat{\boldsymbol{\beta} }$ 是 $\boldsymbol{\beta}$ 的广义最小二乘估计。基于 $\hat{\boldsymbol{\beta} }$，我们可以进一步依据关于 $\boldsymbol{\theta}$ 剖面似然来估计 $\hat{\boldsymbol{\theta}}$。

也就是说，可以这样定义 $\boldsymbol{\theta}$ 和 $\boldsymbol{\beta}$ 的最大似然估计 $\hat{\boldsymbol{\theta}}$ 和 $\hat{\boldsymbol{\beta}}$：

（1）首先，找到 $\hat{\boldsymbol{\theta}}$。令 $\hat{\boldsymbol{\theta}}$ 是 $\Theta$ 中任意能够最大化如下似然的超参数值（推导略）：

$$L(\boldsymbol{\theta}; \mathbf{Y}) = -\frac{1}{2} \log |\boldsymbol{\Sigma} (\boldsymbol{\theta})| - \frac{1}{2} \mathbf{Y}^{\top} \mathbf{P}(\boldsymbol{\theta})\mathbf{Y} \tag{4.1}$$

式中，

$$\mathbf{P}(\boldsymbol{\theta}) = \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta}) - \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta}) \mathbf{X} (\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta})\mathbf{X})^{-1}\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta}) \tag{4.2}$$

（2）然后，将 $\hat{\boldsymbol{\theta}}$ 当作上述的 $\boldsymbol{\theta} _0$ 代入对数似然公式，得到 $\boldsymbol{\beta}$ 的估计：

$$\hat{ \boldsymbol{\beta}} = (\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}( \hat{ \boldsymbol{\theta}})\mathbf{X})^{-1}\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}(\hat{ \boldsymbol{\theta}})\mathbf{Y}$$

这种从对数似然函数中消去 $\boldsymbol{\beta}$ 的方法被称为 profiling（剖面化），而函数 $L(\boldsymbol{\theta}; \mathbf{Y})$ 被称为关于 $\boldsymbol{\theta}$ 的剖面对数似然函数（profile log likelihood function）。

（2）估计 $\hat{\boldsymbol{\theta}}$ 的方法

针对 $\boldsymbol{\theta} \in \Theta$ 最大化 $L(\boldsymbol{\theta}; \mathbf{Y})$ 是一个在 $\mathbb{R}^m$ 的子集 $\Theta$ 上的有约束非线性优化问题，其封闭形式解仅在特殊情况下存在。因此，ML 估计大多采用数值方法。

• 方法 1 : 贪婪搜索，即在在参数空间 $\Theta$ 中进行暴力搜索，这在 $\boldsymbol{\theta}$ 维度较低时比较有效。
• 方法 2 : 迭代算法，采用梯度下降等算法迭代逼近真实解。

在基于梯度的迭代算法中，第 $(k + 1)$ 次迭代的 $\boldsymbol{\theta}^{(k+1)}$ 是通过更新第 $k$ 次迭代 $\boldsymbol{\theta}^{(k)}$ 来计算的：

$$\boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} + \rho^{(k)}\mathbf{M}^{(k)}\mathbf{g}^{(k)}$$

式中的 $\rho^{(k)}$ 是一个标量，$\mathbf{M}^{(k)}$ 是一个 $m × m$ 矩阵，$\mathbf{g}^{(k)}$ 是 $L$ 在 $\boldsymbol{\theta} = \boldsymbol{\theta}^{(k)}$ 处的梯度，即 $\mathbf{g}^{(k)} = \frac{\partial L(\boldsymbol{\theta} ; \mathbf{Y})}{ \partial \boldsymbol{\theta}} |_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(k)}}$ 。$\mathbf{M}^{(k)}$ 和 $\mathbf{g}^{(k)}$ 的矩阵乘积可以被认为是定义搜索方向（相对于第 $k$ 次迭代 $\boldsymbol{\theta}^{(k)}$ ），而 $\rho^{(k)}$ 定义了要采取的步骤的大小方向。

通常与最大化对数似然结合使用的两种梯度算法是 Newton-Raphson 算法 和 Fisher 打分算法 ：

• Newton-Raphson 算法：$\mathbf{M}^{(k)} = (\mathbf{B}^{(k)})^{-1}$，其中 $\mathbf{B}^{(k)}$ 的第 $(i, j)$ 个元素为 $- \frac{\partial^2 L(\boldsymbol{\theta}; \mathbf{Y})}{ \partial \theta_i \partial \theta_j} |_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(k)}}$ 。
• Fisher 评分算法 ： $\mathbf{M}^{(k)} = (\mathbf{B}^{(k)})^{-1}$，其中 $\mathbf{B}^{(k)}$ 是与 $\boldsymbol{\theta}^{(k)}$ 处的 $L(\boldsymbol{\theta}; \mathbf{Y})$ 有关的 Fisher 信息矩阵，即 $\mathbf{B}^{(k)}$ 的第 $(i, j)$ 个元素为 $\mathbb{E} \{ -\frac{\partial^2 L(\boldsymbol{\theta}; \mathbf{Y})}{\partial \theta_i \partial \theta_j} |_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(k)}} \}$。

两种算法中均有 $\rho^{(k)} = 1$ 。因此，Fisher 评分与 Newton-Raphson 相同，只是二阶偏导数被期望所取代。二阶偏导数及其期望的表达式可以在 Mardia 和 Marshall (1984) 中找到。

还有一些非梯度算法。例如，Nelder–Mead 单纯形算法（Nelder 和 Mead，1965 ^[22]）也可能有效。

（3）算法的选择

需要权衡某些因素来决定使用哪一种算法，如：$\boldsymbol{\theta}$ 初始值的选择、收敛准则、协方差函数的参数化形式、包容 $\boldsymbol{\theta}$ 约束（通常是线性不等式约束）的方法等。 Harville (1977 ^[14]) 提供了一些关于此类实施问题的指南。 应该注意的是：初始值的一种可能选择是经典地质统计方法获得的 $\boldsymbol{\theta}$ 估计值。

在许多标准统计问题中，存在唯一的最大似然估计。但在目前情况下，并不能完全保证其存在性或唯一性，甚至不能保证似然函数的所有局部最大值都是全局最大值。实际上，Warnes 和 Ripley (1987 ^[30]) 以及 Handcock (1989 ^[13]) 表明，对应于具有球形协方差函数的平稳高斯过程的网格观测值的似然函数，通常会出现多个峰值，并且这些峰值之间存在很好地分割。 Rasmussen 和 Williams (2006 ^[25]) 展示了具有高斯协方差函数的平稳高斯过程在一条线上的七个不规则间隔位置观测到的双峰似然曲面。

不过，Handcock (1989 ^[13]) 和 Mardia 和 Watkins (1989 ^[21]) 的结果，加上作者的经验表明，对于 Matérn 类协方差函数（例如指数函数）和典型大小的数据集而言，似然函数出现多峰的现象在实践中极为罕见。

无论如何，我们都可以重复从散布比较广泛的多个初值开始进行迭代，以得到更为可靠的结果。

4.3 受限最大似然估计（REML）

此处可参考 Daniel 等 2017年的 《最大似然法与受限最大似然法的比较》 。

（1）最大似然估计的问题

尽管空间协方差参数 $\boldsymbol{\theta}$ 的 ML 估计有理想特性，但它们有一个众所周知的缺点： 估计 $\boldsymbol{\beta}$ 时的 “自由度损失” 会带来偏差（Harville，1977 ^[14]）。即使对于中等大小的样本，在空间相关性很强或 $\boldsymbol{\beta}$ 的维数 $p$ 很大时，此偏差也会很大。

可以通过一种被称为 受限最大似然 (REML) 估计 的变体来显著减少偏差，在某些特殊情况下甚至可以完全消除偏差。

REML 最初由 Patterson 和 Thompson (1971 ^[23], 1974 ^[24]) 提出，并被用于方差分量模型；而 Kitanidis (1983 ^[17]) 则第一个将其用于空间模型。这两种模型现在和最大似然估计一样流行，甚至更流行。

（2）基本思想：误差对照

在 REML 估计中，不直接最大化与观测值相关的似然函数（或等效的对数似然函数），而是 最大化与误差对照相关的似然函数。

误差对照 是观测值的线性组合（ 即 $\mathbf{a^{\top}Y}$ ），其关于所有 $\boldsymbol{\beta}$ 和所有 $\boldsymbol{\theta} \in \Theta$ 的期望为零。

误差对照的性质 ：当向量 $\mathbf{a}$ 和向量 $\mathbf{b}$ 之间线性独立时，两个误差对照 $\mathbf{a^{\top}Y}$ 和 $\mathbf{b^{\top}Y}$ 之间也线性独立。

可以看出，$[\mathbf{I} - \mathbf{X}(\mathbf{X}^{\top} \mathbf{X})^{-1}\mathbf{X}^{\top}]\mathbf{Y}$ 中的任何一组 $n - p$ 个线性独立元素都可以被用于构造误差对照，实际上使用哪组误差对照并没有区别，只要其数量为 $n - p$ 并且独立即可。因为与任何此类集合相关的对数似然函数，与如下受限似然最多相差一个加性常数，且该常数与 $\boldsymbol{\beta}$ 或 $\boldsymbol{\theta}$ 无关：

$$L_{R}(\boldsymbol{\theta}; \mathbf{Y}) =- \frac{1}{2} \log |\boldsymbol{\Sigma} (\boldsymbol{\theta})| - \frac{1}{2} \log |\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta}) \mathbf{X}|- \frac{1}{2} \mathbf{Y}^{\top} \mathbf{P}(\boldsymbol{\theta})\mathbf{Y}$$

上面的受限似然公式中， $\mathbf{P}(\boldsymbol{\theta})$ 由 式（4.2） 定义。可以看到： 受限似然 $L_{R}(\boldsymbol{\theta}; \mathbf{Y})$ 与 式（4.1） 的剖面对数似然 $L(\boldsymbol{\theta}; \mathbf{Y})$ 之间仅相差一个 $- \frac{1}{2} \log |\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta}) \mathbf{X}|$。

$\boldsymbol{\theta}$ 的 REML 估计是能够使受限似然 $L_{R}$ 达到最大值的任意 $\tilde{\boldsymbol{\theta}} \in \boldsymbol{\theta}$ 。该估计可以采用与最大似然估计相同的数值算法。

一旦获得了 $\boldsymbol{\theta}$ 的 REML 估计，就可以得到 $\boldsymbol{\beta}$ 的相应估计，作为其在 $\boldsymbol{\theta} = \tilde{\boldsymbol{\theta}}$ 处计算的广义最小二乘估计，即

$$\tilde{\boldsymbol{\beta}} = (\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1}(\tilde{\boldsymbol{\theta}}) \mathbf{X})^{-1}\mathbf{X}^{\top} \boldsymbol{\Sigma}^{-1} (\tilde{\boldsymbol{\theta}})\mathbf{Y}$$

（3）争论与讨论

REML 是否能有效减少最大似然估计引起的偏差呢？ 两项已发表的模拟研究比较了 REML 和 ML 估计器在空间模型中的性能。

• Zimmerman 和 Zimmerman (1991 ^[40]) 对具有恒定均值和指数协方差函数且没有块的平稳高斯随机场进行的一项早期研究，他们表明，在此情况下，基台和变程参数的 REML 估计量确实比其 ML 方法有更少的偏差，但是基台的 REML 估计值的分布存在较重的上尾。

• 在最近一项更彻底的研究中，Irvine、Gitelman 和 Hoeting（2007 年 ^[15]）对具有常量均值和 exponential-plus-nugget 协方差函数的平稳高斯随机场进行了估计性能比较，后者具有几个不同的有效变程和块金基台比。他们给出了大体相似结果，但对 REML 更为不利；例如，在许多情况下，有效变程参数的 REML 估计的经验分布存在非常重的上尾，以至于该其估计均方误差比 ML 估计大得多。随着有效变程或块金基台比的增加，上尾问题变得更糟。

因此，尽管 REML 似乎确实有效地减少了偏差，但也可能增加估计方差，除非空间依赖性的持续性（由有效变程衡量）很弱。选择 ML 估计还是 REML 估计主要取决于人们认为哪个问题更严重，严重低估或严重高估有效变程和基台参数，会分别对均值参数和预测产生过度乐观或保守的推论。

值得注意的是，这里提到的两项研究都只考虑了均值恒定的随机过程，即 $p = 1$ 。Zimmerman (1986 ^[37]) 证明了 REML 相对于 ML 估计的无块金指数协方差函数的相对性能提高随着 $\boldsymbol{\beta}$ 维数的增加会大大增加，考虑到 REML 存在的理由，这应该不足为奇。尽管如此，在这位作者看来，仍然需要进一步研究。

不过， REML 估计确实有一个明显的优势 。它可用于估计某些非平稳过程的参数，这些过程称为 $k$ 阶高斯固有随机场 (IRF-ks)。此类过程的特征是某些称为广义增量的线性组合的概率结构；更具体地说，高斯 IRF-k 的 $k$ 阶广义增量是均值为零的联合高斯函数和广义协方差函数。广义协方差函数不指定观测值本身的协方差结构，而仅指定广义增量的协方差结构。本质平稳的随机场和 IRF-0 只是同一过程的不同名称；此外，IRF-0 的广义增量和广义协方差函数分别仅为 $\{ Y(\mathbf{s}) - Y(\mathbf{t}):\mathbf{s}, t \in D \}$ 和半变异函数。事实证明，任何高斯 IRF-k 的广义增量与具有 $k$ 阶多项式平均结构的高斯随机场的误差对照一致，因此 REML 估计可以针对此类过程进行。然而，由于观测缺乏完全指定的协方差结构，ML 不能。例如，REML 而不是 ML 可用于估计 IRF-0 的功率半变异函数。这种情况下的 REML 对数似然是

$$L_{R}(\boldsymbol{\theta}; \mathbf{Y}) =- \frac{1}{2} \log |\mathbf{A^{\top} K}(\boldsymbol{\theta}) \mathbf{A}|- \frac{1}{2} \mathbf{Y}^{\top} \mathbf{A}(\mathbf{A^{\top} K}(\boldsymbol{\theta}) \mathbf{A})^{-1} \mathbf{A^{\top} Y}$$

其中 $\mathbf{A}$ 是 $(n - 1) × n$ 矩阵，其第 $i$ 行有第 $i$ 个元素 $1$ ，第 $(i + 1)$ 个元素为 $-1$ ，其他地方为零，$\mathbf{K}(\boldsymbol{\theta})$ 是一个 $n \times n$ 矩阵， 其第 $(i, j )$ 个元素为 $\gamma (\mathbf{s}_i - \mathbf{s}_j ; \boldsymbol{\theta}) = \theta_1 \|\mathbf{s}_i -\mathbf{s}_j \|^{\theta_2}$。

4.4 渐近结果

为了推断高斯随机场的参数，例如，估计标准差和构建置信区间，了解基于似然的估计的渐近特性非常有益。但与大多数其他统计领域相比，对于连续空间过程的统计，有两个（至少）完全不同的渐近框架可以合理地吸引：

• （1） 递增域渐近 （increasing domain asymptotics）：在这种情形中，数据位置之间的距离下界为零，而观测的空间域是无界的，
• （2） 填充渐近 （infill，也称为固定域渐进 fixed domain asymptotics），在此情形下，会在固定且有界的域中进行更密集的观测。

在第 6 章 《空间过程的渐近》 中，Michael Stein 在两个框架内相当详细地给出了一些与空间预测和估计相关的已知渐近结果。下面，我们简要总结一些与估计有关的结论：

• 在递增域渐近框架中，高斯随机场参数的 ML 和 REML 估计是一致且渐近高斯的（Mardia 和 Marshall，1984 ^[20]；Cressie 和 Lahiri，1996 ^[6]）。

• 在填充渐近框架中，可用结果相当有限。而且在此框架中 ML 估计的渐近表现可能有很大不同。例如，对于协方差函数为无块金-指数函数的一维零均值高斯随机场，基台参数和变程参数的 ML 估计存在在填充渐近不一致现象（尽管其乘积与参数的乘积可能渐进一致) (Ying, 1991 ^[33])。Zhang 和 Zimmerman（2005 年 ^[36]）发现， 对于该过程与另一个增加了块金效应的类似过程，在两种渐进框架下 ML 估计都一致的那些参数，其 ML 估计的有限样本近似分布表现也大致相同，但对于估计不一致的那些参数，填充渐近框架提供的有限样本近似表现更好。

因此，建议在可用的情况下采用基于填充渐近结果进行推断。

4.5 假设检验和模型比较

（1）假设检验

对模型参数（ $\boldsymbol{\beta}$ 和 $\boldsymbol{\theta}$ ）进行假设检验，可以被视为 “完整” 模型和 “简化” 模型的对比（对应于替代假设和零假设），并且可以通过似然比来进行检验。

这相当于将 “两个模型在 ML 估计处计算的两倍对数似然函数差” 与 “自由度为 $q$ 的卡方分布的百分位数” 进行比较，其中 $q$ 是两个模型参数空间的维数差。当然，此过程基于以下假设：在零假设下，似然比这个统计量的两倍负对数是渐近分布的随机变量，服从自由度为 $q$ 的卡方分布。在增加域渐近和适当正则性条件下，该假设是合理的。

（2）信息准则

非嵌套模型还可以在似然框架内通过惩罚似然准则进行非正式比较，例如 Akaike 的信息准则：

$$AI C =-2 \log l( \hat{\boldsymbol{\beta}}, \hat{ \boldsymbol{\theta}}) + 2( p + m)$$

或 Schwarz 的贝叶斯信息准则：

$$BIC =-2 \log l(\hat{ \boldsymbol{\beta}}, \hat{ \boldsymbol{\theta}}) + ( p + m) \log n$$

在比较两个模型时，准则值较小的模型被判断为两个模型中较好的那个模型。这两个准则在模型拟合和模型简约性之间做了权衡，只是方式略有不同，其结果是 $BIC$ 比 $AIC$ 更喜欢小的模型。观察到：可以比较具有不同均值函数、不同协方差函数或两者的模型，尽管尚不清楚平等地惩罚均值参数和协方差参数是否合理。

在 REML 框架内，关于 $\boldsymbol{\theta}$ 的假设可以通过比较在 REML 估计中计算的受限对数似然来进行类似检验。为了比较非嵌套模型，可以再次使用惩罚性受限似然准则；在这种情况下，

$$AIC =-2L_{R}(\tilde{\boldsymbol{\theta}}) + 2m \quad \text{and} \quad BIC =-2L_{R}(\tilde{\boldsymbol{\theta}}) + m \log n$$

非常重要的是：此类比较仅在协方差结构不同的模型之间适用。在 REML 框架内，无法对具有不同均值结构的模型进行有效比较，因为两个此类模型的 误差对照 不同。

在贝叶斯框架内进行模型比较的另外一个流行准则是偏差信息准则 (DIC)，

$$DIC = \overline{ D} + p_D$$

（Spiegelhalter、Best、Carlin 和 van der Linde，2002 年 ^[26]）。其中 $\overline{D} = \mathbb{E} [D(\boldsymbol{\beta}, \boldsymbol{\theta})| \mathbf{Y}]$ ，其中 $D(\boldsymbol{\beta}, \boldsymbol{\theta}) =-2 \log l(\boldsymbol{\beta}, \boldsymbol{\theta}; \mathbf{Y}) + c$，其中 $c$ 是一个可以忽略的常数，因为它跨模型相同，$p_D = \overline{D} - D(\overline{ \boldsymbol{\beta}}, \overline{ \boldsymbol{\theta}})$， 其中 $(\overline{ \boldsymbol{\beta}}, \overline{ \boldsymbol{\theta}}) = \mathbb{E} [(\boldsymbol{\beta}, \boldsymbol{\theta})|\mathbf{Y}]$。$\overline{D}$ 是模型拟合度的度量，而 $p_D$ （称为有效参数数）是模型复杂度的度量， $DIC$ 平衡了这两者。 $DIC$ 可以很容易地从马尔可夫链蒙特卡洛 (MCMC) 模拟生成的样本中计算出来，只需通过对样本平均 $D(\boldsymbol{\beta}, \boldsymbol{\theta})$ 和 $(\boldsymbol{\beta}, \boldsymbol{\theta})$ 分别获得 $\overline{D}$ 和 $(\overline{ \boldsymbol{\beta}}, \overline{ \boldsymbol{\theta}})$ ，并在后一个量上计算 $D$ 。

4.6 计算方面的问题

基于似然的估计程序需要大量计算。这主要源于需要重复生成和求逆协方差矩阵，并计算其行列式。无论采用迭代算法更新参数估计，还是进行网格搜索，此类计算负担在实际工作中对于基于似然方法都是一种严重障碍。幸运的是，在某些情况下可以显著减少基于精确似然推断所需要的计算量。

（1） 单一尺度参数的情况

当尺度参数是协方差函数的一个参数时，可以天然减少计算负担。此时只需将协方差矩阵表示为 “尺度参数和另一个与之无关的矩阵之积” ，即对于某些矩阵 $\mathbf{W}(\boldsymbol{\theta}_{-1})$，其中 $\boldsymbol{\theta}_{-1} = (\theta_2, \theta_3,\ldots , \theta_m)^{\top} \in \boldsymbol{\Theta}_{-1}$，协方差矩阵可以被表示为 $\boldsymbol{\Sigma} (\boldsymbol{\theta}) = \theta_1 \mathbf{W}(\boldsymbol{\theta}_{-1})$。这里， $\boldsymbol{\Theta}_{-1}$ 被定义为 $\Theta =\{ \theta_1 > 0, \boldsymbol{\theta}_{-1} \in Θ-1 \}$。

这种情况的一种案例是当随机场 $Y(·)$ 平稳且观测 $\mathbf{Y}$ 同方差的情况，此时 $\boldsymbol{\Sigma} (\boldsymbol{\theta})$ 可以采用上述方式表示，且 $theta_1$ 为观测的公共方差（或半变异函数基台）参数，而 $W(\boldsymbol{\theta}_{-1})$ 为相关矩阵。

在任何此类情况下，很容易证明 $L(\boldsymbol{\theta}; \mathbf{Y})$ 在 $(\hat{ \theta}_1, \hat{ \boldsymbol{\theta}}^{\top}_{-1})^{\top}$ 处取最大值，其中

$$\hat{\theta}_1 = \mathbf{Y}^{\top} \mathbf{Q}(\hat{ \boldsymbol{\theta}}_{-1}) \mathbf{Y}/n$$

而 $\hat{ \boldsymbol{\theta}}_{-1}$ 是使如下似然最大化的任意 $\boldsymbol{\theta}_{-1} \in \boldsymbol{ \Theta}_{-1}$ 值：

$$L_{-1}(\boldsymbol{\theta}_{-1}; \mathbf{Y}) = - \frac{1}{2} \log |\mathbf{W}(\boldsymbol{\theta}_{-1})| - \frac{n}{2} \log(\mathbf{Y}^{\top} \mathbf{Q} (\boldsymbol{\theta}_{-1}) \mathbf{Y})$$

这里，

$$\mathbf{Q}(\boldsymbol{\theta}_{-1}) = \mathbf{W}^{-1}(\boldsymbol{\theta}_{-1}) - \mathbf{W}^{-1}(\boldsymbol{\theta}_{-1}) \mathbf{X}(\mathbf{X}^{\top} \mathbf{W}^{-1}(\boldsymbol{\theta}_{-1}) \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{W}^{-1}(\boldsymbol{\theta}_{-1})$$

这是另一种被应用于 尺度参数 $\theta_1$ （而不是均值参数 $\boldsymbol{\beta}$ ）的剖面化。最大化 $L_{-1}(\boldsymbol{\theta}_{-1}; \mathbf{Y})$ 而不是 $L(\boldsymbol{\theta}; \mathbf{Y})$ 的优势在于：最大化问题的维数从 $m$ 减少到 $m - 1$ 。这可以大大节省计算量。同样的方法样适用于 REML 对数似然函数。

（2）协方差矩阵存在特定结构的情况

当 $\boldsymbol{\Sigma} (\boldsymbol{\theta})$ （或 $W(\boldsymbol{\theta}_{-1})$ ）具有可加速行列式计算和含逆的二次型计算的某种结构时，也会减少基于似然推断的相关计算。

当 $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ 为 $n×n$ 时，所需计算量一般为 $\mathcal{O}(n^3)$ ；不过，对于某些模型和数据位置的空间配置，$\boldsymbol{\Sigma} (\boldsymbol{\theta})$ 中的结构可以显著减少计算量。

例如，当协方差函数的范围相对于观测空间域较小时，协方差矩阵将是 “稀疏的” （即有许多元素的值为零），而稀疏矩阵存在效果很好的计算方法（Barry 和 Pace，1997 ^[3]）。

最后，当协方差函数（或广义协方差函数）为各向同性或可分离时，如果数据位置形成了规则矩形网格，则可能会出现另一种计算可能大幅减少的情况（Zimmerman，1989a ^[38]，1989b ^[39]）。

4.7 近似似然和组合似然

（1）Vecchia 似然分解

除了上一节介绍的一些能够减少似然估计计算负担的精确方法之外，还有一种完全不同的计算策略，那就是用更容易最大化的函数来近似似然函数。

Vecchia (1988) 提出了地统计数据似然的近似值，它基于将观测向量 $\mathbf{Y}$ 划分为子向量 $\mathbf{Y}_1、\ldots 、\mathbf{Y}_b$ 。令 $\mathbf{Y}_{(j)} = (\mathbf{Y}^{\top}_1 ,\ldots , \mathbf{Y}_j^{\top})^{\top}$ ，令 $p(\mathbf{Y}; \boldsymbol{\beta}, \boldsymbol{\theta})$ 表示 $\mathbf{Y}$ 的联合密度，令 $p(\mathbf{Y}_j |\mathbf{Y}_{(j-1)}; \boldsymbol{\beta}, \boldsymbol{\theta})$ 表示给定 $\mathbf{Y}_{( j-1)}$ 时 $\mathbf{Y}_j$ 的条件密度，则（精确）似然可以改写为：

$$p(\mathbf{Y}; \boldsymbol{\beta}, \boldsymbol{\theta}) = p(\mathbf{Y}_1; \boldsymbol{\beta}, \boldsymbol{\theta}) \prod\limits^{b}_{j=2} p(\mathbf{Y}_j | \mathbf{Y}_{( j-1)}; \boldsymbol{\beta}, \boldsymbol{\theta}) \tag{4.3}$$

Vecchia (1988) 建议用 $\mathbf{\mathbf{Y}}_{(j-1)}$ 的子向量 $\mathbf{S}_{(j-1)}$ 来替换每个完整的条件向量 $\mathbf{\mathbf{Y}}_{(j-1)}$，进而得到 式（4.3） 的近似似然，此时计算逆矩阵和行列式的效率要高的多。

由于在此方法中没有划分 $\mathbf{Y}$ 的唯一方法，一个自然的问题是： $\mathbf{Y}_j$ 和 $\mathbf{S}_{(j-1)}$ 应该如何选择？

Vecchia 只考虑了长度为 $1$ 的向量 $\mathbf{Y}_j$ ，按数据位置的两个坐标轴中任一个的值排序，他建议选择 $\mathbf{S}_{(j-1)}$ 以包含最接近 $\mathbf{Y}_j$ 的 $q$ 个观测值，其中 $q$ 远小于 $n$ 。$q$ 的值越小，计算效率越高，但真实似然的近似越粗糙。对于从几个高斯随机场模拟的大小为 $n = 100$ 的数据集，Vecchia 表明取 $q = 10$ 就可以取得很好的近似；在实际工作中，$q$ 的选择与相对于观测点位之间距离的空间依赖性范围有关。

Stein、Chi 和 Welty (2004 ^[27]) 扩展了 Vecchia 的原始提议，允许 $\mathbf{Y}_j$ 的长度不均匀（和不一致），并将其应用于受限似然函数而不是普通似然函数。对于大小为 $n = 1000$ 的模拟高斯数据，他们发现取 $q = 16$ 足以使近似 REML 估计的（统计）效率达到精确 REML 估计的 $80\%$ 到 $95\%$ 范围内。然而，与 Vecchia 不同的是，他们发现在 $\mathbf{S}_{(j-1)}$ 中包含一些远距离观测值对估计有利。当空间依赖性相对于观测的空间域很强时尤其如此，这种情况在 Vecchia 的例子中没有体现过。

（2）组合似然法和伪似然法

一种类似但略有不同的近似似然估计策略是 伪似然 或 组合似然 （Besag，1975 年 ^[4]；Lindsay，1988 年 ^[19]; Varin 等 2011 年的 《组合似然法》）。

这种似然由度多个分量似然函数相乘得到，每个分量似然函数都是有效的条件似然或边缘似然，就好像这些分量似然对应于数据的独立子向量一样。通过这种方式，可以避免必须计算 $n × n$ 协方差矩阵的行列式和逆矩阵，并将这些计算替换为更小的矩阵。这种计算简单性的成本，与 Vecchia/Stein 等人一样，是由于忽略有关组分协方差结构的信息而导致的统计效率损失。

Curriero 和 Lele (1999 ^[7]) 考虑了一种特定的组合似然估计方案。对于本质平稳的高斯随机场，它们根据观测值之间所有成对差异的边缘密度形成复合对数似然函数。忽略常数项，此函数为

$$CL(\boldsymbol{\theta}; \mathbf{Y}) = - \frac{1}{2} \sum\limits^{n-1}_{i=1} \sum\limits_{j>i} \left \{ \log(\gamma (\mathbf{s}_i - \mathbf{s}_j ; \boldsymbol{\theta}) + \frac{(Y(\mathbf{s}_i ) - Y(\mathbf{s}_j ) )^2}{2 \gamma (\mathbf{s}_i - \mathbf{s}_j; \boldsymbol{\theta})} \right \}$$

$CL(\boldsymbol{\theta}; \mathbf{Y})$ 关于 $\boldsymbol{\theta}$ 的最大化产生复合最大似然估计量 (MLE)。请注意，通过仅考虑成对差异，根据 REML 的精神，未知常量均值参数已从该复合对数似然中消除。复合 MLE 在递增域渐近框架内的某些正则性条件下是一致且渐近高斯的，事实证明即使随机场不是高斯场，这种一致性也成立。此外，可以证明（Gotway 和 Schabenberger，2005 ^[11] ，第 171-172 页），这个特定的组合似然估计等效于通过在模型中应用（非线性）加权最小二乘法获得的估计：

$$[Y(\mathbf{s}_i ) - Y( \mathbf{s}_j)]^2 = 2\gamma (\mathbf{s}_i - \mathbf{s}_j ; \boldsymbol{\theta}) + \epsilon_{ij}, \quad \mathbb{E} (\epsilon_{ij}) = 0, \quad \operatorname{Var}(\epsilon_{ij}) = 8[\gamma (\mathbf{s}_i - \mathbf{s}_j ; \boldsymbol{\theta})]^2$$

（3）协方差矩阵锥化

近似似然的第三种策略是协方差逐渐减小（Kaufman、Schervish 和 Nychka，2008 ^[16]）。在这种方法中，对应于空间上较远的观测值对的协方差矩阵的元素被设置为零，因此可以使用上一节中提到的稀疏矩阵求逆和行列式计算算法。假设 $Y(·)$ 是具有线性均值函数和协方差函数 $C_0(h; \boldsymbol{\theta})$ 的各向同性高斯随机场，并设 $C_T(h; \gamma )$ 是一个 “锥化函数” ，也就是一个当 $h \geq \gamma$ 时等于零的各向同性连续相关函数（ $\gamma$ 已知）。然后令

$$C_1(h; \boldsymbol{\theta}, \gamma ) = C_0(h; \boldsymbol{\theta})C_T(h; \gamma ), \quad h \geq 0$$

请注意，$C_1(h; \boldsymbol{\theta}, \gamma)$ 是一个有效的（正定）协方差函数，对于较远的位置与 $C_0$ 不同，但保留与 $C_0$ 相同的方差。此外，Kaufman 等表明如果 $C_0$ 属于协方差函数的 Matérn 族，则可以找到一个锥化函数 $C_T$ ，使其产生的协方差函数 $C_1$ 在原点附近与 $C_0$ 具有相同表现，并且 $C_1$ 对应于与 $C_0$ 相应的零均值高斯的测量。现在，如果 $C_1$ (而不是 $C_0$) 事实上是 $Y(·)$ 的协方差函数，那么相应观测向量 $Y$ 的对数似然函数将为：

$$L_T(\boldsymbol{\theta,\beta}; \mathbf{Y}) =- \frac{1}{2} \log |\boldsymbol{\Sigma} (\boldsymbol{\theta}) \circ \mathbf{T}(\gamma )|- \frac{1}{2} (\mathbf{Y - X} \boldsymbol{\beta})^{\top} [\boldsymbol{\Sigma} (\boldsymbol{\theta}) \circ \mathbf{T}(\gamma )]^{-1}(\mathbf{Y - X}\boldsymbol{\beta})$$

其中 $\mathbf{T}(\gamma)$ 是第 $(i, j)$ 个元素 $C_T$ 的矩阵 $(\|\mathbf{s}_i - \mathbf{s}_j \|; \gamma )$ 和 “ $\circ$ ” 指的是两个矩阵的元素积或 Schur 积。Kaufman 等考虑通过最大化对数似然的近似值获得的 “单锥度估计器” ，以及另一个 “双锥度估计器” ，它基于无偏估计方程的理论，因此偏倚较小。他们通过模拟表明，两个估计器，尤其是 two-taper 估计器，表现得相当好，即使在逐渐变细的情况下也是如此。

（4）谱方法

最后一种策略是使用谱方法来近似似然。对于在规则网格点上观测到的平稳高斯随机场，Whittle (1954 ^[32]) 开发了以下近似似然：

$$L(\boldsymbol{\theta}; \mathbf{Y}) \doteq L_W(\boldsymbol{\theta}; \mathbf{Y}) =- \frac{n}{(2\pi)^2} \sum\limits_{\omega \in F} \{ \log f(\boldsymbol{\omega , \theta}) + I_n(\boldsymbol{\omega}; \mathbf{Y})[ f(\boldsymbol{ \omega} , \boldsymbol{\theta})]^{-1} \}$$

其中 $F$ 是傅里叶频率，$f$ 是随机场的谱密度，$I_n$ 是周期图。使用快速傅里叶变换，可以非常有效地计算 $L_W(\boldsymbol{\theta}; \mathbf{Y})$ ，只需 $\mathcal{O}(n \log_2 n)$ 次操作。 Fuentes (2007 ^[10]) 通过在网格单元上集成随机场 $Y(·)$ 扩展了这种方法，用于不规则间隔的数据。通过模拟，她发现这种方法在统计上与 Stein 等 (2004 ^[27]) 的方法一样有效。 除了在非常短的距离内捕获 $Y(·)$ 的行为外，它的计算效率要高得多。

4.8 非高斯数据的方法

到本章的这一点，我们假设观测值是高斯分布的。如何估计产生非高斯观测的连续空间过程的参数？对于连续但偏斜或有界支持（例如比例）的观测，一种似然是将数据转换为数据更接近高斯的尺度（参见 de Oliveira、Kedem 和 Short (1997 ^[8]) 的贝叶斯实现这个想法）。对于离散观测，例如计数或存在/不存在数据，可以基于一类称为空间广义线性混合模型 (GLMM) 的模型进行分析。我们以对 GLMM 及其参数的基于似然估计的简要讨论来结束本章。第 11 章介绍了其他一些非高斯数据模型和方法。

在经典的广义线性模型 (GLM) 中，对于具有均值 $\boldsymbol{\mu}$ 的独立观测向量，均值的某个函数 $g(\boldsymbol{\mu})$ ，称为链接函数，被假定为固定效应的线性组合，即 $g(\boldsymbol{\mu}) = \mathbf{X}\boldsymbol{\beta}$ 。这些影响和模型的任何分散参数可以通过最大似然或拟似然估计 (Wedderburn, 1974 ^[31])，这取决于是否假设数据分布（来自指数族）或仅假设其第一和第二时刻。 Liang 和 Zeger (1986 ^[18]) 以及 Zeger 和 Liang (1986 ^[34]) 将拟似然方法扩展到 GLM，并从纵向研究中获得序列相关的观测结果，Albert 和 McShane (1995 ^[1]) 以及 Gotway 和 Stroup (1997 ^[12]) 将这种扩展应用于空间数据。然而，空间 GLM 本身并不能提供一种建模和估计空间相关性的方法。 Diggle、Tawn 和 Moyeed (1998 ^[9]) 介绍的空间 GLMM 通过向模型添加空间相关的随机效应来提供这样的机会。在空间 GLMM 中，假定存在具有协方差函数 $C(\mathbf{h}; \boldsymbol{\theta})$ 的潜在零均值平稳高斯随机场 $\{ b(\mathbf{s}):\mathbf{s} \in D \}$ 。有条件地在 $b(·)$ 上，$Y(·)$ 被认为是一个独立的过程，其分布由条件均值 $E\{ Y(\mathbf{s})|b(\mathbf{s}) \}$ 指定，并且对于某些链接函数 $g$ ，

$$g[\mathbb{E}\{ Y (\mathbf{s})|b(\mathbf{s}) \}] = \mathbf{x}_j^T(\mathbf{s})\boldsymbol{\beta} + b(\mathbf{s})$$

空间 GLMM 参数的似然函数由

$$L(\boldsymbol{\beta}, \boldsymbol{\theta}; \mathbf{Y}) = \int_{\mathbb{R}^n} \{ \prod\limits^{n}_{i=1} f_i(Y(\mathbf{s}_i ) | \mathbf{b}, \boldsymbol{\beta}) \} f_{\mathbf{b}}(\mathbf{b}|\boldsymbol{\theta}) d \mathbf{b}$$

给出，其中 $\mathbf{b} = (b(\mathbf{s}_1,\ldots ,b(\mathbf{s}_n))^{\top}$ ，$f_i$ 是给定 $\mathbf{b}$ 的 $Y(\mathbf{s}_i )$ 的条件密度，$f_{\mathbf{b}}$ 是 $\mathbf{b}$ 的多元高斯密度。由于该积分的高维性，通过直接最大化 $L$ 来获得 MLE 是不可行的，因此提出了许多替代方案。Zhang (2002 ^[35]) 提出了一种用于最大化 $L$ 的期望最大化 (EM) 梯度算法。Varin、Host 和 Skare (2005 ^[28]) 和 Apanasovich、Ruppert、Lupton、Popovic (2007 ^[2]) 开发了复合（成对）似然估计方法。Diggle 等人 (1998 ^[9]) 以及 Christensen 和 Waagepetersen (2002 ^[5]) 使用贝叶斯 MCMC 方法获得参数估计。

参考文献

• [1] Albert, P.A. and McShane, L.M. (1995). A generalized estimating equations approach for spatially correlated data: Applications to the analysis of neuroimaging data. Biometrics, 51, 627–638.
• [2] Apanasovich, T.V., Ruppert, D., Lupton, J.R., Popovic, N., Turner, N.D., Chapkin, R.S., and Carroll, R.J. (2007). Aberrant crypt foci and semiparametric modeling of correlated binary data. Biometrics, DOI: 10.1111/j.1541-0420.2007.00892.x
• [3] Barry, R.P. and Pace, R.K. (1997). Kriging with large data sets using sparse matrix techniques. Communications in Statistics — Simulation and Computation, 26, 619–629.
• [4] Besag, J. (1975). Statistical analysis of non-lattice data. The Statistician, 24, 179–195.
• [5] Christensen, O.F. and Waagepetersen, R. (2002). Bayesian prediction of spatial count data using generalized mixed models. Biometrics, 58, 280–286.
• [6] Cressie, N. and Lahiri, S.N. (1996). Asymptotics for REML estimation of spatial covariance parameters. Journal of Statistical Planning and Inference, 50, 327–341.
• [7] Curriero, F.C. and Lele, S. (1999). A composite likelihood approach to semivariogram estimation. Journal of Agricultural, Biological, and Environmental Statistics, 4, 9–28.
• [8] de Oliveira, V., Kedem, B., and Short, D.A. (1997). Bayesian prediction of transformed Gaussian random fields. Journal of the American Statistical Association, 92, 1422–1433.
• [9] Diggle, P.J., Tawn, J.A., and Moyeed, R.A. (1998). Model-based geostatistics (with discussion). Applied Statistics, 47, 299–350.
• [10] Fuentes, M. (2007). Approximate likelihood for large irregularly spaced spatial data. Journal of the American Statistical Association, 102, 321–331.
• [11] Gotway, C.A. and Schabenberger, O. (2005). Statistical Methods for Spatial Data Analysis. Boca Raton, FL: Chapman & Hall.
• [12] Gotway, C.A. and Stroup, W.W. (1997). A generalized linear model approach to spatial data analysis and prediction. Journal of Agricultural, Biological, and Environmental Statistics, 2, 157–178.
• [13] Handcock, M.S. (1989). Inference for spatial Gaussian random fields when the objective is prediction. PhD dis., Department of Statistics, University of Chicago.
• [14] Harville, D.A. (1977). Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association, 72, 320–340.
• [15] Irvine, K.M., Gitelman, A.I., and Hoeting, J.A. (2007). Spatial designs and properties of spatial correlation: Effects on covariance estimation. Journal of Agricultural, Biological, and Environmental Statistics, 12, 450–469.
• [16] Kaufman, C.G., Schervish, M.J., and Nychka, D.W. (2008). Covariance tapering for likelihood-based estimation in large spatial datasets. Forthcoming. Journal of the American Statistical Association, 103, 1545–1555.
• [17] Kitanidis, P.K. (1983). Statistical estimation of polynomial generalized covariance functions and hydrologic applications. Water Resources Research, 19, 909–921.
• [18] Liang, K.-Y. and Zeger, S.L. (1986). Longitudinal data analysis using generalized linear models. Biometrika 73, 13–22.
• [19] Lindsay, B.G. (1988). Composite likelihood methods. Contemporary Mathematics, 80, 221–239.
• [20] Mardia, K.V. and Marshall, R.J. (1984). Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika, 71, 135–146.
• [21] Mardia, K.V. and Watkins, A.J. (1989). On multimodality of the likelihood in the spatial linear model. Biometrika, 76, 289–295.
• [22] Nelder, J.A. and Mead, R. (1965). A simplex method for function minimization. The Computer Journal, 7, 308–313.
• [23] Patterson, H.D. and Thompson, R. (1971). Recovery of inter-block information when block sizes are unequal. Biometrika, 58, 545–554.
• [24] Patterson, H.D. and Thompson, R. (1974). Maximum likelihood estimation of components of variance. In Proceedings of the 8th International Biometrics Conference, 1974, pp. 197–207. Alexandria, VA: The Biometric Society.
• [25] Rasmussen, C.E. and Williams, C.K.I. (2006). Gaussian Processes for Machine Learning. Cambridge, MA: MIT Press.
• [26] Spiegelhalter, D.J., Best, N.G., Carlin, B.P., and van der Linde, A. (2002). Bayesian measures of model complexity and fit (with discussion). Journal of the Royal Statistical Society Series B, 64, 583–639.
• [27] Stein, M.L., Chi, Z., and Welty, L.J. (2004). Approximating likelihoods for large spatial data sets. Journal of the Royal Statistical Society, Series B, 66, 275–296.
• [28] Varin, C., Host, G., and Skare, O. (2005). Pairwise likelihood inference in spatial generalized linear mixed models. Computational Statistics and Data Analysis, 49, 1173–1191.
• [29] Vecchia, A.V. (1988). Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society Series B, 50, 297–312.
• [30] Warnes, J.J. and Ripley, B.D. (1987). Problems with likelihood estimation of covariance functions of spatial Gaussian processes. Biometrika, 74, 640–642.
• [31] Wedderburn, R.W.M. (1974). Quasi-likelihood functions, generalized linear models and the GaussNewton method. Biometrika, 61, 439–447.
• [32] Whittle, P. (1954). On stationary processes in the plane. Biometrika, 41, 434–449.
• [33] Ying, Z. (1991). Asymptotic properties of a maximum likelihood estimator with data from a Gaussian process. Journal of Multivariate Analysis, 36, 280–296.
• [34] Zeger, S.L. and Liang, K.-Y. (1986). Longitudinal data analysis for discrete and continuous outcomes. Biometrics, 42, 121–130.
• [35] Zhang, H. (2002). On estimation and prediction for spatial generalized linear mixed models. Biometrics, 58, 129–136.
• [36] Zhang, H. and Zimmerman, D.L. (2005). Toward reconciling two asymptotic frameworks in spatial statistics. Biometrika, 92, 921–936.
• [37] Zimmerman, D.L. (1986). A random field approach to spatial experiments. PhD dis., Department of Statistics, Iowa State University.
• [38] Zimmerman, D.L. (1989a). Computationally exploitable structure of covariance matrices and generalized covariance matrices in spatial models. Journal of Statistical Computation and Simulation, 32, 1–15.
• [39] Zimmerman, D.L. (1989b). Computationally efficient restricted maximum likelihood estimation of generalized covariance functions. Mathematical Geology, 21, 655–672.
• [40] Zimmerman, D.L. and Zimmerman, M.B. (1991). A comparison of spatial semivariogram estimators and corresponding ordinary kriging predictors. Technometrics, 33, 77–91.