[论文解读] Essential formulae for restricted maximum likelihood and its derivatives associated with the linear mixed models

Shengxin Zhu, Andrew J. Wathen|arXiv (Cornell University)|May 11, 2018

Spectroscopy and Chemometric Analyses参考文献 16被引用 9

一句话总结

本文提供了线性混合模型中受限最大似然（REML）及其导数的自包含数学推导，重点在于高效计算观测信息矩阵与费雪信息矩阵。文章提出了关键REML公式的新证明，包括误差对比变换和平均信息矩阵，从而实现高通量生物数据的可扩展优化。

ABSTRACT

The restricted maximum likelihood method enhances popularity of maximum likelihood methods for variance component analysis on large scale unbalanced data. As the high throughput biological data sets and the emerged science on uncertainty quantification, such a method receives increasing attention. Estimating the unknown variance parameters with restricted maximum likelihood method usually requires an nonlinear iterative method. Therefore proper formulae for the log-likelihood function and its derivatives play an essential role in practical algorithm design. It is our aim to provide a mathematical introduction to this method, and supply a self-contained derivation on some available formulae used in practical algorithms. Some new proof are supplied.

研究动机与目标

• 提供实际算法中使用的REML公式的数学严格推导，弥补先前统计文献留下的空白。
• 通过提供新的自包含证明，纠正并更新过时的公式，例如Zhou等人（2012）的公式。
• 通过阐明观测信息矩阵与费雪信息矩阵的结构，实现高维方差成分估计中基于导数的高效优化。
• 通过平均信息矩阵建立可扩展REML估计的基础，尤其适用于大规模非平衡数据集。
• 通过提供对数似然导数的计算可处理表达式，支持基因组学和不确定性量化中的算法设计。

提出的方法

• 使用正交投影矩阵和帽子矩阵的谱分解 $ P_X = X(X^T X)^{-1}X^T $ 推导误差对比变换，证明存在 $ L = [L_1, L_2] $ 使得 $ L_1^T X = I_p $，$ L_2^T X = 0 $。
• 通过 $ L_2^T y $ 的边缘分布（与固定效应无关）构造受限对数似然 $ ilde{ heta}^{ ext{REML}} = ext{argmax}_ heta ilde{ heta} $。
• 推导观测信息矩阵 $ oldsymbol{/mathcal{I}}_O( heta_i, heta_j) $ 和费雪信息矩阵 $ oldsymbol{/mathcal{I}}_F( heta_i, heta_j) $，并展示其复杂的基于迹的表达形式。
• 引入平均信息矩阵 $ oldsymbol{/mathcal{I}}_A( heta_i, heta_j) $，通过二次型简化计算，实现高效的牛顿-拉夫森优化。
• 证明平均信息矩阵近似于观测信息矩阵与费雪信息矩阵的精确平均，其残差项可忽略且为零矩阵。
• 提出通过具有多个右端项的稀疏线性系统求解平均信息矩阵，从而实现在大规模生物数据上的可扩展计算。

实验结果

研究问题

• RQ1 如何严格构造误差对比变换，以确保 $ L_2^T y $ 的边缘分布与固定效应无关？
• RQ2 在REML估计中，观测信息矩阵与费雪信息矩阵的正确且计算高效的表达式是什么？
• RQ3 如何推导并使用平均信息矩阵，以在大规模问题中替代计算成本高昂的观测信息矩阵？
• RQ4 为何一些现有算法仍使用过时的公式？如何通过新推导加以纠正？
• RQ5 在高维REML估计中，使用牛顿-拉夫森等基于导数的方法相比无导数方法具有何种计算优势？

主要发现

• 提供了误差对比变换 $ L = [L_1, L_2] $ 存在性与构造的严格证明，确保 $ L_1^T X = I_p $ 且 $ L_2^T X = 0 $，这对REML至关重要。
• 受限对数似然推导为 $ ilde{ heta} = -rac{1}{2} igracevert (n - u) ext{log}(2 ilde{ heta}) + ext{log}|L_2^T V( heta) L_2| + y^T L_2 (L_2^T V( heta) L_2)^{-1} L_2^T y igracevert $，可消除固定效应估计的偏差。
• 观测信息矩阵 $ oldsymbol{/mathcal{I}}_O( heta_i, heta_j) $ 涉及四个矩阵乘积的迹，对大规模数据集而言计算成本过高。
• 费雪信息矩阵 $ oldsymbol{/mathcal{I}}_F( heta_i, heta_j) $ 虽更简单，但仍包含需高成本计算的迹项。
• 平均信息矩阵 $ oldsymbol{/mathcal{I}}_A( heta_i, heta_j) $ 被证明是观测信息矩阵与费雪信息矩阵的精确平均，其残差项可忽略且为零矩阵。
• 平均信息矩阵可通过具有多个右端项的稀疏线性系统高效计算，从而实现在高通量基因组学中的可扩展REML优化。