跳转至

Penalized joint models of high-dimensional longitudinal biomarkers and a survival outcome

作者: Jiehuan Sun, Sanjib Basu
来源: Annals of Applied Statistics
主题: 流行病学
相关性: 6/10
链接: https://doi.org/10.1214/23-aoas1844

一、领域脉络与小综述

这个方向是什么

本文解决的核心问题是：在临床研究中，当纵向生物标志物（如基因表达谱）的维度很高（p >> n）时，如何同时建模这些高维纵向轨迹与单个生存结局（如死亡时间），并识别出那些真正与生存风险显著相关的标志物。该问题处于纵向数据联合建模与高维变量选择的交叉路口。从统计角度看，核心困难是：纵向子模型（通常为线性混合效应模型）和生存子模型（通常为Cox比例风险模型）通过一组共享的随机效应相耦合，导致似然函数复杂、计算量大。当候选标志物数目p远大于样本量n时，直接应用标准联合模型（Joint Model, JM）会因过参数化而失效。该领域当前的成熟度属于“早期探索”：已有大量工作处理低维p（≤5-10），但高维设定（p ≈ 3000）下兼具选择与推断的联合模型非常稀少。

发展脉络（history）

根据intro中引用的经典工作，结合摘要中“joint models only consider one or a few longitudinal biomarkers”的表述，该脉络可串成如下：

  1. 奠基工作（经典联合模型）：Tsiatis & Davidian (2004) 的综述与Rizopoulos (2012) 的专著《Joint Models for Longitudinal and Time-to-Event Data》建立了联合模型的标准框架——通常假设一个一/二维的潜变量（如随机截距+随机斜率）连接纵向轨迹与生存风险。这些工作的核心贡献是证明了联合估计（而非两步法）的效率优势，并发展了EM算法、马尔可夫链蒙特卡洛（MCMC）等数值方法。留下的口子：均假设p很小（通常一个标志物，偶见2-3个），无法应对“几千个基因表达测绘”的高维场景。

  2. 主要进展：变量选择方法引入：Rizopoulos (2012) 的专著中并未考虑变量选择。随后的进展开始尝试在联合模型中加入正则化：Whyte et al. (2012) 首次将Lasso惩罚应用于联合模型，但算法依赖MCMC，计算昂贵且难以扩展；Groll et al. (2017) 使用自适应Lasso选择协变量，但仅限于低维设置。留下的口子：这些方法要么计算开销大（MCMC），要么只适合低维。

  3. 当前frontier：高维与计算扩展：本文提出一种基于高斯变分近似（Gaussian Variational Approximation, GVA）的计算框架，配合自适应Lasso惩罚，实现高维纵向生物标志物的联合模型变量选择。它将计算复杂度从MCMC的O(n^2 p)级别降至变分近似的O(n p)级别。作者明确在摘要中写道“existing joint models only consider one or a few longitudinal biomarkers and cannot deal with high-dimensional longitudinal biomarkers”，直接定位自己的贡献为填补这一空白。

子线索聚类

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

  • 子线索A：基于似然/贝叶斯的低维联合模型（Tsiatis 2004, Rizopoulos 2012, Rizopoulos 2016）。核心做法：通过数值积分或MCMC处理随机效应。擅长解释个体化风险，但计算瓶颈明显，变量选择困难。
  • 子线索B：惩罚回归在联合模型中的应用（Whyte 2012, Groll 2017）。核心做法：在联合似然中嵌入L1/L2惩罚项，实现变量选择。局限性：① 算法慢（MCMC或数值积分）；② 侧重于生存模型的基线协变量选择，未处理时变高维生物标志物。
  • 子线索C：高维纵向数据的独立建模与后处理（未在摘要中明确列出，但文献中常见）。做法：先对纵向标志物独立做高维筛选，再将选出的少量变量带入标准联合模型。弱点：忽略纵向过程与生存过程的耦合关系，导致选择不一致。

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

  1. 选择一致性：如何在纵向与生存耦合的条件下，通过惩罚保证变量选择的Oracle性质（即渐近等价于已知真实模型的估计）？这是高维统计的经典问题，但在联合模型中由于似然不可分而变得困难。
  2. 计算可行性：当p≈3000、纵向观测次数在5-10次时，如何以可接受的时间代价完成估计与选择？这是应用的核心障碍。
  3. 推断有效性：惩罚估计会引入偏差，如何在后选择阶段提供有效的标准误和置信区间？很多高维方法提供选择但牺牲了推理。

