这是方差组分估计方法的第二篇博客，介绍ML和REML两种方法及一些很早之前的算法。

ML

在最大似然方法中，我们需要新增分布假设，一般我们都采用正态假设，因此模型假设为 (这里假设所有随机效应的协方差矩阵均为对角矩阵)

其中 是一个 并且秩为 的矩阵。 是一个 的正定矩阵。为了方便书写，我们记 ， ，因此协方差矩阵变成 (假设随机向量之间互不相关)

随机向量 的联合分布密度函数，即似然函数为

其对数为

对他求未知参数 的偏导数如下（证明略）

令上面两式为0，我们设 ，因此我们有

其中 ，因为下式成立

因为

故有 ML 的估计方程为

这个式子还可以换一种形式，对于左手项：

这里的 就是矩阵 所有元素的平方和。

同理右手项可以写为

因此，上面的式子可以写为

很显然这个式子不是 的线性方程组，左手项的 和右手项的 均包含 。根据此式子，我们显然找不到闭式解，还需要采用迭代的方式求解。我们可以设置 的一个初始值，然后得到 和 ，此时上面的方程组就变成了线性方程组，就可以求解得到一组新的 。

似然法的一个缺点式，当数据量比较小时，其估计值时有偏的。而且随着固定效应水平数目的增加和随机效应水平数目的减少，其偏差会愈加严重。

估计值的方差

由于 ML 估计值必须通过迭代求解得到，因此不可能求得估计值的准确的抽样方差，但当样本数目较大时，可以求得估计值的近似方差（称为大样本方差）。这个方差是借助参数的信息矩阵来求得的，参数 的信息矩阵用 表示，其定义为

其为二阶导的数学期望的负值。

Searle (1970) 证明，当样本数目足够大时，参数 的 ML 估计值 的方差近似为：

注意这式子本身并不需要 ML 估计值 。

ML 估计值的信息矩阵写为下式，其中 表示 ，其余式子类似。

因为我们有

因此，推导得到

同时我们有

因此，推导得到

然后再计算期望，我们已知 ，即 。我们再利用一个公式 , 其中 是一个非随机的矩阵。

证明过程如下，将 改成为

因此我们得到：

因此信息矩阵为

因此估计值的协方差矩阵为

The Hartley-Rao form

Hartley 和 Rao 在 1967 年提出了基于 矩阵的似然函数，其中 矩阵定义为

因此 矩阵和 矩阵具有相同的形式，不过 矩阵中的 变成了1，而 变成了 ，这意味着此时需要估计的参数为 ， 和 (这里的 就是上面的 ，随机效应数目)，此时我们将 单独拎了出来。

我们从之前的公式 出发进行推导，首先，当 时，此时 ， ，因此我们有

代入 ，得到

因此

然后我们继续看其它方差组分 () 的式子，对于公式 ，左右乘以 并对 的式子求和得到：

代入 ， 以及 ，得到

化简得到

等价于

根据上面的推导公式，方括号内的值为 0，因此我们得到 的计算公式

对于公式 ，化简得到

汇总一下，Hartley 和 Rao 方法中设计的方程组为

EM 算法

混合详细模型如下：

对于公式 ，我们在等式两边同乘以 ，并求总和，得到：

这里 为 矩阵维度，即表型数目。

我们再看右手项，在之前 ML 的推导中，我们知道 ，因此我们有

因此，我们有

其它方差组分公式的证明过程如下

首先我们有

这里我们进一步提出两个矩阵 和 矩阵，第一个矩阵如下，其中 ，即所有随机效应水平总数。

第二个矩阵是 矩阵的变体， 是将 矩阵中的所有 替换为 1，而其它 ( ) 替换为0 （注意这里的 矩阵是一个对角矩阵），因此我没有

举个例子，假设 ，那么我们有

根据 公式，我们可以推导出

对于其它方程组分我们有 的左手项，我们有

其中 是 矩阵中第 个子矩阵。

右手项可以写为

因为

因此右手项可以继续写为

因此我们将 改写为下式

因此我们得到了两个迭代公式

或者

这里 是第 个随机效应的水平数目； 是第 个随机效应的预测值。

同时使用公式

额外证明

1. 矩阵存在，因为 的行列式等于 。

   不会证明。

2. 

   证明：

   因为

   

   其逆矩阵为

   

   左右均乘以 得到

   

   这里 。

   其中 是一个正定矩阵，因为其二次型 ，因此其逆矩阵 也是正定矩阵，其对角线元素均大于 0。因为 是一个对角矩阵，根据上面的公式 ，因此 矩阵的对角线元素也均大于 0， 因此 矩阵正定，因此 。

