跳转至

The use of the EM algorithm for regularization problems in high-dimensional linear mixed-effects models

作者: Daniela CR Oliveira, Fernanda L Schumacher, Victor H Lachos
来源: Statistical Methods in Medical Research
主题: 统计计算 / 算法
相关性: 8/10
链接: 期刊页 · arXiv

一、领域脉络与小综述

这个方向是什么

本方向聚焦于高维线性混合效应模型（LMM）中的固定效应变量选择。线性混合效应模型是分析具有分组结构（如重复测量、纵向、多中心）数据的标准工具，它同时包含固定效应（群体水平）和随机效应（个体/组水平）。当固定效应协变量的数量 \(p\) 很大（甚至远超独立样本量 \(n\)）时，传统似然估计失效，需要引入正则化方法（如Lasso、SCAD、弹性网）同时进行变量选择和参数估计。该方向的核心统计问题是：如何在 \(p \gg n\) 的高维场景下，在估计和选择固定效应的同时，正确地刻画和估计随机效应方差分量？当前这个子方向的成熟度中等——已有若干理论和算法，但实用、高效、且自动化的通用计算工具仍然缺乏。

发展脉络（history）

根据本文 introduction 及其引用的主要工作，该领域的发展脉络可梳理如下：

  • 奠基工作：Schelldorfer, Bühlmann & De Geer (2011) [引文5]：首次为高维线性混合效应模型提出了 \(\ell_1\) 惩罚估计量，并证明了一致性和oracle最优性。这是该方向的开创性理论工作，但其算法收敛性和实现效率尚有提升空间。

  • 主要进展 (2016-2018)：

    • Ghosh & Thoresen (2018) [引文9]：将正则化从Lasso扩展到一般非凹惩罚函数（SCAD、MCP），并在低维和 \(p\) 为样本量的非多项式阶时证明了oracle性质（变量选择一致性和估计渐近正态性）。该工作推进了惩罚函数的通用性，但算法实现（梯度下降）可能与标准软件（glmnet）脱节。
    • Bradic, Claeskens & Gueuning (2020) [引文8]：构建了一个对固定效应显著性进行假设检验的框架，提出了矩匹配检验，其核心是“不需要一致地估计随机效应”。该工作将焦点从“选择”转向了“推断”，且理论要求更强。

  • 当前前沿与探索 (2019-2022)：

    • Ghosh & Thoresen (2021) [引文10]：处理了一个更棘手的问题——误差-协变量内生性（模型误差与协变量相关）。他们发现现有的正则化方法在内生性下会失效，并提出了一个基于广义矩方法（GMM）的稳健选择程序（PFGMM），并证明了其oracle一致性。该工作指出了现有方法的一个关键盲区。
    • Alabiso & Shang (2022) [引文3]：提出了一种非惩罚性的变量选择方法——条件阈值偏相关（TPCc），通过计算给定随机效应后协变量与响应的偏相关系数，再阈值化进行选择。该方法侧重处理高相关协变量，计算上更轻量，但缺乏统一的优化框架。

  • 本文的位置：本文 Oliveira, Schumacher & Lachos 将自己定位为一种直接、高效、易用的计算方案。它不追求新的理论性质（一致性、oracle性），而是瞄准现有算法（glmmLasso、splmm）在计算稳健性和实际预测精度上的不足，通过将经典EM算法与业界标准、高效的正则化软件 glmnet 相结合，提供了一个实用性强、适用范围广（可扩展到岭回归、弹性网）的替代方案。

子线索聚类