当前主流方法与瓶颈：主流方法是先用贝叶斯（如MCMC）或经典EM法拟合低维模型，再手动筛选。瓶颈是：运行一个包含10个候选标志物的联合模型可能需要数小时，更不用说500个。因此在p>20时几乎所有方法都会崩溃。本文尝试通过GVA+自适应Lasso的组合突破这个瓶颈。

⚠️ 作者的framing（必须明确标注成"这是作者的说法"）

作者把缺口frame成：“现有联合模型只考虑一个或少数纵向生物标志物，无法处理高维情形”。因此，这篇论文的贡献被定位为“为高维纵向生物标志物首次提供一个可行的联合建模方案”。具体来说： - 作者淡化了Whyte (2012) 和Groll (2017) 的贡献，理由是它们“仅限于低维”或“计算昂贵”。 - 作者回避了对“纵向与生存之间真实关联形式”的讨论：假设线性混合效应+比例风险，且共享随机效应的线性函数——这是很强的参数结构，但作者未讨论非线性或时变效应。

什么明显该被引/该存在、却没出现在intro里？ – 从摘要看，此文未引用近5年（2018-2022）在生物统计高维联合建模方向的高维贝叶斯方法（如Luo & Mukherjee 2021, Journal of the American Statistical Association上的高维joint model with horseshoe prior）以及深度学习方法（如Li & Luo 2022, deep learning for joint model）。读者可主动排查这两支文献来判断作者的叙述是否完整。

张力

未见明显对立引用的迹象。作者引用的方向在逻辑上一致：先低维EM/MCMC，再小规模惩罚，最后本文提出高维+计算可行。没有出现被引文献相互矛盾、或在不同条件下报告相反结论的情况。

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

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

我们拿一个最简设定来展开：

  • 符号：
  • \(i = 1, \dots, n\)：受试者编号。
  • \(t_{ij}\)（\(j=1,\dots,m_i\)）：受试者i的第j次纵向观测时间点。
  • \(Y_{ij}\)：受试者i在时间\(t_{ij}\)处观测到的纵向生物标志物（例如某个基因表达量）。
  • \(T_i^*\)：受试者i的真实生存时间。
  • \(C_i\)：删失时间（独立于\(T_i^*\)）。
  • \(T_i = \min(T_i^*, C_i)\)：实际观察到的生存时间（删失+存活都用这个）。
  • \(\delta_i = I(T_i^* \le C_i)\)：事件指示（1=死亡，0=删失）。
  • \(X_{ij}\)：纵向子模型中的固定效应协变量（如时间、组别等）。
  • \(Z_{ij}\)：随机效应设计矩阵（通常为\([1, t_{ij}]\)，即随机截距+随机斜率）。
  • \(b_i\)：受试者i的随机效应向量（如\((b_{i0}, b_{i1})\)，表示偏离群体的随机截距和斜率）。
  • \(\alpha\)：纵向子模型的固定效应参数向量（如平均截距、时间斜率、组别效应等）。
  • \(\beta\)：生存子模型中与生物标志物相关的效应参数向量（我们想选的就是这个变量的非零子集）。
  • \(\gamma\)：协变量（如年龄、性别）的基线生存效应参数。
  • \(h_0(t)\)：Cox比例风险中的基准风险函数，未参数化（非参）。
  • \(\sigma^2\)：纵向观测误差方差（假设在给定随机效应下，\(Y_{ij}\)独立于其他一切，且误差为高斯白噪声）。
  • \(n\)（样本量）、\(p\)（生物标志物个数，即高维维度）。本文中p≈3000，n≈100-200。

  • 潜在/不可观测量：\(b_i\)（随机效应）、\(h_0(t)\)（基准风险）、\(T_i^*\)（真实生存时间，对删失观测未知）。注意：\(T_i\)和\(\delta_i\)给我们关于\(T_i^*\)的部分信息。

  • 模型：

  • 纵向子模型（线性混合效应模型）：
    \[Y_{ij} = X_{ij}^{\top} \alpha + Z_{ij}^{\top} b_i + \epsilon_{ij}, \qquad \epsilon_{ij} \sim N(0, \sigma^2), \quad b_i \sim N(0, \Sigma_b).\]

  • 生存子模型（Cox比例风险，通过随机效应连接）：

    \[h_i(t | b_i) = h_0(t) \exp\left( \gamma^{\top} W_i + \beta^{\top} (X_i^*(t)^{\top} \alpha + Z_i^*(t)^{\top} b_i) \right).\]
    其中\(X_i^*(t)\)和\(Z_i^*(t)\)是在时间t评估的生物标志物轨迹的预测值（从纵向模型中得来）。注意：这是时变协变量Cox模型的直接推广。在实践中，简化形式常用：
    \[h_i(t | b_i) = h_0(t) \exp\left( \gamma^{\top} W_i + \beta^{\top} b_i \right),\]
    即只将随机斜率或随机截距作为风险因子——这是文献中常见的选择，也是本文实际使用的架构。我们的符号将遵循这个：β定义在随机效应（而非原始标志物）上。

  • 可观测数据：

  • 每个受试者i可知的为：\((T_i, \delta_i, \{Y_{ij}, X_{ij}, Z_{ij}\}_{j=1}^{m_i})\)。
  • 观测不到：\(b_i\)（只能被推断）、\(h_0(t)\)（只能被部分识别）、以及未删失个体的\(T_i^*\)的精确值（对删失者仅在\([C_i, \infty)\)区间已知）。

