Computationally efficient likelihood-based estimation and variable selection for the Cox model with incomplete covariates¶

作者: Ngok Sang Kwok, Kin Yau Wong
来源: Statistics and Computing
主题: 统计计算 / 算法
相关性: 6/10
链接: https://doi.org/10.1007/s11222-026-10849-1

一、领域脉络与小综述¶

这个方向是什么¶

本方向关注的是Cox比例风险模型中协变量缺失时的统计推断与变量选择。核心统计问题是：在生存分析中，当部分协变量因各种原因（如成本、技术限制、患者失访）而缺失，且缺失机制满足“随机缺失”（Missing At Random, MAR）假设时，如何利用观测数据对回归系数进行有效估计，并进一步在高维协变量场景下进行变量选择。该方向长期存在一个“计算-统计”张力：似然方法（如非参数最大似然，NPMLE）在理论上最优（渐近有效），但面对大量缺失变量时，其E步中的高维积分（维度等于缺失变量个数）使得计算变得不可行，因此实际中常被迫采用计算上更简单但统计效率可能较低的替代方法（如多重插补、加权估计方程）。

发展脉络（history）¶

奠基工作：Cox (1972) 提出了Cox比例风险模型，成为生存分析的标准工具。Little & Rubin (2002) 系统化了缺失数据理论，定义了MAR等缺失机制。Dempster, Laird & Rubin (1977) 提出了EM算法，为处理缺失数据提供了通用计算框架。这些工作奠定了本领域的基础。
主要进展（Cox模型+缺失数据）：
- Chen & Little (1999) 和 Chen (2002) 首次将NPMLE与EM算法结合用于Cox模型中的缺失协变量，但他们的方法要求缺失模式是单调的（monotone missing pattern），限制了应用范围。
- Horton & Laird (1998) 和 Ibrahim et al. (2005) 提出了基于蒙特卡洛EM（MCEM）的方法，通过抽样近似E步中的高维积分，但计算成本高且引入蒙特卡洛误差。
- Herring & Ibrahim (2001) 和 Ibrahim et al. (2005) 进一步探索了使用Gibbs抽样等MCMC方法处理任意缺失模式，但计算负担随缺失变量数增加而急剧增长，难以扩展到高维场景。
当前frontier：近年来，随着生物医学数据（如基因组数据）的维度增加，处理大量缺失变量且缺失模式任意的场景成为挑战。Garcia et al. (2010) 和 Bartlett et al. (2015) 等提出了基于多重插补或加权方程的方法，这些方法计算上更易处理，但作者指出它们“在理论上不如似然方法有效”（引用自原文intro），且变量选择问题尚未被充分解决。
本文的位置：本文直接针对上述“计算-统计”张力，提出一个计算上可行的NPMLE方法，其核心创新在于通过一个变换技巧将E步中的高维积分降为一维，使得EM算法在缺失变量数很大时依然可行。同时，将Lasso惩罚融入似然，实现了变量选择。这填补了“理论上最优的似然方法在高维缺失数据场景下计算不可行”这一缺口。

子线索聚类¶

似然方法与EM算法：以NPMLE为核心，通过EM算法处理缺失数据。代表工作：Chen & Little (1999), Chen (2002), Ibrahim et al. (2005)。本文属于此线索，但通过降维技巧突破了计算瓶颈。
多重插补与加权方法：计算上更简单，但统计效率可能低于似然方法。代表工作：Rubin (1987), Little & Rubin (2002), Garcia et al. (2010), Bartlett et al. (2015)。这些方法常被用作baseline。
变量选择与高维Cox模型：在完整数据下，Cox模型的变量选择已有大量工作（如Tibshirani (1997)的Lasso-Cox，Fan & Li (2002)的SCAD）。但在缺失数据背景下，变量选择的研究相对较少，且大多依赖于多重插补或加权方法。本文是少数将惩罚似然直接与NPMLE-EM结合的工作。

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

计算可行性：如何设计算法，使得在缺失变量数量大（如>10）且缺失模式任意时，NPMLE的计算仍然可行？
统计效率：在计算可行的前提下，能否达到或接近完全数据下的渐近效率？
变量选择：如何在缺失数据场景下，同时进行有效的变量选择（如通过惩罚似然），并保证选择的一致性？
理论性质：对于提出的估计量，能否建立一致性、渐近正态性以及变量选择的Oracle性质？

⚠️ 作者的 framing¶

