知经 KNOWECON
登录 注册

知经百科 / X

限制最大似然

限制最大似然 (REML)

限制最大似然（Restricted Maximum Likelihood，简称 REML），又称残差最大似然，由统计学家 Patterson 和 Thompson 于 1971 年提出，是方差分量估计和混合效应模型中最为核心的参数估计方法之一。与标准的极大似然估计（ML）不同，REML 通过线性变换将数据投影到固定效应的误差对比空间（error contrast space），在此变换后的空间中对方差参数进行似然推断，从而消除固定效应估计所带来的自由度损失对方差分量估计的偏差。该方法在生物统计、数量遗传学、计量经济学的面板数据分析和分层线性模型中广泛应用，是现代随机效应模型估计的基准方法。

动机：标准ML的方差分量偏差

理解 REML 的动机，必须从标准极大似然估计在方差分量估计中的系统性偏差入手。考虑线性混合模型：

y=Xβ+Zu+ϵ\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{u} + \boldsymbol{\epsilon}

其中 y\mathbf{y} 为 n×1n \times 1 响应向量，X\mathbf{X} 为固定效应设计矩阵，β\boldsymbol{\beta} 为 p×1p \times 1 固定效应向量，Z\mathbf{Z} 为随机效应设计矩阵，u∼N(0,G)\mathbf{u} \sim N(0, \mathbf{G}) 为随机效应，ϵ∼N(0,R)\boldsymbol{\epsilon} \sim N(0, \mathbf{R}) 为误差项。令 θ\boldsymbol{\theta} 表示 G\mathbf{G} 和 R\mathbf{R} 中包含的所有未知方差参数。

标准的 ML 估计通过同时最大化关于 β\boldsymbol{\beta} 和 θ\boldsymbol{\theta} 的联合似然函数获得。然而，ML 估计量 θ^ML\hat{\boldsymbol{\theta}}_{\text{ML}} 在有限样本下是有偏的——它系统性地低估方差分量。这一偏差的根源在于：ML 在最大化似然时"消耗"了 pp 个自由度用于估计固定效应 β\boldsymbol{\beta}，却未对此进行任何校正，相当于将估计的 β^\hat{\boldsymbol{\beta}} 当作已知真实值处理。一个经典例子是正态分布方差的最大似然估计 σ^ML2=1n∑i=1n(yi−yˉ)2\hat{\sigma}^2_{\text{ML}} = \frac{1}{n}\sum_{i=1}^{n}(y_i - \bar{y})^2，其分母为 nn 而非 n−1n-1，恰好低估了方差。在更复杂的混合模型中，这一偏差会传导至所有方差分量，导致假设检验中的Ⅰ类错误率膨胀和置信区间覆盖率不足。

数学原理：误差对比变换

REML 的核心思想是找到一个线性变换，将原始观测 y\mathbf{y} 转化为一组误差对比（error contrasts），即其分布不再依赖于固定效应 β\boldsymbol{\beta}，仅依赖于方差参数 θ\boldsymbol{\theta}。具体地，寻找矩阵 An×(n−p)\mathbf{A}_{n \times (n-p)} 使得 A′X=0\mathbf{A}'\mathbf{X} = \mathbf{0} 且 A′A\mathbf{A}'\mathbf{A} 满秩。经过变换：

z=A′y=A′Xβ+A′Zu+A′ϵ=A′Zu+A′ϵ\mathbf{z} = \mathbf{A}'\mathbf{y} = \mathbf{A}'\mathbf{X}\boldsymbol{\beta} + \mathbf{A}'\mathbf{Z}\mathbf{u} + \mathbf{A}'\boldsymbol{\epsilon} = \mathbf{A}'\mathbf{Z}\mathbf{u} + \mathbf{A}'\boldsymbol{\epsilon}

由于 A′X=0\mathbf{A}'\mathbf{X} = \mathbf{0}，固定效应被完全消除。变换后的数据 z\mathbf{z} 服从分布：

z∼N(0,A′(ZGZ′+R)A)\mathbf{z} \sim N\left(\mathbf{0}, \mathbf{A}'(\mathbf{Z}\mathbf{G}\mathbf{Z}' + \mathbf{R})\mathbf{A}\right)