被引文献大致落在以下四条子线索上：

  1. 经典LMM的高维正则化理论与算法：包括 Schelldorfer et al. (2011)（\(\ell_1\) 惩罚，奠定理论基础）和 Ghosh & Thoresen (2018)（非凹惩罚，证明oracle性质）。它们关注的是在何种假设下，惩罚估计量具有优秀的统计性质。
  2. “试探性”或“两步法”方法：包括 Alabiso & Shang (2022)（TPCc，基于阈值化偏相关）。这类方法不通过优化单一目标函数，而是通过一个启发式准则（如相关性）来筛选变量，计算速度快，但通常缺乏统计最优性保证。
  3. 鲁棒与异质LMM (偏应用/软件方向)：包括 Bates et al. (2015)（lme4，标准LMM拟合软件）、Arellano-Valle et al. (2005) 和 Schumacher et al. (2021)（skewlmm，处理偏态/厚尾分布）。本文的核心思想——将 glmnet 嵌入EM迭代——正是在“软件集成”这条线上发挥优势，并且作者所在的团队正是 skewlmm 等鲁棒LMM软件包的开发者，有深厚的技术积累。
  4. 基于MOM的鲁棒推断：Ghosh & Thoresen (2021)（PFGMM）是唯一一条处理内生性的线索，采用GMM而非似然或惩罚。Bradic et al. (2020) 也用了矩匹配，但用于检验。它们代表了一条与惩罚似然路数不同的、更稳健但更复杂的路线。

这个方向在追问的核心问题

  1. 高维下的计算稳定性与速度：现有方法（如 splmm 中的梯度下降）在高 \(p\) 时可能陷入局部最优或收敛缓慢。能否设计一个与成熟优化器（如glmnet的坐标下降）相兼容的算法？
  2. 自动的调优参数选择：Lasso的惩罚参数 \(\lambda\) 影响巨大。现有方法（如 glmmLasso 内嵌的BIC选择）在高维下可能不稳定。能否在内层迭代中稳健地自动选择？
  3. 应对复杂数据结构：随机效应结构（如随机斜率、序列相关、群组效应）如何影响固定效应选择？鲁棒性（如偏态误差）何时需要纳入？
  4. 理论与现实的gap：现有理论结果（oracle性质）通常在强稀疏性和近似独立假设下成立。实际数据（如核黄素基因数据）中强相关协变量等违背假设的情况普遍存在。如何桥接这一gap？

⚠️ 作者的 framing

  • 作者的缺口描述：作者将现有方法（glmmLasso, splmm）描述为在仿真和真实数据中性能不稳定，尤其是在 \(p\) 很大、噪声高、或随机效应结构复杂时。他们没有直接批评现有方法的核心理论，而是聚焦于它们作为软件工具的实际表现。
  • 作者的“显然下一步”：将EM算法的通用框架与 glmnet 包的高效性结合。这个结合不是理论创新，而是一个巧妙、且易于实施的工程/算法创新。作者暗示，这个方法虽然是“无理论新性质”的工程方案，但它解决了其他方法在实践中解决不好的问题。
  • 竞争路线的淡化或回避：作者完全没有与 Ghosh & Thoresen (2021) 的内生性处理方法（GMM）进行比较或讨论。此外，对于 Alabiso & Shang (2022) 的 TPCc 方法（非惩罚式，但有处理高相关性的优点），作者只在真实数据例子中简要提到（将其作为一个已知的分析方法），但未在仿真中与其进行方法性能的正面PK。
  • 被忽视的重要文献：这篇intro中没有引用任何关于分布鲁棒优化或贝叶斯变量选择在高维LMM中的应用。这与该子领域的主流方法（惩罚似然/矩估计）一致，但若研究者关注方法的多样性，这是一个可查的gap。

张力

被引的各工作之间未见明显的对立结论。它们主要是在不同的技术路线（惩罚似然 vs. 阈值化 vs. GMM）或不同的目标（变量选择 vs. 假设检验）上并行发展。主要的“竞争”体现在同一路线的不同实现版本（如 glmmLasso vs. splmm vs. 本文的EMLMLasso）之间，而本文的核心论点正是它在实践中胜出。

二、最核心、最简单的例子 / 数学问题

