这是方差组分估计方法的第四篇博客，介绍 EM-REML 算法。

混合线性模型的一般式子

我们将混合线性模型写为

其中 和 的联合分布为

其中 ， ， 。这里我们假定这三个矩阵均为正定矩阵。

这里的 和 是与随机效应 和 相关的方差组分的向量，而标量参数 和 是与残差和全部方差组分的缩放系数。因此我们需要估计的方差参数为 和 。

因为 的分布为

其中 是一个正定矩阵。

这里的标量参数 和 均是用于缩放的参数，我们发现如果我们同时包含这两个参数，那么这两个参数的值便均无法确定，因为此时 ， 。实际上我们至少会将其中一个参数设定为 1 。

按照育种模型的一般思路，这里我将 和 均设为 1， 此时 ，我们需要估计的参数为 。

假设我们有 个随机效应，即 ，因此此时随机效应的协方差矩阵可以写为

此时的方差组分 也可以进行相应的分块，例如 ，因此我们有 。

但是对于残差，我们通常还是认为所有数据来自与一个数据集，即残差是均质的，其协方差矩阵不用分块。

推导 C 逆矩阵

假设 矩阵满秩，则混合线性模型的系数矩阵 可逆。

我们将原始的系数矩阵分块为

根据舒尔补公式，系数矩阵逆矩阵可以写为

在之前的推导中，我们已经得到

根据附录1，我们有 ，因此其它分块子矩阵推导如下

因此我们得到

推导 P 矩阵

之前我们得到了两个 P 矩阵的式子，如下（ 为满足 的一个矩阵）

这里我们推导第三个式子

其中 ， 为系数矩阵。

证明过程

根据下面的附录1，我们有

我们可以将其写为

其中 。

我们也可以用 来表示

我们同样可以用 来表示 矩阵的元素

因此系数矩阵逆矩阵可以表示为

因此，我们有

因此

证明完毕

EM-REML

在之前的DF-REML 推导中，我们得到 REML 的一个似然函数的式子如下：

对于某个参数 进行求导，如下：

我们定义 ，并且使用下面的求导公式

因此我们得到

因此我们得到 REML score 式子：

下面我们要考虑 REML score 对于 和 的具体式子。

对于一个与 相关的方差组分 ， 推导如下

其中 。

如果模型中只有一个随机组分，则有：

对于残差的方差组分 则有

因此我们得到对于 的 REML score 式子为：

我们定义 是一个 的矩阵（ 为所有固定因子和所有随机因子的水平数目， 是第 个随机因子的水平数目 ），并且满足 ，因此 中只有对应 在 中的位置存在一个单位矩阵，其它元素全是 0 （因为矩阵向量乘积可以理解为矩阵的列的线性组合，因此 每一列只有一个对应 的列号的元素为1，其它为0 ）。

我们将系数矩阵 写为

其中

我们对上面式子进一步推导，得到

其中 可以写为：

这里 就是 中对应着随机效应的部分，由于我们可以有多个随机效应 ，因此 可以进一步分块为（其中 为随机效应数目）：

其中 就是对应着第 个随机效应 的部分。

同样的道理，我们可以将 矩阵分块为：

因此我们得到

因为 中只有对应 在 中的位置存在一个单位矩阵，其它元素全是 0 ，因此我们得到

因此， 可以进一步推导得到

对第二个式子推导前，我们先用 得到随机效应 的式子，根据BLUP方法，我们有

同时我们有

因此，我们得到

我们同时得到

因此，对右手项推导得到

因此，我们得到对于 的 REML score 式子可以写为：

如果模型中只有一个随机效应，则有

对于残差的方差组分 ，我们将 代入 得到

我们对其进行推导得到：

对于第二个式子，首先我们推导 ，如下

因此，我们得到

因此我们得到关于 的 REML score 式子如下

单性状模型的迭代公式

对于单性状模型（假设所有随机因子之间均不存在相关），我们将 设为 0， 推导一下，我们得到：

整理一下，我们得到（其中 为 的水平数目，或者说是 的列数）

一个收敛更快的迭代公式为

联合使用之前推导得到的残差计算公式(其中 是表型数目， 为 的秩，只适合单性状模型)，我们便得到了完整的 EM-REML 估计公式

这里我们还有残差的另外两种推导公式，如下

将 设为0，我们得到

因此可以得到第一个式子，我们对其中 推导如下

对于 ，我们对其推导得到

因此，我们得到

因此我们得到

代入上式，我们得到

因此可以得到第二个式子。

多性状模型的迭代公式

随机效应

对于多性状模型，假设只有一个随机因子（其实多个因子同理，只是这样不用写 矩阵的下标了，符号更加清晰），按照性状排序， 则有

