最大似然估计（ML）和REML

A.F. Zuur et al., Mixed Effects Models and Extensions in Ecology with R, p116-119, chr 5.6

当利用混合效应建模时，你会遇到诸如REML和ML这样的词汇。不像线性回归模型，就算你不知道背后的数学理论也可以照样使用。但是在混合效应建模中，你必须懂得一些相关的数学知识。所以REML是什么意思？它能干什么？
第一个问题比较简单，REML限制性最大似然估计英文首字母的缩写。但是对于第二个问题，大多数教材在这点上就变得相当专业，或者解释得比较粗略，只是提到它是矫正自由度的一个神秘方法【ML没有考虑到估计固定效应带来自由度的损失，造成参数的低估】^[1]。而我们则选择尝试更为详细地解释它，所以需要利用矩阵代数的知识。但是为了理解REML，首先需要理解最大似然估计的原理，我们从这儿开始。如果你不熟悉矩阵代数，或者如果这一部分对数学水平的要求太高，我们仍建议你跳过这一部分。
我们首先回顾用于线性回归的最大似然，然后给出REML是如何用于矫正方差估计量的。
假设有一个线性回归模型
$Y_{i}=\alpha+\beta \times X_{i}+\varepsilon_{i}$ ，其中 $\varepsilon_{i}\sim N(0,\sigma^2)$ 。模型中有3个未知参数 $\alpha$ 、 $\beta$ 和 $\sigma$ 。为了简便，令 $\boldsymbol{\theta}=(\alpha, \beta, \sigma)$ 。普通最小二乘是估计 $\boldsymbol{\theta}$ 的一个方法，它给出 $\boldsymbol{\theta}$ 中每一个元素的表达式。利用线性回归获取方差估计的表达式是：
$\hat{\sigma}^{2}=\frac{1}{n-2} \sum_{i=1}^{n}\left(Y_{i}-\hat{\alpha}-\hat{\beta} \times X_{i}\right)^{2} \quad\quad(5.14)$
我们给参数加了一个帽子^，表示它是估计值， $n$ 是观测值的个数。可以证明 $\hat{\sigma}$ 是 $\sigma$ 的无偏估计，意味着 $E[\hat{\sigma}]=\sigma$ 。现在让我们看一看最大似然估计方法。假设 $Y_i$ 服从正态分布，其密度函数为：
$f_{i}\left(Y_{i}, X_{i}, \alpha, \beta, \sigma\right)=\frac{1}{\sigma \sqrt{2 \pi}} \mathrm{e}^{\frac{\left(V_{i}-\alpha-\beta \times X_{i}\right)^{2}}{2 \sigma^{2}}}\quad\quad(5.15)$

因为我们也假设了 $Y_i$ 是独立的，可以将 $Y_{1}, Y_{2}, \dots, Y_{n}$ 的联合密度函数写成单个密度曲线 $f_{1}, f_{2}, \dots, f_{n}$ 乘积的形式。这个乘积就叫做似然函数L。它是数据和 $\boldsymbol{\theta}$ 的一个函数。问题是如何选择 $\boldsymbol{\theta}$ 使L最大。为了简化，对L取自然对数，得到如下的log似然方程式：
$\ln L\left(Y_{i}, X_{i}, \alpha, \beta, \sigma\right)=-\frac{n}{2} \ln 2 \pi-\frac{n}{2} \ln \sigma^{2}-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{n}\left(Y_{i}-\alpha-\beta \times X_{i}\right)^{2}\quad\quad(5.16)$
我们需要最大化这个式子，问题就变成了对每一个参数偏导，令偏导数=0并求解。因为我们很容易计算，这些偏导数=0的式子称之为封闭解。对于广义线性混合模型我们将会看到开放解，意思是参数没有直接的解。
这里没有给出 $\alpha$ 和 $\beta$ 估计量的式子，但对于方差我们得到：
$\hat{\sigma}^{2}=\frac{1}{n} \sum_{i=1}^{n}\left(Y_{i}-\hat{\alpha}-\hat{\beta} \times X_{i}\right)^{2}\quad\quad(5.17)$
注意这个式子与我们利用普通最小二乘得到的式子（5.14）非常相似。实际上，受到因子 $(n-2) / n$ 的影响，利用最大似然得到的方差估计量是有偏的（回归分析为什么误差方差中自由度是n-2？）。如果线性回归模型含有p个解释变量，那么偏度是 $(n-p) / n$ 。最大似然是有偏的的原因是它忽略了截距和斜率也被估计的事实。所以我们需要更好的ML估计量，而这正是REML所做的事情。

