这是方差组分估计方法的第三篇博客，介绍非求导算法 (deritative free REML, DF-REML) 算法。非求导算法顾名思义即不计算导数，对于不同的参数直接计算似然函数值，从中找到使得似然函数最大的一组参数。不过我这里主要是介绍如何进一步推导似然函数。

推导 REML 似然函数

这里只是进一步推导REML的似然函数。

对于表型 ，其对数似然函数如下（ 为表型数目）

Verbyla (1990) 提出一个非奇异矩阵 ， 其中 和 分别为 和 的矩阵（ 为表型数目， 为 的列数，假设 满秩，如果不满秩可以提取其中一个满秩的子矩阵 ），并且满足

我们对表型 进行线性转换得到：

转换后的数据服从分布：

根据概率统计公式，我们有 ，因此其对数似然函数满足

首先对于 ，我们有

对于 ，我们有

根据正态分布的条件分布公式

我们得到

其协方差部分可以进一步推导：

因此，我们得到

其中 。

因此 的对数似然值可以写成

因为 ，因此左右两边的对数项也相等，得到：

化简得到 (因为 )

因此 REML （ ）的对数似然函数可以写为（前面两项是常数项，因为不包含参数）

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

因此，我们得到

对于 ，我们进行推导得到，其中 为 的维度。

因此 。

因此我们得到

代入 得到

左右乘以 -2 ，得到

忽略常数项，得到

ln|C|和 y’Py 的计算

这里有两种方法可以计算 和 ，第一种是是使用高斯消去法，我们建立混合线性模型方程组的增广矩阵如下：

利用高斯消去法消除其对角线下方的元素，将其变成一个上三角矩阵，此时最后一行的对角线元素就等于 。因为我们可以用等价的方式，将 吸收近 ，具体方式就是第二行子矩阵减去 乘以第一行（类似于元素级别的高斯消元，根据矩阵乘法的定义，这个行为是增广矩阵行之间的线性组合，和高斯消去法原理一样），得到

其中右下角元素进一步推导得到（根据之前 REML 公式的推导，我们得到 ）。

高斯消元后， 就是其前 n 个对角线元素的对数值之和，即

第二种方法是使用 Cholesky 分解，我们对系数矩阵进行分解，得到

通过三角方程组求解，我们可以很容易地求解，得到 。

通过推导，我们得到

此时我们可以通过 计算 。对于 ，我们有

如果我们对增广矩阵 进行 Cholesky 分解，则有

第一条证明略，第二条证明如下，根据舒尔补公式，我们有

根据Cholesky 分解，我们有

因此我们得到

DF-REML

非求导算法 (deritative-free, DF) 顾名思义即不计算导数对于不同的参数直接计算似然函数值，从中找到使得似然函数最大的一组参数。比如最简单的算法是将整个参数空间划分成 个网格，每个网格中取一个点，然后计算一次似然函数值，找出最大值。

目前DFREML 中推荐的方法主要是 Simplex 法，其基本思想如下，假设需要估计的方差组分数目为 ，那么我们首先选择 组方差组分初始值，每一组初始值就相当于 维空间中的一个点，这里也称为一个 Simplex ，然后我们通过比较这 组方差组分的似然函数值，淘汰其中最差的一个点，并通过某种方式产生一个新的点加入其中，如此不停迭代下去，不断向最大值运行，最终收敛于最大值。

当需要估计的方差组分数目较少时，DF-REML 算法较为适用，反之其收敛较慢。

附录1：行列式性质

对于矩阵 , 其中 可逆，则存在

第一个等式不用证明，我们证明第二个等式。

证明完毕。

参考文献

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

2. Verbyla A P. A conditional derivation of residual maximum likelihood[J]. Australian Journal of Statistics, 1990, 32(2): 227-230.

3. 分块矩阵的行列式