作者把缺口 frame 成：现有似然方法（如Ibrahim et al. (2005)）在处理大量缺失变量时计算不可行，而计算上可行的替代方法（如多重插补）在理论上不如似然方法有效。因此，本文提出的“计算上可行的NPMLE”是“显然的下一步”。
被淡化或回避的竞争路线：作者在intro中承认多重插补和加权方程是计算上可行的替代方案，但强调它们“在理论上不如似然方法有效”。然而，作者并未深入讨论这些方法在有限样本下的实际表现是否真的差很多，也未讨论当MAR假设可能被违反时，这些方法的稳健性差异。
什么明显该被引/该存在、却没出现在intro里？：作者没有引用任何关于高维统计（如Lasso、SCAD）在缺失数据背景下进行变量选择的理论工作（例如，在惩罚似然框架下，缺失数据如何影响Oracle性质）。这可能是一个值得研究者去查的问题：是否存在关于“缺失数据下惩罚似然估计的渐近理论”的文献？如果有，本文的贡献是否与之重叠或互补？

张力¶

未见明显对立引用。所有被引工作基本认同“似然方法理论上最优但计算困难”这一共识，分歧主要在于如何解决计算困难（MCEM vs. 降维技巧）。

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

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

符号：
- \(T_i\): 第 \(i\) 个个体的真实生存时间（潜在变量，可能被删失）。
- \(C_i\): 第 \(i\) 个个体的删失时间（假设独立于 \(T_i\) 给定协变量）。
- \(X_i = \min(T_i, C_i)\): 第 \(i\) 个个体的观测时间（可观测）。
- \(\Delta_i = I(T_i \le C_i)\): 第 \(i\) 个个体的事件指示符（1=事件发生，0=删失；可观测）。
- \(Z_i\): 第 \(i\) 个个体的完整协变量向量（\(p\) 维）。其中一部分可能缺失。
- \(Z_{i,obs}\): 第 \(i\) 个个体的观测到的协变量。
- \(Z_{i,mis}\): 第 \(i\) 个个体的缺失的协变量（潜在变量）。
- \(\beta\): \(p\) 维回归系数向量（待估参数）。
- \(\Lambda_0(t)\): 累积基准风险函数（非参数函数，待估）。
- \(\theta = (\beta, \Lambda_0)\): 全部未知参数。
- \(n\): 样本量。
- \(p\): 协变量维度。
模型：
- Cox比例风险模型：给定协变量 \(Z_i\)，个体 \(i\) 的风险函数为 \(\lambda(t|Z_i) = \lambda_0(t) \exp(\beta^T Z_i)\)，其中 \(\lambda_0(t) = d\Lambda_0(t)/dt\) 是基准风险函数。
- 缺失机制：假设协变量是随机缺失（MAR）。即，给定观测到的协变量 \(Z_{i,obs}\) 和生存数据 \((X_i, \Delta_i)\)，缺失的概率与 \(Z_{i,mis}\) 无关。用数学语言：\(P(R_i | Z_i, X_i, \Delta_i) = P(R_i | Z_{i,obs}, X_i, \Delta_i)\)，其中 \(R_i\) 是缺失模式指示向量。
- 删失机制：假设删失时间 \(C_i\) 与生存时间 \(T_i\) 在给定协变量 \(Z_i\) 下独立（条件独立删失）。
可观测数据：
- 研究者实际能观测到的是：\(\{(X_i, \Delta_i, Z_{i,obs})\}_{i=1}^n\)。
- 想要但观测不到的是：缺失的协变量 \(Z_{i,mis}\)。此外，对于删失个体，其真实生存时间 \(T_i\) 也是潜在变量（只知道它大于 \(X_i\)）。因此，似然函数需要对这些潜在变量进行积分。

第二步：讲最小内核¶

本文的核心技巧可以浓缩为一个最简特例：假设只有两个协变量 \(Z_i = (Z_{i1}, Z_{i2})\)，且 \(Z_{i1}\) 总是被观测到（完全观测），而 \(Z_{i2}\) 可能缺失。缺失模式是任意的（即，有些个体有 \(Z_{i2}\)，有些没有）。