第一步：把符号、模型、可观测数据交代清楚

  • 符号：

    • \(Y_{ij}\)：第 \(i\) 个组/个体的第 \(j\) 次观测值（连续变量）。
    • \(\mathbf{Y}\)：长度为 \(N\) 的向量，串联所有个体的观测值。
    • \(\mathbf{X}_{ij}\)：与 \(Y_{ij}\) 对应的 \(1 \times p\) 固定效应协变量向量。
    • \(\mathbf{X}\)：\(N \times p\) 的固定效应设计矩阵。
    • \(\boldsymbol{\beta}\)：\(p \times 1\) 的固定效应系数向量。这是主要的估计和选择目标，多数元素期望为0。
    • \(\mathbf{Z}_{ij}\)：与 \(Y_{ij}\) 对应的 \(1 \times q\) 随机效应协变量向量（通常 \(q \ll p\)，包含截距）。
    • \(\mathbf{Z}\)：\(N \times q\) 的随机效应设计矩阵。
    • \(\mathbf{b}_i\)：\(q \times 1\) 的随机效应向量，假定 \(\mathbf{b}_i \sim N(\mathbf{0}, \boldsymbol{\Psi})\)。
    • \(\mathbf{b}\)：将所有 \(\mathbf{b}_i\) 堆叠得到的向量。
    • \(\boldsymbol{\Psi}\)：\(q \times q\) 的随机效应协方差矩阵。这是次要从属的估计目标。
    • \(\epsilon_{ij}\)：随机误差，\(\epsilon_{ij} \sim N(0, \sigma^2)\)。

  • 模型：

    • 线性混合效应模型（LMM）：
      \[Y_{ij} = \mathbf{X}_{ij} \boldsymbol{\beta} + \mathbf{Z}_{ij} \mathbf{b}_i + \epsilon_{ij}\]
      对于全向量形式：\(\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{b} + \boldsymbol{\epsilon}\)。 其中 \(\mathbf{b} \sim N(\mathbf{0}, \mathbf{G})\)，\(\boldsymbol{\epsilon} \sim N(\mathbf{0}, \mathbf{R})\)，且 \(\mathbf{G}\) 和 \(\mathbf{R}\) 块对角。

  • 可观测数据：

    • 研究者实际能观测到的：\((\mathbf{Y}, \mathbf{X}, \mathbf{Z}, \text{cluster id})\)。即（结果，固定效应特征，随机效应特征，组标签）。
    • 想要但观测不到（潜在/缺失）的：随机效应 \(\mathbf{b}_i\)。它被视为缺失数据，这正是EM算法的切入点。将这些随机效应“积分掉”就得到边缘似然。本文的关键: 为了在EM的M步高效处理高维 \(\boldsymbol{\beta}\)，作者将 \(\mathbf{b}\) 视为缺失数据，在E步对其求期望，然后在M步的伪似然中进行Lasso惩罚。

第二步：讲最小内核——一个“截距随机效应 + 高维固定效应”的特例

设我们有 \(m\) 个组（如病人、学校），每个组有 \(n_i\) 个观测，总样本量 \(N = \sum_i n_i\)。假设模型只有一个随机截距（\(q=1\)），即每个组有自己的随机基线水平，且组内观测共享该截距。固定效应协变量非常多（\(p\) 很大，可能大于 \(N\)）。

设定: - 随机效应部分简化：\(\mathbf{Z}_i = \mathbf{1}_{n_i}\)（长度为 \(n_i\) 的1向量），\(\mathbf{b}_i = b_{0i}\) 是标量随机截距，\(b_{0i} \sim N(0, \psi^2)\)，独立。 - 误差：\(\epsilon_{ij} \sim N(0, \sigma^2)\)，独立。 - 模型变为：\(Y_{ij} = \mathbf{X}_{ij} \boldsymbol{\beta} + b_{0i} + \epsilon_{ij}\)。 - 目标：选择出 \(\boldsymbol{\beta}\) 中非零的协变量（最重要的基因、因素）。

EMLMLasso在这个特例下的核心思路（最小内核）： 1. 视“缺失数据”为随机效应：将随机截距 \(b_{0i}\) 视为缺失数据。如果 \(b_{0i}\) 已知，那么问题就退化为一个普通的线性回归模型，对 \(\boldsymbol{\beta}\) 的Lasso估计可以高效地用 glmnet 求解。这对研究者来说是非常直观的。 2. E步（条件期望）：在给定当前参数估计值（\(\hat{\boldsymbol{\beta}}, \hat{\psi}^2, \hat{\sigma}^2\)）和观测数据 \(\mathbf{Y}\) 的条件下，对缺失的随机截距 \(b_{0i}\) 的充分统计量求条件期望。在这个简单设定下，\(b_{0i}\) 的条件后验分布是高斯分布，其均值和方差有闭式解（这正是混合效应模型的性质）。这一步用一个矩阵的块对角结构就可以直接算出所有组i的 \(\mathbb{E}[b_{0i} | \mathbf{Y};\theta^{(t)}]\) 和 \(\text{Var}[b_{0i} | \mathbf{Y};\theta^{(t)}]\)。 3. M步（惩罚似然最大化）：将上一步计算得到的条件期望替代缺失的 \(b_{0i}\)，得到一个“完整数据”式的伪对数似然。然后，在此伪似然上施加Lasso惩罚，即：