3. 

   证明：根据上面的公式，我们有

   

   因此我们可以看出每一轮迭代产生的 均为正数，即均在参数空间中。（不知道如何证明 大于 0 ）

REML

限制性最大似然 (restricted (residual) maximum likelihood, REML) ，由 Patterson 和 Thompson 在 1971 年提出。我们强调 REML 方法的原因是在标准的线性模型中，样本方差 也是 REML 的估计值。同时，REML 是一种一般的方法，比如 REML 在不平衡数据中也可以使用。在某些平衡数据中，REML 可能有解析解 (闭式解)。REML 同时是最佳二次无偏估计值。

REML 的思想是对数据 进行最大似然估计，而不是 ，这里 是我们人为挑选的，以使得 的分布只包含方差组分，而不包含 。为了实现这一点，矩阵 需要满足 ，因此 。为了方便，我们同时要求 是一个满秩矩阵。我们同时要求 尽可能包含最多的关于方程组分的信息，因此 必须有最大的行数。

定理： 是一个满秩矩阵，满足 ，但是同时拥有最大的行数，因此它是一个 的矩阵 ( 为 的秩)。进一步地说， 矩阵的形式必须为 ，这里 是一个 的满秩矩阵。

证明： 的行 满足等式 ，转置得到 （注意这里出现了符号重叠问题，这里的 是第 行元素组成的列向量）。根据定理， 如果线性方程组 是相容的，那么所有可能的解为 ，这里 是某个特定的广义逆， 为所有可能的值组成的向量。因此，这里的方程组的解为 （书中有笔误，书里是 ，这里又有符号重叠问题，这里的 就是上面定理中的 ，为所有可能的 的向量） ，因此 ，也就是说 的行 是 的行的所有可能的线性组合。

我们知道 ，并且 ，因此 是一个幂等矩阵，因此 ，因此 的行的张成空间的维度为 ，因此最多有 个线性无关的向量 ，也就是说 矩阵的行数最大为 。

因为 ，那么一定存在一个 的满秩矩阵 ，使得 矩阵可以表示为 ，其中 矩阵的每一行均是 的行的线性组合，并且由于 矩阵满秩， 因此选择 矩阵的标准是必须使得 矩阵的行之间线性无关。根据广义逆的性质，我们知道 是 的一个广义逆，因此我们可以这里的 可以选择为 ，因此此时 矩阵可以表示为 。

满足要求的 矩阵有无数个，但是这不影响我们使用。同时我们注意到校正固定效应的残差为 ，因此

的每个元素都是所有残差的线性组合，这也是 residual maximum likelihood 这个名称的由来。

而 的分布见下面的定理。

定理：在混合线性模型中，假设 ，其中 ，那么

证明很简单，因为 满足正态，所以 也满足正态，同时根据 易得 的均值为 , 其方差为 ，得证。

因此 的分布只与 个方差组分有关。为了估计方差组分，REML 的下一步是针对这些方差组分最大化 的似然值。

对似然函数求偏导，使之为零，我们得到下面的定理。

定理：在上面的模型中， 个方差组分 的估计值满足下面的方程组，其中

证明：因为 ， 的对数似然值为

根据矩阵求导的性质，假设 是一个 的非奇异矩阵，其元素 是关于标量 的函数。我们定义 是一个 的矩阵，其元素为 。那么存在

和

利用这两条性质，我们求 对每一个 的偏导数，得到

将这个式子设为0，得到上面的结果。

证毕。

有一点很有意思，根据二次型的期望公式 ，上面式子的右手项的期望 (这里设期望公式里的 是 ， ) 正好是左手项。

SEARLE (1979) 证明

证明如下：

和 是对称幂等矩阵，对称证明略，幂等证明如下， 。我们已知 ，因此 并且 。因此

是一个对称幂等矩阵，幂等证明如下（其实直接根据对称幂等矩阵的性质即可证明，如果 对称幂等，那么 对称幂等 ）：

因此

根据幂等矩阵的性质，我们知道幂等矩阵的秩等于迹，因此 ，说明 （只有零矩阵的秩为 0）。因此，我们有

因为 是一个正定矩阵，因此总是存在 矩阵使得 。那么由于 ，那么我们知道 ，之前适应 和 的结果均适用于 和 。我们替换 中的MP逆矩阵，我们得到 。将 和 带入得到

因此 （注意下面的 是 MP 逆 ）

得证。

其中 ，将该式带入，得到

注意最大似然法的估计公式为 ，因此 REML 方法相比与 ML 方法就是将左手项的 替换为 。

我们有

因此，我们有

因而，对于 ，我们有

这就是 REML 的估计方程，它也必须通过迭代求解。

根据这个定理，我们可以得到 个方程，而我们总共有 个未知数。在某些情况下，这些方程组可以进一步简化从而得到闭式解。但是在大多数情况下，这些方程组是无法直接求解的。