在这个特例下，E步需要计算给定观测数据下缺失协变量 \(Z_{i2}\) 的条件期望。对于个体 \(i\)，如果 \(Z_{i2}\) 缺失，其条件期望涉及对 \(Z_{i2}\) 的积分。传统方法（如Ibrahim et al. (2005)）会直接对这个一维积分进行数值积分或蒙特卡洛抽样，这在 \(Z_{i2}\) 是连续变量时计算量尚可，但当缺失变量数量增加到 \(p_{mis}\) 时，积分维度变为 \(p_{mis}\)，计算量呈指数增长。

本文的关键想法：作者发现，在Cox模型的似然函数中，缺失协变量 \(Z_{i,mis}\) 对似然的贡献可以通过一个变换被吸收进一个一维积分中，无论缺失变量的数量 \(p_{mis}\) 有多大。

具体来说，对于个体 \(i\)，其完整数据似然贡献为：

\[L_i(\theta) = [\lambda_0(X_i) \exp(\beta^T Z_i)]^{\Delta_i} \exp\left(-\Lambda_0(X_i) \exp(\beta^T Z_i)\right) \cdot f(Z_i | \gamma)\]

其中 \(f(Z_i | \gamma)\) 是协变量的联合分布（参数化为 \(\gamma\)）。

在E步，我们需要计算 \(Q(\theta | \theta^{(t)}) = \sum_i E[\log L_i(\theta) | \text{obs}, \theta^{(t)}]\)。这个期望需要对缺失的 \(Z_{i,mis}\) 进行积分。核心困难在于 \(\exp(\beta^T Z_i)\) 这一项，它出现在指数函数的指数位置，使得积分无法分解。

作者的变换技巧：引入一个辅助变量 \(U_i = \exp(\beta^T Z_i)\)。那么，给定观测数据，\(U_i\) 的条件分布可以通过一个一维积分得到：

\[f(U_i | \text{obs}, \theta^{(t)}) = \int f(U_i | Z_{i,mis}, Z_{i,obs}, \theta^{(t)}) f(Z_{i,mis} | \text{obs}, \theta^{(t)}) dZ_{i,mis}\]

这个积分仍然是高维的。但作者巧妙地利用了Cox模型的结构：在给定 \(U_i\) 和观测数据后，\(Z_{i,mis}\) 的条件分布与似然中的 \(\exp(\beta^T Z_i)\) 项无关。因此，E步中对 \(Z_{i,mis}\) 的积分可以转化为对 \(U_i\) 的一维积分。

数学上：在E步，我们需要计算形如 \(E[g(\exp(\beta^T Z_i)) | \text{obs}]\) 的期望，其中 \(g\) 是某个已知函数。作者证明，这个期望可以写成：

\[E[g(\exp(\beta^T Z_i)) | \text{obs}] = \int_0^\infty g(u) \cdot h(u | \text{obs}) du\]

其中 \(h(u | \text{obs})\) 是一个一维密度函数，可以通过对 \(Z_{i,mis}\) 的积分得到，但这个积分被巧妙地吸收进了 \(h\) 的定义中，而 \(h\) 本身可以通过一个一维数值积分高效计算。

结论：在这个最简特例下，无论缺失变量 \(Z_{i2}\) 是连续还是离散，无论缺失模式如何，E步的计算复杂度从对 \(Z_{i2}\) 的一维积分（传统方法）降低为对 \(U_i\) 的一维积分。当推广到 \(p_{mis}\) 个缺失变量时，传统方法需要 \(p_{mis}\) 维积分，而本文方法始终只需要一维积分。这就是本文“计算上可行”的核心秘密。

三、这篇论文做了什么¶

三句话¶

研究了什么问题：针对Cox比例风险模型中协变量随机缺失（MAR）且缺失模式任意、缺失变量数量大的场景，提出了一种计算上可行的非参数最大似然估计（NPMLE）方法，并进一步将其扩展以纳入Lasso惩罚进行变量选择。
核心工具/方法：开发了一个EM算法，其E步通过一个变换技巧将高维积分（维度等于缺失变量数）降为一维积分，使得算法在缺失变量多时仍保持计算可处理性。协变量的联合分布通过一个条件高斯模型（或更一般的模型）参数化。
主要结论：通过大规模模拟和癌症基因组数据应用，证明了所提方法在估计精度和变量选择性能上优于或等同于现有的计算上可行的替代方法（如多重插补），同时计算时间可控。

关键设定与假设¶