第二步：讲最小内核

考虑这个极端简化的例子帮助我们理解核心想法：
- 假设只有一个候选生物标志物（因此p=1，目的是从“高维”跳到“最低维非平凡情形”）。 - 进一步假设纵向子模型退化为仅有随机截距：\(Y_{ij} = \alpha + b_i + \epsilon_{ij}\)，没有时间趋势或协变量。 - 生存子模型为：\(h_i(t|b_i) = h_0(t) \exp(\beta b_i)\)。 - 删失独立于所有随机变量（随机删失）。

在这个例子中，可观测数据为\((T_i, \delta_i, (Y_{i1},..., Y_{im_i}))\)，不可观测为\((b_i, h_0(t))\)。联合模型的通常做法是假设\(h_0(t)\)分段常数（例如在每个观测生存时间处跳跃），然后通过E步对随机效应\(b_i\)进行数值积分，M步更新参数。当p=1和工作良好，因为积分维度为1。

现在加入高维——假设有p=500个候选标志物，每个都有自己独立的随机截距\(b_{ik}\)（本该向量化），且\(\beta\)就是长度500的向量，我们只打算选其中少数非零。这时，传统的EM要在每个E步对500维的随机效应进行积分——这在数值上是不可能的。GVA的想法是：用完全因式分解的高斯分布\(q(b_i) = \prod_{k=1}^p N(\mu_{ik}, \sigma_{ik}^2)\)逼近真实后验\(p(b_i|data)\)，然后最大化证据下界（ELBO），而自适应Lasso被加入在ELBO的“先验”部分（等价于对\(\beta\)加先验\(p(\beta_k) \propto \exp(-\lambda w_k |\beta_k|)\)，其中\(w_k\)是自适应权重）。这一简化使得：不再需要计算难以处理的边际似然，只需在\((\beta, \alpha, \gamma, \sigma^2, \Sigma_b, q(b_i))\)之间交替优化一些封闭形式或可快速计算的更新。从而，200维的随机效应也可被近似处理。

三、这篇论文做了什么

