方差组分估计方法三之DF-REML

这是方差组分估计方法的第三篇博客,介绍非求导算法 (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. 分块矩阵的行列式

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2023 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信