REML的工作如下：有线性回归模型 $Y_{i}=\alpha+\beta \times X_{i}+\varepsilon_{i}$ 可以写成 $Y_{i}=\mathbf{X}_{i} \times \boldsymbol{\beta}+\varepsilon_{i}$ 。这是简单的矩阵形式， $\mathbf{X}_{i}=\left(1 X_{i}\right)$ ， $\boldsymbol{\beta}$ 的第一个元素是截距，第二个元素是原始的 $\beta$ 。正态性假设意味着
$Y_{i} \sim N\left(\mathbf{X}_{i} \times \beta, \sigma^{2}\right)\quad\quad(5.18)$
用ML估计量的问题是我们不得不估计式子5.18中 $\boldsymbol{\beta}$ 中的截距和斜率。显然，如果没有 $\boldsymbol{\beta}$ ，就能解决问题。为了消除 $\boldsymbol{\beta}$ ，可以找到一个维度 $n \times (n – 2)$ 的特殊矩阵 $\mathbf{A}$ ，特殊指的是“与 $\mathbf{A}^\prime$ 正交”，然后用这个矩阵乘以 $\mathbf{Y}=\left(Y_{1}, \dots,Y_{n} \right)^{\prime}$ 之后再用ML估计。正交指的是如果 $\mathbf{A}$ 与 $\mathbf{X}$ 相乘，结果是0。因此，我们得到 $\mathbf{A}^{\prime} \times \mathbf{Y}=\mathbf{A}^{\prime} \times \mathbf{X} \times \boldsymbol{\beta}+\mathbf{A}^{\prime} \times \boldsymbol{\varepsilon}=\mathbf{0}+\mathbf{A}^{\prime} \times \boldsymbol{\varepsilon}=\mathbf{A}^{\prime} \times \boldsymbol{\varepsilon}$ 。现在 $\mathbf{A}^{\prime} \times \mathbf{Y}$ 的分布是
$\mathbf{A}^{\prime} \times \mathbf{Y} \sim N\left(\mathbf{0}, \sigma^{2} \times \mathbf{A}^{\prime} \times \mathbf{A}\right)$
而不再依赖 $\boldsymbol{\beta}$ 。那么对 $\mathbf{A}^{\prime} \times \mathbf{Y}$ 进行似然估计就会得到 $\sigma^{2}$ 的无偏估计量（5.14）。现在我们讨论REML如何应用到混合线性模型。我们的起点是边际模型
$\mathbf{Y}_{i} \sim N\left(\mathbf{X}_{i} \times \boldsymbol{\beta}, \mathbf{V}_{i}\right) \quad \mathbf{V}_{i}=\mathbf{Z}_{i} \times \mathbf{D} \times \mathbf{Z}_{i}^{\prime}+\mathbf{\Sigma}_{i}$
故事又重新开始，如之前，我们可以写一个略微不同的log似然式子。未知参数是 $\boldsymbol{\beta}$ 和 $\mathbf{D}$ 及 $\mathbf{\Sigma}_{i}$ 中的元素，依然用 $\boldsymbol{\theta}$ 表示。似然函数：
$\ln L\left(\mathbf{Y}_{i}, \mathbf{X}_{i}, \boldsymbol{\theta}\right)=-\frac{n}{2} \ln 2 \pi-\frac{1}{2} \sum_{i=1}^{n} \ln \left|\mathbf{V}_{i}\right|-\frac{1}{2} \sum_{i=1}^{n}\left(\mathbf{Y}_{i}-\mathbf{X}_{i} \times \boldsymbol{\beta}\right)^{\prime} \times \mathbf{V}_{i}^{-1} \times\left(\mathbf{Y}_{i}-\mathbf{X}_{i} \times \boldsymbol{\beta}\right)$
$|\mathbf{V}_{i}|$ 是 $\mathbf{V}_{i}$ 的行列式。对 $\boldsymbol{\beta}$ 求偏导并=0解方程。如之前讨论的例子，得到的参数是有偏的，因此我们需要REML。
总之，REML，就是用一个特殊的矩阵乘以Y，这样X×β就消去。然后用ML估计得到的参数估计子就是无偏的，并且与特定的矩阵相乘无关。因此， $\boldsymbol{\beta}$ 的REML估计子与ML的估计子不同。如果相对于观测值的个数，固定协变量的个数很少，就没有太大的不同，相反有许多的固定协变量，情况就大不相同。

Lindstrom MJ, Bates DM (1988) Newton—Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data. J Am Stat Assoc 83:1014–1022. ↩

最大似然估计（ML）和REML

推荐阅读更多精彩内容