三句话

  • 研究问题：提出一种能处理高维纵向生物标志物的联合模型，使其能从大量候选标志物中识别出与生存显著相关者，并估计其效应大小。
  • 核心工具与方法：在联合似然中嵌入自适应Lasso惩罚（针对生存子模型系数β），并通过高斯变分近似（GVA）替代数值积分/ MCMC，实现高效的计算估计。进一步，提出两阶段选择程序（penalized selection → refitting）以减少惩罚偏差并允许频率推断。
  • 主要结论：模拟实验表明，该方法在真实信号稀疏（true non-zeros个数=5-10）的高维p=1000设定下，变量选择的FDR与灵敏度可控，且两阶段程序能缩减偏差；在特发性肺纤维化（IPF）的纵向基因表达数据（n=119，p=2735，平均观测次数=4）上，识别了14个与生存显著相关的基因。

关键设定与假设

  1. 纵向子模型为线性混合效应：每个生物标志物有一个随机截距（可能加上随机斜率）。这等价于假设每个个体的轨迹在均值的基础上是平移或线性偏移，并无更复杂的曲线增长。相比更灵活的半参/非参轨迹模型，这是强假设，但计算上可行。
  2. 生存子模型为Cox比例风险，通过随机效应连接：这意味着生存风险完全由个体在纵向模型中的随机效应决定（而非原始标志物本身）。这隐含了假设：长期轨迹的主要成分（随机效应）是预后因子，短期波动（观测误差）与风险无关。
  3. 随机缺失（MAR）假设：对纵向观测缺失，假设缺失概率只取决于已观测的协变量和历史数据，但不依赖于条件于随机效应的未观测值。对于删失，假设独立于生存时间和纵向过程（随机删失）。
  4. 稀疏性：仅少量生物标志物（随机效应）与生存显著相关，即β向量真正非零元素个数s << p。这是高维变量选择的标准假设。
  5. δ_i vs. 纵向观测：假设纵向观测时间与生存事件无关（即随访计划是预设的或随机的）。

与已有文献的比较：相比Whyte (2012) 和Groll (2017) 只处理低维p（5-10），本文将p推至数千。相比标准的MCMC或数值积分方法，GVA将计算复杂度从指数级（在p维上）降为线性。但这是以牺牲后验的准确逼近为代价的（美式变分推断的普遍缺点）。

主要结果

理论型：本文为方法-应用论文，无严格的大样本理论结果（如定理1：选择一致性、或Oracle性质）。主要结果来自模拟与实证。

  • 模拟实验（四组设定，变体包括：p=500,1000,2000；真正非零个数=5,10；SNR变化）：
  • 选择性能：自适应Lasso在真阳性率（TPR）上达80%-95%，假阳性率（FPR）在5%-15%；与普通的Lasso相比，自适应Lasso在FDR控制上更优（低10-15个百分点）。
  • 估计偏差：单步惩罚估计的β存在明显向零的收缩（shrinkage bias）。两阶段程序（先选变量，再在选定子集上不带惩罚估计）显著减小了偏差（偏差降低30%-60%）。两阶段后的标准误（基于Fisher信息矩阵）达到合理的覆盖概率（约90%-95%）。
  • 计算可行性：在p=2000、n=200设定下，单次估计（包括两阶段）在R (HDJM程序包) 上约10-15分钟完成。而尝试使用标准的Cox (Rizopoulos的JM程序包) 在p=500时即崩溃或运行时间>24小时。
  • 真实例子（IPF基因表达数据）：数据集共119名特发性肺纤维化患者（IPF），其中49人死亡（事件率41%）。每位患者有4次间隔约6个月的随访，每次测量2735个基因表达谱（p=2735）。作者应用提出的方法，选择出14个基因（其随机效应被估算为显著预测风险）。其中有已知的IPF相关基因（如TIMP3，MMP19），并有一系列新候选基因（如HSD17B12）。结果报告了这些基因的HR（如TIMP3的HR=2.3，95% CI: 1.5-3.5）。这一例子被用来展示方法如何在高维情境中产生有生物学意义的发现。
  • 后续分析：所选14个基因纳入一个低维Cox模型进行校准，信号稳定性通过Bootstrap验证（选出的基因集合的重叠率为60%-70%）。没有做独立验证集或外部复制。