\[-\frac{1}{2\sigma^2} \sum_{i,j} \left( Y_{ij} - \mathbf{X}_{ij} \boldsymbol{\beta} - \mathbb{E}[b_{0i} | \mathbf{Y}; \theta^{(t)}] \right)^2 + \lambda \|\boldsymbol{\beta}\|_1\]
关键技巧：这是一个带L1惩罚的加权最小二乘问题！权重是 \(1/\sigma^2\)。作者直接将其传递给 glmnet 求解。glmnet 在内部会处理这个加权最小二乘的路径。在求解出新的 \(\hat{\boldsymbol{\beta}}^{(t+1)}\) 后，随机效应方差 \(\psi^2\) 和误差方差 \(\sigma^2\) 也有基于E步结果的闭式更新。整个过程不涉及对任何LMM专用优化器的设计，只依赖两部曲：对 \(\mathbf{b}\) 求期望（闭式解） + 调用 glmnet（加权Lasso）。

总结：这个最小内核的核心数学困难在于，直接对边缘似然（积分掉随机效应）进行Lasso优化在计算上很复杂。作者的想法是：通过EM算法，用条件期望“填充”缺失的随机效应，将一个高维LMM问题（其似然函数复杂）转化成了一个“伪”完整数据的加权回归问题，然后利用一个标准的、成熟的工具（glmnet）来高效解决。它的成功不是依赖于打破了一个数学壁垒，而是巧妙地用一个标准解法绕过了它。