我们一般会用一些迭代的方法来进行估计 (Rao 1997 pp. 104 – 105, McCulloch and Searle 2001, pp. 265 – 269)。我们注意到上面定理中的 个方程组成的方程组可以写为

这里 ，这里 是一个非奇异的 的矩阵，其中的元素 ，然后 是一个 的向量，其中的元素 。

我们注意到 这个式子貌似更加复杂，因为 和 中都含有未知数 。然而，这个方程组对于我们从一个初始值 逐步迭代来估计的方式而言很有用。我们在第 步可以用 来计算 和 ，那么 ，然后一直迭代直到收敛。

估计值的方差

首先我们证明

我们已知

求二阶导得到l

这里利用了两个性质： 和 ，证明过程参考 Linear models in statistics 。

再求期望得到信息矩阵， 这里利用二次型的期望公式

因此估计值的协方差矩阵为

这和 ML 的式子格式相同，只是将 替换成了 。

EM 算法

混合详细模型如下：

对于公式 ，我们在等式两边同乘以 ，并求总和，得到：

我们先看左手项 (下面的 是表型数目， 为 的秩)

我们再看右手项，在之前 ML 的推导中，我们知道 ，因此我们有

因此，我们有

对于某个随机因子的方差组分 ，其迭代公式推导如下。

首先我们需要证明一个引理，因为 ，如果我们将 中设置 ，将新的式子定义为 ，得到：

那么我们有下面的引理。

引理：

证明：因为我们有 ，代入 得到

因此我们有

注意这个公式和 有着相似的格式，只是其中的 替换成了

舒尔补公式说明如下，如果一个非奇异分块矩阵中，如果 可逆，则其逆矩阵可以表示为

其中矩阵 称为 的舒尔补，Marsaglia 和 Styan (1974a,b) 给出了两个重要公式，分别为（第二个式子可以理解为将 替换成了 ）。

我们再看公式 的左手项 ，之前在 ML 的 EM 算法推导中，我们得到

这里 是 的第 个子矩阵。

因为 和 的差别就是将 替换成了 ， 因此我们得到：

这里 是 中将 替换成了 的子矩阵，我们标记为

因为 REML 公式和 ML 公式的右手项不变，因此我们得到：

或者

通过上面这两个公式我们得到了 REML 的 EM 算法。在计算时，我们先给出一组初始的 值，由混合模型方程组求出相应的 和 ，以及相应的系数矩阵逆矩阵元素；再由上面的两个公式求出一组新的 值，如此迭代下去，知道两次迭代得到的 值得差异小于给定阈值时，迭代收敛。

而且根据上面两个公式，可以得到 和 都比大于0，因此 EM 算法一定满足参数空间的要求。

EM 算法

EM 算法的全名为 expectation-maximization ，因为它在计算条件期望值和最大化简化的似然函数值间不断跳转。EM算法只会计算估计值，而不会计算估计值的方差。

EM 算法用于当数据增加时会简化问题难度的最大似然估计，应用 EM 算法的核心在于决定什么是完整数据（观测数据+未观测到的数据）。实际数据（观测数据）在 EM 算法中一般称为不完整的数据。因此在方差组分估计中，我们将表型 视为不完整的数据，完整数据为表型 加上未观测的随机效应 () 。

如果我们知道随机效应 的值，我们可以用其平方和均值估计随机效应的方差，该式为基于完整数据的最大似然估计值。

但是，实际上我们并不知道随机效应 的值，但是 EM 算法给了我们一种基于 估计 的方法。首先我们采用一组参数的初始值，我们之后基于 计算 的条件期望（E 步, expectation step），然后利用其条件期望值使用最大似然法得到一组新的参数估计值（M 步， maximization step）。然后不断循环这两步直到收敛。

EM 算法的一个重要特征是，因为 EM 算法对完整数据执行最大似然估计，因此每一次迭代得到的估计值均在参数空间中。

联合分布

因为我们需要计算给定 条件下 的条件期望，因此我们需要 和 的联合分布。

我们有下面的式子

其中 , , 。

我们有

同时我们有

因此 和 的联合分布为 ，其中

因此我们有概率密度函数

其中

根据舒尔补，我们有公式

我们推导得到

同时根据舒尔补公式

我们得到

因此

因此，将上面的式子代入基于完全数据的对数似然函数（概率密度函数），得到

其中第二步推导是因为 。

因此根据完全数据的对数似然函数，对 求偏导，得到

使该式为 0 ，我们可以轻松得到 ML 估计值

对 求偏导，计算过程如下

因此

使之为 ，得到

即