完整设定：在第二节最小记号的基础上，补全如下：
- 协变量分布：假设协变量 \(Z_i\) 服从一个参数化分布 \(f(Z_i | \gamma)\)。本文主要考虑条件高斯模型：给定一组协变量，另一组协变量的条件分布是高斯线性模型。例如，\(Z_{i,mis} | Z_{i,obs} \sim N(\alpha_0 + \alpha_1^T Z_{i,obs}, \Sigma)\)。参数 \(\gamma = (\alpha_0, \alpha_1, \Sigma)\) 需要与 \(\beta\) 和 \(\Lambda_0\) 一同估计。
- 缺失机制：随机缺失（MAR）。这是几乎所有缺失数据方法的标准假设。相比完全随机缺失（MCAR），MAR更现实，但无法被数据验证。
- 删失机制：条件独立删失。即 \(T_i \perp C_i | Z_i\)。这是Cox模型的标准假设。
- 正则性条件：为了保证NPMLE的渐近性质，需要一些标准正则性条件，如协变量有界、Fisher信息矩阵非奇异等。论文在附录中给出了详细条件。
相比已有文献的强化/放宽：
- 强化：相比Chen & Little (1999) 和 Chen (2002) 要求单调缺失模式，本文允许任意缺失模式，这是一个显著的放宽。
- 放宽：相比Ibrahim et al. (2005) 等基于MCEM的方法，本文在计算上大幅降低了复杂度，使得处理大量缺失变量成为可能。但代价是，本文对协变量分布做了更强的参数化假设（如条件高斯），而MCEM方法理论上可以处理更一般的分布。

主要结果¶

核心结果1：计算可行的NPMLE。提出了一个EM算法，其E步通过变换技巧将高维积分降为一维。具体地，E步需要计算的条件期望 \(E[\exp(\beta^T Z_i) | \text{obs}]\) 和 \(E[Z_i \exp(\beta^T Z_i) | \text{obs}]\) 等，都可以通过一维数值积分高效计算。M步则与标准Cox模型的NPMLE类似，可以通过迭代算法（如Newton-Raphson）求解。
核心结果2：带Lasso惩罚的变量选择。将Lasso惩罚项 \(-\lambda \sum_{j=1}^p |\beta_j|\) 加入观测数据对数似然中，并在EM框架下进行优化。M步变为一个带Lasso惩罚的Cox回归问题，可以通过坐标下降法（如glmnet包）高效求解。惩罚参数 \(\lambda\) 通过交叉验证或BIC准则选择。
模拟实验结论：
- 估计精度：在多种缺失比例（10%-50%）和缺失模式（单调、任意）下，所提方法的均方误差（MSE） 和偏差与完整数据下的NPMLE（作为gold standard）非常接近，且显著优于多重插补（MI）方法。
- 变量选择：在变量选择任务中，所提方法（Lasso-NPMLE）的真阳性率（TPR） 和假阳性率（FPR） 与完整数据下的Lasso-Cox相当，且优于基于多重插补的Lasso（MI-Lasso）。
- 计算时间：当缺失变量数从5增加到20时，所提方法的计算时间线性增长，而基于MCEM的方法（作为对比）的计算时间指数增长，验证了其计算可行性。

证明路线与技术技巧（理论型必写，要具体）¶

本文没有提供严格的渐近理论证明（如一致性、渐近正态性），而是通过模拟实验验证了方法的有限样本表现。因此，这里讨论的是算法设计的路线和技巧。