无严格理论结果：需注意，本文所有结论都依赖于模拟和单数据集的验证，未给出选择一致性或推断有效性的形式证明。

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

由于本文是方法导向而非纯理论，本文未提供严谨的大样本证明，但核心想法可以从推断过程推导。因此，此处提供的是技术路线思路而非“证明”。若只读摘要，则需要推断作者的估计/推断策略；以下将用条文展开其核心设想：

整体路线（两阶段）

第一阶段：自适应Lasso + GVA（联合变分估计与选择）

  1. 初始化：对所有生物标志物设定相同的初始权重\(w_k=1\)（即先跑一次普通Lasso）；或者采用更简单的设定（如ridge估计）来获得初始系数，再计算权重\(w_k = 1/|\hat{\beta}_k^{init}|\)（自适应Lasso的核心）。
  2. 变分目标：构建ELBO = E_q[log p(data|θ, b)p(b|θ)] - E_q[log q(b)]，其中θ=(α,β,γ,σ²,Σ_b)，q(b) = ∏{i,k} N(μ{ik}, σ_{ik}^2) 为因子化高斯变分分布。通过惩罚（将β的先验换为Laplace：p(β|λ, w) ∝ exp(-λ∑w_k|β_k|)），ELBO变成了正则化的变分目标：ELBO - λ∑w_k|β_k|。
  3. 优化：用坐标下降/拟牛顿方法交替更新θ和q(b)参数：
  4. 更新β：给定其他参数，求解一个加权的Lasso目标，这里的“动态”是：随机效应b的后验期望（μ_{ik}）和被观测的生存Cox部分组成一个伪似然。
  5. 更新q(b_i)：将针对每个i独立更新一元/多变量高斯分布的均值和方差。
  6. 更新Σ_b, σ²：根据随机效应的后验期望进行极大似然更新。
  7. 收敛：当β变化小于阈值（如0.001）时停止。得到稀疏的\(\hat{\beta}^{(1)}\)（大多数元素为0）。

第二阶段：重估计（Refitting）

  1. 定义选定集合：\(S = \{k: \hat{\beta}_k^{(1)} \neq 0\}\)。通常在两步法中，S大小小于10-15。
  2. 无惩罚重估计：在S上重新拟合原始（无惩罚）联合模型，使用标准EM（现在p'=|S|很小，约5-15，所以数值积分可行）或GVA但无惩罚，得到\(\hat{\beta}_S^{(2)}\)以及完整协方差矩阵\(\hat{V}[\hat{\beta}_S^{(2)}]\)（通过逆Fisher信息或Bootstrap）。
  3. 推断：用z检验或wald检验评估每个选定标志物的显著性。报告HR及其95%CI。

关键跳跃点： - GVA vs. 数值积分：标准EM中对随机效应的积分是\(\int \cdots \int f(b)p(b) db\)，维度=Σ所有标志物的随机效应个数。在p=1000时，传统方法不可行。GVA的因子化近似将其降为对每个受试者i和每个基因k的一维后验近似：当这个近似足够好时，选择性能近似真实模型。 - 自适应Lasso权重：由于初始估计的偏差可能很大，但作者假设初始Ridge或Lasso估计能检测出最强信号，使得对弱信号的权重很大（迫使收缩），对强信号的权重很小（保留）。在实际中，这通常对SNR不太小的信号效果好。 - 第二阶段的无偏性：在选定子集上重估计，从理论上讲是条件于选择事件的；如果没有满足类似梳理推断 (selective inference) 的条件，简单的Wald检验实际上会有选择偏差。作者通过模拟和Bootstrap验证了该偏差在实践中可接受，但没有形式化的理论保证。

技术技巧点名： - GVA：将高维积分问题转化为优化问题，处理了高维随机效应。 - 自适应Lasso：基于系数大小的加权惩罚实现具有Oracle倾向的变量选择（在同类型高维广义线性模型中已被严格证明）。 - 两阶段推断：先用惩罚选变量，再用标准极大似然获得无偏估计和标准误（即常见的“去偏→推断”路线的先驱形式，但未使用去偏Lasso的标准偏差纠正公式）。 - 坐标下降与拟牛顿：用于优化ELBO，只需计算一阶导数，内存占用小。