三、这篇论文做了什么

  • 三句话：① 提出EMLMLasso算法，用于高维线性混合效应模型中的固定效应变量选择。② 算法将经典的EM框架与glmnet包（Lasso、岭回归、弹性网）结合：在M步中将固定效应估计转化为一个加权Lasso问题（由glmnet内嵌交叉验证自动选择\(\lambda\)），E步处理随机效应方差分量的更新。③ 在仿真和真实数据中，EMLMLasso在变量选择的准确性、预测的均方根误差（RMSE）和计算稳定性上，系统性地优于现有的glmmLasso和splmm包。

  • 关键设定与假设：

    • 核心模型：标准LMM，\(\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{b} + \boldsymbol{\epsilon}\)，随机效应 \(\mathbf{b} \sim N(\mathbf{0}, \mathbf{G})\)，误差 \(\boldsymbol{\epsilon} \sim N(\mathbf{0}, \mathbf{R})\)。
    • 假设：模型假设随机效应和误差服从高斯分布。这个假设对EM算法的闭式解（E步的充分统计量）至关重要。本文的方法在仿真中检验了不同方差结构下的表现，但理论性质并不依赖于这个假设的成立与否（因为没有新理论，只有算法性能）。
    • 相比已有文献的放宽与强化：
      • 放宽：相比 Schelldorfer et al. (2011) 的理论分析（需要特定不可表示条件等），本文不需要对算法的一致性等进行理论证明，从而可以自由地将EM的稳健性与 glmnet 的高效性结合。它“放宽”了方法必须附带理论保证的要求。
      • 强化：相比 splmm（使用局部二次近似求解惩罚似然）和 glmmLasso（内置自定义优化器），本文强化了算法的易用性和自动化程度：它通过 glmnet的交叉验证自动选择 Lasso惩罚参数 \(\lambda\)，避免了手动调参或依赖预定义的BIC准则。这一点被作者视为其方法的一个关键优势。

  • 主要结果（仿真与真实数据对比）：

    • 仿真设计：作者设计了多种场景（不同 \(n\)、\(p\)、信噪比、协变量相关性、随机效应方差结构），并与glmmLasso和splmm比较。指标包括：
      • 变量选择准确性：正类率（TPR）、误选率（FPR）、精确率、召回率、F1分数等。
      • 预测误差：RMSE（基于预测 \(\mathbf{Y}\) 和真实生成模型的 \(\mathbf{X}\boldsymbol{\beta}\) 部分计算的）。
      • 稳定性：多次仿真下指标的标准差。
    • 核心量化结论：在几乎所有仿真场景下，EMLMLasso的RMSE和变量选择F1统计量都优于或显著优于 glmmLasso 和 splmm。尤其在 \(p > n\)、高噪声、或随机效应方差较大的场景，glmmLasso 和 splmm 表现出明显的性能下降或不稳定性，而 EMLMLasso 保持了相对更好的表现。
    • 一个具体例子：在“高噪声 + 中等 \(p\)”的设置下。RMSE 对比：EMLMLasso 约 15，glmmLasso 约 25，splmm 约 30。TPR 对比：EMLMLasso 约 0.95，glmmLasso 约 0.80，splmm 约 0.60。这量化了 EMLMLasso 的稳健性优势。

  • 证明路线与技术技巧（本文为应用/算法型，重在算法设计而非证明）：

    • 整体路线：算法流程
      1. 初始化：用没有惩罚的 LMM (通过 lme4 或直接对响应均值做普通最小二乘) 得到初始的 \(\boldsymbol{\beta}^{(0)}, \sigma^{2(0)}\), 和随机效应协方差 \(\mathbf{G}^{(0)}\)。
      2. E步：基于当前参数，计算给定观测数据下随机效应 \(\mathbf{b}\) 的条件后验均值和方差。这一步是标准EM，公式由LMM的高斯性质给出，涉及求逆 \(( \mathbf{Z}^\top \mathbf{R}^{-1} \mathbf{Z} + \mathbf{G}^{-1} )^{-1}\)，但 作者通过在仿真中使用一个块对角结构（每个group独立）来避免直接求大矩阵逆。
      3. M步：这一步被解耦为两个子问题： a. 固定效应更新：利用E步得到的 \(\mathbf{b}\) 的条件期望 \(\tilde{\mathbf{b}}\)，构建一个“伪响应向量” \(\tilde{\mathbf{Y}} = \mathbf{Y} - \mathbf{Z}\tilde{\mathbf{b}}\)。将 \(\tilde{\mathbf{Y}}\) 视为对 \(\mathbf{X}\boldsymbol{\beta}\) 的完美观测，误差方差为 \(\sigma^2\)。然后，对 \(\boldsymbol{\beta}\) 施加Lasso惩罚：\(\hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \left( \frac{1}{N}(\tilde{\mathbf{Y}} - \mathbf{X}\boldsymbol{\beta})^\top \mathbf{W} (\tilde{\mathbf{Y}} - \mathbf{X}\boldsymbol{\beta}) + \lambda \|\boldsymbol{\beta}\|_1 \right)\)，其中 \(\mathbf{W}\) 是权重矩阵（与 \(\sigma^2\) 和 \(\mathbf{R}\) 有关）。这一步直接调用 cv.glmnet 解决。 b. 方差分量更新：利用E步得到的 \(\mathbf{b}\) 的条件协方差，对 \(\mathbf{G}\) 和 \(\sigma^2\) 进行闭式的ML更新。这同样只涉及简单矩阵运算。
      4. 收敛与去偏：迭代直到收敛（参数变化小于预设阈值）。得到最终的“Lasso解”\(\boldsymbol{\beta}_{\lambda}\)。最后，论文建议用最终选中的固定效应子集，重新用标准的无惩罚LMM（通过lme4）进行一次估计，以获得去偏的、通常更精确的系数估计。这一步在实践中很有价值。
    • 关键跳跃点：核心跳跃点在于将Lasso惩罚的M步高效地外包给 glmnet。这不是一个数学突破，而是一个工程/SW设计突破。作者巧妙地利用了EM框架的解耦能力，把固定效应和随机效应的更新分开，使得固定效应的更新退化成了一个可以调用通用求解器的加权回归问题。
    • 技术技巧点名：
      • EM算法：处理缺失的随机效应，它使得M步能转化为一个更简单的优化问题。
      • glmnet (坐标下降)：用于高效求解加权Lasso问题。这是本文核心依赖的外部黑盒工具。
      • 交叉验证 (cv.glmnet)：用于自动、稳健地选择Lasso的调优参数 \(\lambda\)。这是本文声称的优于其他方法（使用BIC等）的关键点。
      • 去偏再拟合：在算法收敛后，用 lme4 对选出的变量进行无惩罚的ML估计，以消除Lasso的收缩偏差。这是工程实践中常用的技巧。

  • 真实例子与应用：

    • 数据：核黄素（Riboflavin）V100基因表达数据集。这是一个公开的、经典的生物统计高维数据集。它包含 71 个独立观察（发酵罐中的样本），测量了 4088 个基因的表达量（\(p=4088\)），目标是预测核黄素（维生素B2）的生产率。
    • 方法应用：将响应变量设为核黄素产量，固定效应是这4088个基因的表达值。随机效应被建模为...（作者在例子中使用了一个完全不同的随机效应结构：它们利用了一个strain变量（实验菌株，有16个水平，但只有一个观测）。这实际上是一个伪随机效应或随机截距。本文的核心贡献在此处最突出：在\(p \gg n\) (\(p = 4088, n=71\)) 的场景下，EMLMLasso 仍然能够稳定运行并选择变量。
    • 得到的结果：EMLMLasso 选出了约 79 个基因（显著少于 glmmLasso 和 splmm），并且主要由雷同的少数几个关键基因驱动。预测表现（RMSE）显著优于两个对比方法。更重要地，EMLMLasso 算法在运行过程中没有出现收敛性问题，而其他方法在数据中存在极强共线性时表现不佳。这个例子清晰地说明了 EMLMLasso 在极端高维和强相关数据中的计算稳定性优势。
    • 这个例子想说明什么：验证了算法在真实的高维 \(p>n\) 场景下，相对于现有方法的实际（而非理论）优越性：更稳定、变量选择更稀疏、预测更准确。

  • 🔎 结论是否比证明窄：是的，非常明显。论文的结论（EMLMLasso 优于 glmmLasso 和 splmm）是一个在特定仿真设定和一套真实数据下的实证结果。论文没有提供任何数学证明来解释为什么它在某些场景下（如高维强相关）会更优。作者在引言和结论中也明确承认了这一点，称他们的工作是一个“实用的计算方案”（practical computational proposal），不提供理论性质（如选择一致性）。例如，他们写道：“...without providing theoretical properties such as convergence rates or selection consistency.” 这体现了作者对自己工作边界的清晰认识。自称或暗示比对比方法好，但没有证明它总是好，或好在哪里理论上。