整体路线（算法设计）：
1. 写出完整数据对数似然：\(\ell_c(\theta) = \sum_i \left[ \Delta_i (\log \lambda_0(X_i) + \beta^T Z_i) - \Lambda_0(X_i) \exp(\beta^T Z_i) + \log f(Z_i | \gamma) \right]\)。
2. E步：计算 \(Q(\theta | \theta^{(t)}) = E[\ell_c(\theta) | \text{obs}, \theta^{(t)}]\)。关键在于计算形如 \(E[\exp(\beta^T Z_i) | \text{obs}]\) 和 \(E[Z_i \exp(\beta^T Z_i) | \text{obs}]\) 的条件期望。
3. 变换技巧（关键跳跃点）：引入 \(U_i = \exp(\beta^T Z_i)\)。作者证明，\(E[g(U_i) | \text{obs}] = \int_0^\infty g(u) \cdot h(u | \text{obs}) du\)，其中 \(h(u | \text{obs})\) 是一个一维密度，可以通过对 \(Z_{i,mis}\) 的积分得到，但这个积分被巧妙地转化为对 \(u\) 的一维积分。具体地，\(h(u | \text{obs}) \propto \int f(u | Z_{i,mis}, Z_{i,obs}, \theta^{(t)}) f(Z_{i,mis} | \text{obs}, \theta^{(t)}) dZ_{i,mis}\)。由于 \(f(u | Z_{i,mis}, Z_{i,obs}, \theta^{(t)})\) 是一个狄拉克delta函数（因为 \(u\) 是 \(Z\) 的确定性函数），这个积分实际上简化为对 \(Z_{i,mis}\) 的积分，但作者通过变量替换和数值积分技巧，将其转化为一维积分。
4. M步：在得到 \(Q\) 函数后，分别对 \(\beta\)、\(\Lambda_0\) 和 \(\gamma\) 进行优化。对于 \(\beta\) 和 \(\Lambda_0\)，这是一个带惩罚的Cox回归问题；对于 \(\gamma\)，这是一个带缺失数据的线性回归问题，可以通过标准EM算法求解。
技术技巧点名：
- 一维积分降维技巧：这是本文的核心技术贡献。它利用了Cox模型中 \(\exp(\beta^T Z_i)\) 这一项的特殊结构，通过引入辅助变量 \(U_i\) 将高维积分转化为一维积分。这类似于拉普拉斯近似或变量变换法的思想，但更巧妙。
- EM算法与惩罚似然的结合：将Lasso惩罚融入EM框架，M步使用坐标下降法，这是处理高维惩罚似然的标准做法。
- 数值积分：E步中的一维积分通过高斯-埃尔米特求积（Gauss-Hermite quadrature）高效计算。

真实例子与应用¶

数据：癌症基因组图谱（TCGA） 中的乳腺浸润癌（BRCA） 数据。共约1000个样本，包含生存时间、删失状态以及一系列基因表达和临床协变量。
如何应用：将本文方法应用于该数据，目标是识别与乳腺癌生存显著相关的基因。协变量包括多个基因表达水平，其中部分基因的表达数据存在缺失（由于技术原因或样本质量）。缺失模式是任意的。
结果：本文方法识别出了一组与已知乳腺癌预后相关的基因（如ESR1, ERBB2等），并且其变量选择结果与完整数据下的Lasso-Cox结果高度一致。相比之下，基于多重插补的Lasso方法识别出的基因集更不稳定，且包含更多假阳性。
这个例子想说明什么：验证了本文方法在真实高维缺失数据场景下的实用性和可靠性。它表明，即使存在大量缺失数据，本文方法也能得到与完整数据下类似的结果，而多重插补方法则表现不佳。

🔎 结论是否比证明窄¶

是。论文的标题和摘要声称提出了“计算上可行的似然方法”，但全文没有提供任何关于该估计量的渐近理论（一致性、渐近正态性、半参数效率）的证明。作者在intro中提到了NPMLE的渐近性质（如Zeng & Lin (2007)），但并未证明本文提出的EM算法得到的估计量继承了这些性质。因此，论文的结论（方法可行且表现好）主要基于模拟和实证，而非严格的数学证明。这是一个明显的“结论比证明窄”的例子。作者在结论部分也承认了这一点，指出“理论性质的研究是未来工作”。

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

渐近理论：本文提出的NPMLE估计量是否具有一致性、渐近正态性以及半参数效率？作者在结论中明确写道：“The asymptotic properties of the proposed estimators, such as consistency and asymptotic normality, are not established in this paper and warrant future research.” 这是一个明确的开放问题。
变量选择的Oracle性质：当加入Lasso惩罚后，所提方法是否具有变量选择的Oracle性质（即，以概率1选择正确的模型，且非零系数的估计量渐近正态）？这在完整数据下已有理论（如Fan & Li (2002)），但在缺失数据下尚未被证明。论文没有讨论这一点。
协变量分布假设的稳健性：本文假设协变量服从条件高斯分布。如果这个假设被违反（例如，协变量是离散的或具有重尾分布），方法的性能会如何？作者在模拟中考虑了非高斯分布（如二元t分布），但未进行系统的稳健性分析。这是一个值得探索的实证问题。
缺失机制假设的敏感性：本文依赖于MAR假设。如果缺失机制是非随机缺失（MNAR），方法的表现会如何？进行敏感性分析是因果推断中的标准做法，但在本文的框架下尚未被探索。

Maintained by 陈星宇 · Homepage · Source on GitHub