方差组分估计常用方法介绍

方差组分估计应该是遗传评估中最难的一部分,资料也不是很好找,尽量整理了一下。

ML

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

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

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

其对数为

对他求未知参数 的偏导数如下(证明略)

令上面两式为0,我们设 $ \frac{\partial \mathbf{V}}{\partial \sigma_{i}^{2}} = \mathbf{V}_{i}$ ,因此我们有

其中 ,因为下式成立

因为

故有 ML 的估计方程为

根据此式子,我们显然找不到闭式解,还需要采用迭代的方式求解。

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

REML

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

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

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

证明 的行 满足等式 ,转置得到 (注意这里出现了符号重叠问题,这里的 是第 行元素组成的列向量)。根据定理, 如果线性方程组 是相容的,那么所有可能的解为 ,这里 是某个特定的广义逆, 为所有可能的值组成的向量。因此,这里的方程组的解为 (书中有笔误,书里是 ,这里又有符号重叠问题,这里的 就是上面定理中的 ,为所有可能的 的向量) ,因此 $\mathbf{k}_{i}^{\prime}=\mathbf{c}^{\prime} \left(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\right) $ ,也就是说 的行 的行的所有可能的线性组合。

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

因为 $\mathbf{k}_{i}^{\prime}=\mathbf{c}^{\prime} \left(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\right) $ ,那么一定存在一个 的满秩矩阵 ,使得 矩阵可以表示为 ,其中 矩阵的每一行均是 的行的线性组合,并且由于 矩阵满秩, 因此选择 矩阵的标准是必须使得 矩阵的行之间线性无关。根据广义逆的性质,我们知道 的一个广义逆,因此我们可以这里的 可以选择为 ,因此此时 矩阵可以表示为

满足要求的 矩阵有无数个,但是这不影响我们使用。同时我们注意到校正固定效应的残差为 $\hat{\boldsymbol{\varepsilon}} = (\mathbf{I}-\mathbf{H}) \mathbf{y} $ ,因此

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

的分布见下面的定理。

定理:在混合线性模型中,假设 ,其中 ,那么

证明很简单,因为 满足正态,所以 也满足正态,同时根据 易得 的均值为 , 其方差为 $ \mathbf{K^{\prime} \mathbf{V}} \mathbf{K}$ ,得证。

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

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

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

证明:因为 的对数似然值为

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

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

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

证毕。

有一点很有意思,根据二次型的期望公式 ,上面式子的右手项的期望 (这里设期望公式里的 是 $\mathbf{K^{\prime}y} $ , ) 正好是左手项。

SEARLE (1979) 证明 (这篇文献没找到)

证明如下:

是对称幂等矩阵,对称证明略,幂等证明如下,$(\mathbf{K K}^{+})^{2} = \mathbf{K K}^{+} \mathbf{K K}^{+} = \mathbf{K K}^{+} $ 。我们已知 ,因此 并且 。因此

是一个对称幂等矩阵,幂等证明如下(其实直接根据对称幂等矩阵的性质即可证明,如果 对称幂等,那么 对称幂等 ):

因此

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

因为 是一个正定矩阵,因此总是存在 矩阵使得 。那么由于 ,那么我们知道 ,之前适应 的结果均适用于 和 $ \mathbf{V}^{-1/2} \mathbf{X}$ 。我们替换 中的MP逆矩阵,我们得到 。将 和 $ \mathbf{V}^{-1/2} \mathbf{X}$ 带入得到

因此 (注意下面的 是 MP 逆 )

得证。

其中 ,将该式带入,得到

注意最大似然法的估计公式为 $\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{V}{i}\right)
=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}
{i} \hat{\mathbf{P}} {\mathbf{y}} $ ,因此 REML 方法就是将左手项的 替换为

我们有

因此,我们有

因而,对于 ,我们有

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

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

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

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

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

EM 算法

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

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

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

因此,我们有

同时,我们有 (SEARLE, 1979;缺证明)

基于相同公式,另一种迭代公式为

这里 是第 个随机效应的水平数目; 是第 个随机效应的协方差矩阵; 是系数矩阵的广义逆矩阵中第 的随机效应对应的对角子矩阵; 是第 个随机效应的预测值。

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

而且根据上面两个公式,可以得到 都比大于0(我看不出来),因此 EM 算法一定满足参数空间的要求。

AI 算法

首先我们看 Newton-Raphson 算法,其迭代公式为:

其中 是需要最大化的原函数; 函数关于参数 的一阶偏导数向量; 函数关于参数 的二阶偏导数矩阵,即 Hessian 矩阵,举例如下

由于计算 矩阵的计算量较大,scoring 方法采用以下迭代公式,即采用 矩阵替换 矩阵。 矩阵是 矩阵的数学期望,即 矩阵则称为 Fisher 信息矩阵

所谓平均信息算法,就是将期望信息 () 和观察信息 () 平均起来,得到一个平均信息矩阵

此时相应的迭代公式为

在前面的推导中,我们已知 REML 方法的似然函数为

举个例子,假设模型中随机效应只有加性效应( )和残差 (),此时 Searle (1970) 证明 矩阵为

其中

矩阵为

其中左上角元素推导过程如下,其余元素同理(下面有一处应该有误,$ \boldsymbol{P X}= \boldsymbol{0} $ 不成立)。

此时平均信息矩阵为

此时在 AI 算法中,不涉及迹函数的运算,计算难度相对低一点。

对于迭代公式中用到的一阶偏导数 ,Johnson 和 Thompson (1995) 年证明存在下式

综上所述,可用如下步骤来实现 AI 算法,我们先给出一组初始的 值,由混合模型方程组求出相应的 ,以及相应的系数矩阵逆矩阵元素;再计算出 矩阵和一阶偏导数向量,再求出一组新的 值,如此迭代下去,知道两次迭代得到的 值得差异小于给定阈值时,迭代收敛。

DF 算法

非求导算法 (deritative-free, DF) 顾名思义即不计算导数对于不同的参数直接计算似然函数值,从中找到使得似然函数最大的一组参数。

参考文献

  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. Searle S R. An overview of variance component estimation[J]. Metrika, 1995, 42(1): 215-230.
  6. Notes on variance components estimation-A detailed account of maximum likelihood and kindred methodology
  7. 张沅, 张勤, 家畜育种. 畜禽育种中的线性模型[M]. 北京农业大学出版社, 1993.
  8. Searle S R, Casella G, McCulloch C E. Variance components[M]. John Wiley & Sons, 2009.
  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2022 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信