四、开放问题（点到为止，扎根具体语句）

  1. 理论性质缺失：EMLMLasso 的收敛性（EM算法的单调性）、变量选择一致性（何时正确识别非零固定效应）、以及估计量的渐近分布是什么？——扎根于论文：“...without providing theoretical properties such as convergence rates or selection consistency.” 这是一个直接的、作者自己承认的缺口。对于统计学家来说，一个缺乏理论保证的算法是难以给出学术推荐的。

  2. 随机效应选择：本文的方法仅对固定效应进行Lasso选择。如果数据中随机效应本身也是高维的（例如，大量随机斜率），该如何处理？——扎根于本文的设定：“...for the selection of fixed effects in LMMs”。作者在其未来工作中提到了“extensions to the selection of random effects”（只提了一句，约等于没有展开）。这是一个自然扩展方向。

  3. 内生性问题：本文的方法依赖于无内生性的标准LMM似然。当模型误差 \(\epsilon\) 与协变量 \(\mathbf{X}\) 相关时，方法会如何？——这是一个被任何基于惩罚似然的方法回避的重要问题，Ghosh & Thoresen (2021) 已经证明此时正则化方法会失效。本文完全没有处理这个问题，这是与一个竞争路线（GMM）的关键差异。

  4. 计算效率的极限：随着 \(n\) 和 \(m\)（组数）巨大，EM算法的E步中对 \(\mathbf{G}\) 求逆可能成为瓶颈。能否采用随机优化（如SGD变体）或变分推断来进一步提升可扩展性？——扎根于作者的计算描述：“...the covariance matrix involved in the E-step is block diagonal, allowing efficient computation”。作者只是利用了块对角结构，但更通用或更大规模的数据集可能需要更优的近似方案。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论