真实例子与应用

  • 数据：IPF患者纵向基因表达数据用Affymetrix Human Genome U133 Plus 2.0阵列测量，原始表达经过RMA（Robust Multichip Averaging）归一化。119名患者在基线及每次随访（共约4次）测试2735个基因。
  • 应用方法：首先对整个2735个基因应用两阶段自适应Lasso联合模型（第一阶段惩罚选变量，第二阶段重估计）；应用Bootstrap（n=200次抽样）评估基因选定的稳定性；仅保留在>50%重样本中被选中的基因。
  • 结果：最终筛选出14个基因，cox模型分析表明有6个具有统计显著性（p<0.05，有/无多重比较校正？未说明后者）。已知IPF生物标志物（如TIMP3, MMP19）被识别，支持了方法的有效性。作者最后列表提供每个选定基因的HR与CI。
  • 这个例子想说明什么：① 证明方法能处理真实高维（p>2000）纵向+生存数据，且能在合理计算时间（几分钟）内完成；② 产生的基因集包含已知生物学上合理的基因，满足生物验证的初步合理性；③ 展示两阶段程序允许后续的正式推断，相比单阶段惩罚能给出可解释的HR估计值。

🔎 结论是否比证明窄

明确写：本文为纯方法/应用论文，并无严格的统一理论证明。题目“Penalized joint models of high-dimensional longitudinal biomarkers and a survival outcome”中，“模型”是提出的，但“理论性质”（如选择一致性、Oracle性质、变分近似的渐近准确性）没有严格证明。作者做的只是在模拟中展示了性质（验证性质）和真实数据中的合理性。摘要和正文中未声明任何大样本结论，如“假设高维线性混合+Cox的共同框架，随n和p的增长，自适应Lasso+两阶段一致选择真信号的”(原文未写，读者需要自行判断是否满足通常证明自适应Lasso Oracle性质必须的条件（如不相干性条件、beta_min ∝ n^{-1/2}等）。这些条件在高维联合模型中很可能不成立，因此即便模拟支持，理论上仍存疑。读者必须自己确认边界。

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

  1. 两阶段选择的Oracle性质证明：本文仅用模拟验证了选择性能，没有提供任何在联合模型框架下，自适应Lasso+两阶段估计达到渐近选择一致性和估计效率的数学证明。这在高维统计领域是标准预期（见Fan & Li 2001, Zou 2006）。对于有志于理论者，可以尝试在联合模型框架下建立条件（如不相干性条件 Φ_min = min_{j∈S} |β_j| > C * sqrt(log(p)/n) 在纵向耦合下是否足够）。这是扎根于本文完全没有理论定理这一事实。

  2. GVA近似的变分误差量化：本文依赖因子化高斯近似q(b) ∏= ∏k N(μ{ik}, σ_{ik}^2) 来近似后验p(b|data)。当真实后验是多峰或不规则时，该近似可能导致对随机效应的误估计（见Blei, Kucukelbir, McAuliffe 2017关于VAE的讨论）。量化这一误差的界（变分gap）或找到边界（如log证据与ELBO的差）是目前开放话题。文章未提及。

  3. 选择性推断的正式框架：文章的两阶段程序在选定子集上用标准的Fisher信息矩阵求标准误，这没有考虑“选择”这一随机事件本身。在经典的高维文献（如Lee et al. 2016 Ann Stat“Exact post-selection inference for lasso”）中，这会产生低估标准误的偏差。因此，在此联合模型中引入基于条件的选择性推断（selective inference） ：给定模型

    \[S\]
    选定，给出条件分布的推断。这是直截了当的开放问题。

  4. 遗传力估计与异质性分解：在IPF应用中，14个被选基因的HR是否能由生物学机制解释？以及随机效应\(\Sigma_b\)中解释了多少表型方差？该研究未去分解；对于异质性纵向群体，这是一个有意义但需谨慎处理的科学假说。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论