其中 就是相应的方差组分的矩阵形式，如下所示，即反过来 ，其中 为性状数目。

对于随机因子 ，假设对于每个性状的随机因子水平数目均为 ，则 也可以表示为一个 的矩阵 ，其中每一列代表一个性状的所有水平，如下

反过来则有 。

对于 ，我们有公式

因为 相应 位置的元素 就是 ，所有易得 ，其中 。

因此，我们推导得到。

最后一个式子中，我们设 。

代入 得到

我们需要对这个式子进一步简化，我们知道 是一个对称矩阵，因此 也是一个对称矩阵， 我们将 中的元素表示为 ( ) 。

首先我们考虑 这一项，当 时，矩阵 只有 位置元素为1，其它位置元素均为 0 。因此 只有第 列元素为 的第 列元素，其它元素为 0 。因此，我们得到

当 时, 矩阵 只有 和 位置为 1，其它位置元素为 0 。因此 其第 列为 的第 列， 其第 列为 的第 列，其它位置元素为 0 。因此我们得到

即我们得到

我们现在考虑第三项，如下

其中我们设 。

当 时，矩阵 只有第 列元素为 的第 列元素 ，其它位置元素为 0 。因此，我们得到

当 时, 矩阵 其第 列为 的第 列， 其第 列为 的第 列，其它位置元素为 0 。因此，我们得到 (因为 对称)

综上，我们有

我们现在考虑第二项，我们将 按照性状分块为 个 的子矩阵，如下

处于简化的目的，我们这里考虑 的情况（这个不影响推导出来的结果），根据矩阵乘法，当 时，我们有

因此，我们得到

利用迹函数的性质，我们得到

当 时, 我们有

因此

因此我们得到

我们对需要求和的第一项交换 和 ，并且我们知道 , ，因此我们继续推导得到

综上，我们将 改为一般的 ，我们得到

现在我们定义 是一个 的矩阵，其元素为 ，根据矩阵乘法性质（假设 均为 的矩阵，则 ），我们发现

我们定义 是 的元素，即 。

因此，我们得到

将上面的式子代入 ，进一步简化得到

将之设为0， 对 求解得到 EM 推导公式

使用矩阵格式，我们得到

左右均乘以 ，得到

因此，我们得到 EM 推导公式

因为 ， 可以进一步分块为：

因此 。因此上面的 EM 推导公式可以写为

这个式子和单性状模型的推导公式正好具有相同的格式。

残差

按照性状排序， 此时我们有下式成立，其中 为表型行数， 为需要估计的协方差组分的矩阵。

形式如下（这里我们使用 ，与上面的随机效应格式相同 ）。

对于随机因子 ，同上， 也可以表示为一个 的矩阵 ，其中每一列代表一个性状的所有残差，如下

反过来则有 。

对于 ，我们有公式

易得 ，其中 。

因此，我们推导得到

最后一个式子中，我们设 。

代入 得到

我们需要对这个式子进一步简化，我们知道 是一个对称矩阵，因此 也是一个对称矩阵， 我们将 中的元素表示为 ( ) 。

首先我们考虑第一项 ，同上，我们可以得到

我们现在考虑第三项，如下

其中我们设 。

当 时，矩阵 只有第 列元素为 的第 列元素 ，其它位置元素为 0 。因此，我们得到

当 时, 矩阵 其第 列为 的第 列， 其第 列为 的第 列，其它位置元素为 0 。因此，我们得到 (因为 对称)

综上，我们有

对于第二项，首先我们对 进行推导如下

因此，第二项推导得到

对于第二项的第一项进行推导（其余项同理可得），假设 ，则 可以分解为

类似与随机效应，我们可以得到（推导过程略）

因此，我们得到（同上，我们有 ）

现在我们定义 是一个 的矩阵，其元素为 ，根据矩阵乘法性质（假设 均为 的矩阵，则 ），我们发现

我们定义 是 的元素，即 。

因此，我们得到

将上面的式子代入 ，进一步简化得到

将之设为0， 对 求解得到 EM 推导公式

使用矩阵格式，我们得到

左右均乘以 ，得到

因此，我们得到

因为 ， 可以进一步分块为：

因此 。因此上面的 EM 推导公式可以写为

如果我们设 , , ，则我们有

则 的 EM 推导公式可以写为

附录1: V逆相关推导

根据舒尔补公式 ， 我们可以得到

两边左乘 得到

因此我们得到

转置一下，得到

附录2：P 矩阵对方差组分的导数

我们有

证明如下：

参考文献

1. Knight E. Improved iterative schemes for REML estimation of variance parameters in linear mixed models[D]. , 2008.