其分布仅依赖于 θ\boldsymbol{\theta}，与 β\boldsymbol{\beta} 无关。

基于 z\mathbf{z} 的似然函数称为限制似然函数（restricted likelihood），其对数形式为：

ℓR(θ;y)=−12log⁡∣V∣−12log⁡∣X′V−1X∣−12(y−Xβ^)′V−1(y−Xβ^)\ell_R(\boldsymbol{\theta}; \mathbf{y}) = -\frac{1}{2}\log|\mathbf{V}| - \frac{1}{2}\log|\mathbf{X}'\mathbf{V}^{-1}\mathbf{X}| - \frac{1}{2}(\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}})'\mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}})

其中 V=ZGZ′+R\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}' + \mathbf{R} 为边际方差-协方差矩阵，β^=(X′V−1X)−1X′V−1y\hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{X}'\mathbf{V}^{-1}\mathbf{y} 为广义最小二乘（GLS）估计量。与标准 ML 的对数似然相比，REML 多出了惩罚项 −12log⁡∣X′V−1X∣-\frac{1}{2}\log|\mathbf{X}'\mathbf{V}^{-1}\mathbf{X}|，该惩罚项恰好校正了由于估计 β\boldsymbol{\beta} 而损失的自由度，是 REML 消除偏差的关键所在。

值得强调的是，REML 变换矩阵 A\mathbf{A} 的选取并不唯一——任何满足 A′X=0\mathbf{A}'\mathbf{X} = \mathbf{0} 的满秩矩阵均产生相同的限制似然函数（最多相差一个乘法常数），因此 REML 估计具有变换不变性。常见的构造方式包括利用 X\mathbf{X} 的零空间基底或使用残差形成矩阵 I−X(X′X)−1X′\mathbf{I} - \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'。

REML估计量的统计性质

REML 估计量具有一系列优良的统计性质，使其成为方差分量估计的首选方法。

无偏性。在平衡设计（balanced design）下，REML 估计量是方差分量的最小方差无偏估计（MVU），即其期望精确等于真实方差参数值。在非平衡设计下，REML 虽然不再严格无偏，但其偏差远小于 ML，且随样本量增大迅速收敛至零。

一致性。随着样本量和组数（cluster 数）的增加，REML 估计量依概率收敛于真实参数值。一致性在固定效应维度 pp 固定、组数趋于无穷的标准渐近框架下成立。当 pp 也随样本增长时（高维情形），需借助Profile REML或修正的 REML 方法保持一致性。

渐近正态性。在适当正则条件下，REML 估计量具有渐近正态分布：

θ^REML∼⋅N(θ,IREML−1(θ))\hat{\boldsymbol{\theta}}_{\text{REML}} \stackrel{\cdot}{\sim} N\left(\boldsymbol{\theta}, \mathbf{I}_{\text{REML}}^{-1}(\boldsymbol{\theta})\right)

其中 IREML\mathbf{I}_{\text{REML}} 为 REML 的Fisher信息矩阵，通常通过期望信息矩阵（Expected Information）或观测信息矩阵（Observed Information）估计。这一性质为Wald检验和置信区间的构造提供了理论基础。

不变性。REML 估计量对固定效应 β\boldsymbol{\beta} 的重新参数化保持不变。换言之，无论采用何种对比编码（contrast coding）方案设定 X\mathbf{X}，方差分量的 REML 估计结果均一致。

REML与ML的比较

REML 与 ML 的选择取决于推断目标。

  • 方差分量估计：REML 显著优于 ML。REML 校正了自由度损失，提供近乎无偏的方差分量估计；ML 在有限样本下系统性低估方差分量，在小样本或复杂随机效应结构下偏差尤为严重。
  • 固定效应推断：ML 和 REML 均可用于固定效应估计，因为固定效应的 GLS 估计量 β^\hat{\boldsymbol{\beta}} 在两者的框架下形式一致。但对于固定效应的假设检验，REML 下的方差估计更准确，因此t检验和F检验的Ⅰ类错误率更接近名义水平。
  • 模型选择：此处歧见较大。若使用似然比检验（LRT）比较具有不同固定效应结构的模型，必须使用 ML 而非 REML。原因在于 REML 对数据进行了消除固定效应的变换——当两个模型的固定效应结构不同时，它们对应的变换矩阵 A\mathbf{A} 不同，限制似然函数定义在不同的尺度上，不可直接比较。因此，LRT 比较不同固定效应结构时应使用 ML 似然；比较不同随机效应结构（相同固定效应）时，REML 和 ML 的 LRT 均可使用，但 REML 通常更稳健。信息准则（如 AIC、BIC）用于 REML 下的模型比较时须格外谨慎，建议使用 ML 版本的 AIC/BIC 进行固定效应部分的模型选择。