这个公式等价于从 中剔除除了残差之外的随机效应，再应用普通最小二乘法得到的估计值。

为了完善 EM 算法的迭代步骤，我们还需要知道 和 在给定 时的条件期望。我们可以直接使用多变量正态分布的公式

进行推导得到

因此我们得到

根据二次型期望公式 ，我们可以推导得到

ML 的 EM 算法1

我们现在可以形成 ML 的 EM 算法的正式描述，其中上标 为第 m 次迭代得到的估计值

step 0: 给定一组初始值 和 ， 。

step 1 (E 步): 计算 和 这两个统计量的条件期望，将它们标记为 和 。

step 2 (M 步): 基于完整数据，采用最大似然法，得到

step 3: 如果达到了收敛标准，则将 ， ；不然 m 递增1，返回 step 1 。

ML 的 EM 算法2

如果我们使用 的 ML 估计值 ，则 的 ML 估计值为 。

Laird (1992) 建议在 EM 算法的迭代过程中不计算 和 ，而只是收敛时计算一次 ，推导过程如下，如果上面式子中的 均采用 的公式计算，即

而在计算 其中用到了 ，这正好是 ，其中 为：

因此此时在每个迭代过程中，我们不需要得到 ，我们只需要计算 。

这个算法其实严格来说不是 EM 算法，但是和 EM 算法差距很小，具体过程如下：

step 0: 给定一组初始 ， 。

step 1 (E 步): 计算 的条件期望，标记为 。

step 2 (M 步): 基于完整数据，采用最大似然法，得到

step 3: 如果达到了收敛标准，则将 ， ；不然 m 递增1，返回 step 1 。

EM 算法和 ML 方程组的等价性

我们现在考虑EM 算法和 ML 方程组的等价性，首先当 EM 算法收敛时，我们有 ， ，利用第一个等式，根据 ML 的 EM 算法1， 我们得到

根据广义逆性质我们有 ，因此化简得到

等式左右两侧乘以 , 根据广义逆性质 , 我们得到

这就是 ML 方程组中的一个式子。

类似地，对于 ，我们有

只要 ，我们可以简单得到

这是 ML 方程组中的另外一个式子。

REML 的 EM 算法

我们只需要将 ML 的 EM 算法中的矩阵和向量进行以下替换：

然后我们就可以从 ML 的 EM 算法2 推导得到 REML 的 EM 算法。

其中会应用到下面的式子

我们先推导一下ML 的 EM 算法2式子中的 矩阵，将其中的 和 矩阵进行替换，得到：

因此ML 的 EM 算法2式子中的 ，因此其在 REML 的算法中仍为 不变。

因此，我们得到 REML 的 EM 算法步骤如下：

step 0: 给定一组初始 ， 。

step 1 (E 步): 计算 的条件期望，标记为 。

step 2 (M 步): 基于完整数据，采用最大似然法，得到

step 3: 如果达到了收敛标准，则将 ；不然 m 递增1，返回 step 1 。

我们可以看到 REML 的 EM 算法与 ML 的 EM 算法只有一点不同，就是将公式中的 替换成了 。

参考文献

1. Rencher A C, Schaalje G B. Linear models in statistics[M]. John Wiley & Sons, 2008.
2. Patterson H D, Thompson R. Recovery of inter-block information when block sizes are unequal[J]. Biometrika, 1971, 58(3): 545-554.
3. Hofer A. Variance component estimation in animal breeding: a review[J]. Journal of Animal Breeding and Genetics, 1998, 115(1‐6): 247-265.
4. Misztal I. Reliable computing in estimation of variance components[J]. Journal of animal breeding and genetics, 2008, 125(6): 363-370.
5. 张沅, 张勤, 家畜育种. 畜禽育种中的线性模型[M]. 北京农业大学出版社, 1993.
6. Searle S R. Large sample variances of maximum likelihood estimators of variance components using unbalanced data[J]. Biometrics, 1970: 505-524.
7. Searle S R. Notes on variance components estimation-A detailed account of maximum likelihood and kindred methodology[J]. Technical Report BU-673-M, 1979.
8. Searle S R. An overview of variance component estimation[J]. Metrika, 1995, 42(1): 215-230.
9. Searle S R, Casella G, McCulloch C E. Variance components[M]. John Wiley & Sons, 2009.
10. Matsaglia G, PH Styan G. Equalities and inequalities for ranks of matrices[J]. Linear and multilinear Algebra, 1974, 2(3): 269-292.
11. Marsaglia G, Styan G P H. Rank conditions for generalized inverses of partitioned matrices[J]. Sankhyā: The Indian Journal of Statistics, Series A, 1974: 437-442.
12. Bard, Y. (1974). Nonlinear Parameter Estimation. Academic Press, New York.