计算实现与算法

REML 的实际计算通常不通过显式构造变换矩阵 A\mathbf{A}，而是通过数值优化算法最大化限制似然函数。最常用的方法包括：

Fisher Scoring 算法。基于期望信息矩阵进行迭代更新，收敛速度快且稳定性好，是多数统计软件的默认方法。迭代公式为：

θ(t+1)=θ(t)+IREML−1(θ(t))⋅S(θ(t))\boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} + \mathbf{I}_{\text{REML}}^{-1}(\boldsymbol{\theta}^{(t)}) \cdot \mathbf{S}(\boldsymbol{\theta}^{(t)})

其中 S(θ)\mathbf{S}(\boldsymbol{\theta}) 为 REML 的得分函数。

EM算法（期望最大化）。将随机效应 u\mathbf{u} 视为缺失数据，通过 E 步计算 u\mathbf{u} 的条件期望和条件方差，M 步更新方差分量估计。EM 算法收敛较慢但数值稳定，特别适用于复杂分层模型和广义线性混合模型（GLMM）的方差分量估计。实践中常采用ECME算法（Expectation/Conditional Maximization Either）加速收敛。

Newton-Raphson 算法。直接使用观测信息矩阵进行优化，收敛速度最快但对初值敏感。在高维问题中，计算 Hessian 矩阵的成本较高。

平均信息算法（Average Information, AI）。由 Gilmour、Thompson 和 Cullis 于 1995 年提出，使用期望信息矩阵和观测信息矩阵的平均值替代 Hessian，兼具 Fisher Scoring 的稳定性和 Newton-Raphson 的速度，是ASReml等专业遗传评估软件的核心算法。

主流统计软件对 REML 均有完善支持：R 语言的 \texttt{lme4} 包（\texttt{lmer} 函数默认使用 REML）、\texttt{nlme} 包的 \texttt{lme} 函数、SAS 的 \texttt{PROC MIXED}（默认 REML）、Stata 的 \texttt{mixed} 命令（默认 REML），以及 Python 的 \texttt{statsmodels} 中的 \texttt{MixedLM} 类。

扩展与前沿

REML 的经典框架已延伸至多个方向。REML 得分检验允许在无需拟合完整备择模型的情况下检验随机效应的显著性，特别适用于全基因组关联研究（GWAS）中海量遗传标记的快速筛选。贝叶斯 REML将先验分布引入方差分量，结合马尔可夫链蒙特卡洛（MCMC）方法进行后验推断，在数量遗传学的育种值预测中尤为流行。惩罚 REML（Penalized REML）在高维数据（p≫np \gg n）下通过对固定效应施加LASSO或岭回归惩罚，同时实现变量选择和方差分量估计。此外，稳健 REML通过替换正态假设为t分布或偏正态分布抵御离群值和重尾分布的影响。

记忆口诀与注意事项

记忆口诀：ML 估方差偏低，REML 做个投影——消掉固定效应再算，无偏一致最优选；比固定效应结构用 ML，比随机效应用 REML，混淆两者失之毫厘谬以千里。

使用 REML 时需注意：①REML 应用于比较不同固定效应结构的模型时必须切换至 ML；②REML 不允许直接比较不同变换空间下的似然值，如不同响应变量变换后的模型；③方差分量估计可能触及参数空间的边界（如方差估计为零），此时需使用约束优化或边界校正技术；④小样本下 REML 的无偏性依赖于正态性假设，偏离正态时仍需借助稳健方法或Bootstrap进行校正。

返回百科索引