跳转至

FastJM: An R Package for Efficient Implementation of Semiparametric Joint Models for Longitudinal and Survival Data

作者: Shanpeng Li, Emily Ouyang, Ace Isabel Mejia-Sanchez, Xinping Cui, Gang Li
主题: 统计计算 / 算法
相关性: 7/10
链接: https://arxiv.org/abs/2608.14127

一、领域脉络与小综述

这个方向是什么

这个子方向是纵向数据与生存数据的半参数联合建模，其根本的统计问题是：如何同时建模一个（或多个）随时间重复测量的生物标志物（纵向过程）和一个（或多个）时间至事件结局（生存过程），并刻画两者之间的潜在关联。该方向的核心挑战在于：纵向过程通常含有测量误差和不规则观测，而生存过程可能包含竞争风险（多个互斥的事件类型）和右删失。当前成熟度较高，已有多个R包实现，但计算可扩展性——特别是面对大规模电子健康档案（EHR）和生物银行数据时——仍是主要瓶颈。

发展脉络（history）

  1. 奠基工作：Wulfsohn and Tsiatis (1997) 提出了通过共享随机效应连接纵向和生存过程的联合建模框架，奠定了该方向的基础。Henderson et al. (2000) 和 Tsiatis and Davidian (2004) 进一步推广了这一框架，并系统阐述了其统计性质。Elashoff et al. (2008) 将其扩展到竞争风险设定。

  2. 主要进展与软件实现：早期R包如 JM (Rizopoulos, 2010b) 和 joineR (Hickey et al., 2018a) 实现了单个生物标志物的联合模型，支持竞争风险。JMbayes (Rizopoulos, 2016) 提供了贝叶斯MCMC实现。JSM (Xu et al., 2020) 引入了变换生存模型。这些包主要聚焦于单个纵向生物标志物。

  3. 多变量与复杂结构：近年来，软件支持扩展到多变量纵向过程，包括 joineRML (Hickey et al., 2018b, MCEM)、gmvjoint (Murray and Philipson, 2022, 2023, 近似EM)、JMbayes2 (Rizopoulos and Taylor, 2024, 贝叶斯MCMC)、rstanarm (Goodrich et al., 2024, Stan-based) 和 INLAjoint (Rustand et al., 2024, INLA)。FlexVarJM (Courcoul et al., 2025) 通过混合效应位置-尺度模型引入了异质性个体内（WS）变异性。

  4. 当前前沿与本文位置：当前前沿是计算可扩展性与模型灵活性的平衡。本文（FastJM）定位为：在频率主义EM框架下，通过定制的线性扫描算法解决半参数联合建模中非参数基线风险函数更新的计算瓶颈（从O(n²)降至O(n)），并支持三种模型类（单生物标志物、多生物标志物、异质性WS变异性）。此外，本文还引入了一个landmark多变量联合建模框架，以支持时变潜在关联结构，同时保留计算优势。

子线索聚类

  • 线索一：单生物标志物联合模型。代表：JM, joineR, JSM, JMbayes。核心是处理一个纵向结果与一个生存结局的关联，通常通过共享随机效应实现。FastJM的jmcs()属于此簇，但通过线性扫描算法提升了计算效率。

  • 线索二：多变量纵向联合模型。代表：joineRML, gmvjoint, JMbayes2, rstanarm, INLAjoint。核心是同时建模多个相关纵向生物标志物与生存结局。FastJM的mvjmcs()属于此簇，采用正态近似（Murray and Philipson, 2022）避免高维张量积求积。

  • 线索三：异质性WS变异性联合模型。代表：FlexVarJM。核心是允许纵向结果的残差方差随受试者和访视时间变化，并评估其与事件风险的关联。FastJM的JMMLSM()属于此簇，采用自适应Gauss-Hermite求积。

  • 线索四（本文新增）：Landmark多变量联合建模。这是本文提出的扩展，将landmark分析与多变量联合模型结合，支持时变潜在关联结构（如当前值、潜在过程当前值），同时保留线性扫描算法的计算优势。

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

  1. 如何高效计算非参数基线风险函数？ 这是半参数联合模型的主要计算瓶颈，因为需要反复评估风险集。现有方法（如直接实现）复杂度为O(n²)。FastJM通过线性扫描算法将其降至O(n)。

  2. 如何处理高维随机效应带来的数值积分困难？ 当生物标志物数量增加时，随机效应维度上升，张量积求积变得不可行。FastJM的mvjmcs()采用正态近似（Murray and Philipson, 2022）来避免高维积分。

  3. 如何实现时变潜在关联结构（如当前值关联）而不牺牲计算效率？ 传统的共享随机效应参数化假设时间不变关联。FastJM通过landmark多变量联合建模框架，在保留线性扫描算法优势的同时，支持时变关联。

  4. 如何评估动态预测性能？ 需要交叉验证的时变准确度指标（如AUC、Brier Score）和时不变一致性统计量（C-index）。FastJM提供了统一的接口。

⚠️ 作者的 framing

作者将缺口 frame 成：“计算可扩展性仍是拟合复杂联合模型的核心障碍”（原文："computational scalability remains a central obstacle to fitting complex joint models"）。他们通过线性扫描算法和模型特定的E步近似来解决这个问题，使FastJM成为“显然的下一步”——一个在频率主义框架下兼顾计算效率和模型灵活性的统一软件。

被淡化或回避的竞争路线： - 贝叶斯方法（如JMbayes2, rstanarm, INLAjoint）被提及但未深入比较。作者暗示频率主义EM框架在计算上更具可扩展性，但未讨论贝叶斯方法在不确定性量化或处理复杂随机效应结构方面的潜在优势。 - 完全参数化基线风险（如使用样条）未被讨论。FastJM坚持非参数基线风险，这虽然增加了计算负担（即使有线性扫描），但提供了更大的灵活性。

什么明显该被引/该存在、却没出现在intro里？ - 没有引用关于联合模型渐近理论的工作（如Zeng and Cai, 2005; Zeng et al., 2005），尽管这些工作为半参数MLE的推断提供了理论基础。FastJM的标准误估计正是基于这些理论（profile似然方法）。 - 没有引用关于计算-统计权衡或算法复杂度分析的文献。线性扫描算法的O(n)复杂度是一个计算贡献，但作者没有将其置于更广泛的统计计算理论背景下讨论。

张力

未见明显对立引用。不同包在计算方法（频率主义EM vs. 贝叶斯MCMC vs. INLA）和模型假设（共享随机效应 vs. 时变关联）上存在差异，但作者将其视为互补而非对立。FastJM的landmark扩展正是为了弥合共享随机效应和时变关联之间的差距。

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

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

符号： - i = 1, ..., n: 受试者索引。 - n: 样本量。 - K: 竞争风险事件类型数。 - eTi: 受试者i的真实事件时间。 - eDi ∈ {1, ..., K}: 受试者i的真实事件类型。 - Ci: 受试者i的独立、非信息性删失时间。 - Ti = min(eTi, Ci): 观测到的事件或删失时间。 - Di = eDi * I(eTi ≤ Ci): 观测到的事件类型（0表示删失）。 - Yi(t): 受试者i在时间t的纵向生物标志物值。 - tij, j = 1, ..., ni: 受试者i的第j次观测时间。 - ni: 受试者i的重复测量次数。 - Xi(t): 纵向子模型中固定效应的设计向量（可能含时间）。 - Zi(t): 纵向子模型中随机效应的设计向量。 - β: 纵向子模型的固定效应系数向量（参数）。 - bi: 受试者i的随机效应向量（潜在变量）。 - σ²: 纵向子模型的残差方差（参数）。 - Wi: 生存子模型中的基线协变量向量。 - γk: 对事件类型k的基线协变量回归系数向量（参数）。 - λ0k(t): 事件类型k的未指定基线风险函数（非参数）。 - αk: 连接纵向和生存过程的关联参数向量（参数）。 - Ψ: 所有模型参数和基线风险函数的集合。

模型（以最简单的单生物标志物模型为例）： - 纵向子模型（线性混合效应模型）： Yi(t) = Xi(t)ᵀβ + Zi(t)ᵀbi + εi(t), 其中 εi(t) ~ N(0, σ²), bi ~ N(0, Σ)。 - 生存子模型（原因别Cox比例风险模型）： λik(t | Wi, bi) = λ0k(t) * exp(Wiᵀγk + biᵀαk), 对于 k = 1, ..., K。 - 连接机制：两个子模型通过共享随机效应 bi 连接。αk 量化了随机效应对事件类型k的风险的影响。

可观测数据： - 纵向数据：{Yi(tij), tij, Xi(tij), Zi(tij)} 对于 i = 1, ..., n, j = 1, ..., ni。这是研究者实际能观测到的。 - 生存数据：{Ti, Di, Wi} 对于 i = 1, ..., n。这是研究者实际能观测到的。 - 潜在/不可观测量： - 随机效应 bi（只能通过假设和EM算法推断）。 - 基线风险函数 λ0k(t)（非参数，通过EM算法估计）。 - 真实事件时间 eTi 和类型 eDi（当被删失时）。 - 测量误差 εi(t)。

第二步：讲最小内核

最简特例：考虑一个受试者 i，只有一个纵向生物标志物，且只有两种竞争风险（K=2）。纵向子模型只包含随机截距（Zi(t) = 1, bi = bi0），没有随机斜率。生存子模型中没有基线协变量（Wi 为空）。那么模型退化为： - 纵向：Yi(t) = β0 + β1*t + bi0 + εi(t), εi(t) ~ N(0, σ²), bi0 ~ N(0, σ²_b)。 - 生存：λik(t | bi0) = λ0k(t) * exp(αk * bi0), k = 1, 2。

核心思路：在这个特例下，EM算法的E步只需要对一维随机效应 bi0 进行积分，可以用Gauss-Hermite求积。M步的关键是更新非参数基线风险 λ0k(t)。直接实现需要为每个不同的事件时间重新计算风险集，复杂度为O(n²)。线性扫描算法的核心想法是：将受试者按观测事件时间 Ti 排序，然后通过一次从后向前的扫描，递归地计算风险集上的累积和，从而将更新基线风险的复杂度从O(n²)降至O(n)。具体来说，对于事件类型k，在排序后，λ0k(t) 的Breslow估计器可以写成： Λ̂0k(t) = Σ_{l: t_kl ≤ t} [d_kl / Σ_{r ∈ R(t_kl)} exp(α̂_k * b̂_r0)]， 其中 t_kl 是第l个类型k的事件时间，d_kl 是该时间的事件数，R(t_kl) 是风险集。线性扫描通过预先计算并存储 Σ_{r ∈ R(t_kl)} exp(α̂_k * b̂_r0) 的累积和，避免了每次重新求和。

为什么成立：因为生存子模型中的协变量（这里是 b̂_r0）是时间独立的，所以风险集上的和可以写成从当前时间到最大时间的累积和。通过一次排序和一次扫描，所有需要的和都可以在O(n)时间内得到。这个特例抓住了论文计算贡献的本质：利用时间独立协变量和事件时间排序，将风险集计算从O(n²)降为O(n)。

三、这篇论文做了什么

三句话

  1. 研究了什么问题：本文介绍了R包 FastJM，用于高效实现纵向数据与生存数据的半参数联合模型，解决了大规模生物医学数据（如EHR和生物银行）中联合模型拟合的计算瓶颈。
  2. 核心工具/方法：在EM算法框架下，采用定制的线性扫描算法来高效更新非参数基线风险函数（从O(n²)降至O(n)），并结合模型特定的E步近似（伪自适应Gauss-Hermite求积、正态近似、自适应Gauss-Hermite求积）来处理不同复杂度的随机效应结构。
  3. 主要结论：FastJM 提供了三种模型类（单生物标志物、多生物标志物、异质性WS变异性）的统一频率主义工作流，包括模型拟合、推断、可视化、动态预测和预测性能评估。此外，还引入了一个landmark多变量联合建模框架，以支持时变潜在关联结构，同时保留计算优势。

关键设定与假设

在第二节最小记号的基础上，补全完整设定：

  • 共享随机效应参数化：所有三个核心模型（jmcs, mvjmcs, JMMLSM）都采用共享随机效应参数化，即纵向和生存过程通过受试者特定的随机效应 bi（或 θi）连接。这假设了时间不变的关联结构。
  • 条件独立性：给定协变量和随机效应，纵向过程和事件过程是独立的。这是联合模型的标准假设，使得似然可以分解。
  • 非信息性删失：删失时间 Ci 独立于事件时间和纵向过程，给定协变量和随机效应。
  • 纵向子模型：所有模型假设连续纵向生物标志物，使用高斯混合效应子模型。测量误差独立同分布于 N(0, σ²)（或 N(0, σ²_i(t)) 对于异质性WS模型）。
  • 生存子模型：使用原因别Cox比例风险模型，基线风险函数 λ0k(t) 完全未指定（非参数）。
  • 竞争风险：支持K个互斥的事件类型。
  • 相比已有文献的放宽/强化：
  • 放宽：mvjmcs() 通过正态近似避免了高维张量积求积，放宽了对低维随机效应的依赖。
  • 强化：JMMLSM() 引入了异质性WS变异性，这是对标准同方差假设的强化。
  • 强化：线性扫描算法将计算复杂度从O(n²)降至O(n)，这是对计算效率的显著强化。

主要结果

本文是应用/方法型论文，核心结果是软件实现和算法贡献，而非新的统计理论。主要结果包括：

  1. 线性扫描算法：将M步中更新非参数基线风险函数的计算复杂度从O(n²)降至O(n)。同样，标准误估计（profile似然方法）的复杂度也从O(n³)降至O(n)。这是本文最核心的计算贡献。
  2. 三种模型类的实现：
  3. jmcs(): 单生物标志物 + 竞争风险。E步使用伪自适应Gauss-Hermite求积（Rizopoulos, 2012）。
  4. mvjmcs(): 多生物标志物 + 竞争风险。E步使用正态近似（Murray and Philipson, 2022），避免高维积分。
  5. JMMLSM(): 单生物标志物 + 异质性WS变异性 + 竞争风险。E步使用自适应Gauss-Hermite求积。
  6. Landmark多变量联合建模框架：这是本文的方法论贡献。它结合了landmark分析和多变量联合模型，支持两种时变潜在关联结构：
  7. latAsso = "present": 当前值参数化，Aig(s) = mig(s)（潜在均值轨迹在landmark时间s的值）。
  8. latAsso = "presentlp": 潜在过程当前值参数化，Aig(s) = Zig(s)ᵀbig（随机效应贡献在landmark时间s的值）。 该框架保留了线性扫描算法的计算优势，因为关联函数在landmark时间s处固定，不随时间变化。
  9. 统一的工作流：提供了从数据模拟、模型拟合、诊断、动态预测到预测性能评估（交叉验证的AUC、C-index、Brier Score、MAEQ）的完整接口。

证明路线与技术技巧

本文是软件论文，没有传统意义上的定理证明。但EM算法的推导和线性扫描算法的正确性可以视为“证明”。

线性扫描算法的逻辑主干： 1. 问题：在M步中，更新非参数基线风险 Λ0k(t) 需要计算每个不同事件时间 t_kl 的风险集和 S(t_kl) = Σ_{r ∈ R(t_kl)} exp(η_r)，其中 η_r 是线性预测器。直接计算每个 S(t_kl) 需要遍历所有 r，总复杂度O(n²)。 2. 关键观察：如果协变量是时间独立的（如基线协变量和随机效应），那么风险集 R(t) 是单调递减的（随着t增加，风险集缩小）。因此，S(t) 可以写成从当前时间到最大时间的累积和。 3. 算法： - 将所有受试者按观测事件时间 Ti 升序排序。 - 初始化一个累积和变量 cumsum = 0。 - 从最大事件时间开始，向前扫描（或从最小开始，向后扫描）： - 对于每个事件时间 t_kl，S(t_kl) 等于当前累积和加上所有事件时间等于 t_kl 的受试者的 exp(η_r)。 - 更新累积和。 - 这样，每个受试者只被访问一次，总复杂度O(n)。 4. 正确性：这本质上是Breslow估计器的标准计算技巧，利用了风险集的嵌套结构。FastJM 将其扩展到竞争风险设定，并为每个事件类型分别维护累积和。

技术技巧点名： - 线性扫描算法：用于M步中基线风险更新和标准误估计。核心是排序和累积和。 - 伪自适应Gauss-Hermite求积（jmcs）：在EM迭代前，使用初始线性混合模型拟合的后验众数和曲率来重新中心和缩放求积节点，减少所需节点数。 - 正态近似（mvjmcs）：在每个EM迭代中，用后验众数处的多元正态分布近似随机效应的后验分布，从而解析计算所需条件矩，避免高维数值积分。 - 自适应Gauss-Hermite求积（JMMLSM）：在每个E-step，根据当前后验众数和曲率重新中心和缩放求积节点，以处理非高斯后验分布。 - 矩生成函数（MGF）近似：在landmark模型的M-step中，用于近似 E[exp(Ai(s)ᵀαk)] 等期望，其中 Ai(s) 是随机效应的线性函数。利用正态近似下MGF的解析形式。 - Profile似然标准误：通过逆经验Fisher信息矩阵估计参数标准误，其中得分向量通过EM算法收敛后的期望完整数据对数似然的导数近似。

真实例子与应用

本文使用模拟数据来演示所有功能，没有使用真实数据。这是论文的一个明确局限。

  • 数据：使用包自带的模拟函数（simJMdata, simmvJMdata, simJMWSVdata）生成数据。
  • 场景：
  • 单生物标志物（jmcs）：n=1000，两个竞争风险，纵向有随机截距和斜率。
  • 多生物标志物（mvjmcs）：n=5000，三个生物标志物，每个有不同的随机效应结构（截距、截距+斜率、截距+斜率+二次项），两个竞争风险。
  • 异质性WS变异性（JMMLSM）：n=1000，两个竞争风险，纵向有随机截距和斜率，方差子模型包含固定效应和随机截距。
  • Landmark多变量（mvjmcs with landmark）：n=5000，三个生物标志物，landmark时间s=5。
  • 结果：展示了参数估计、标准误、诊断图、动态预测（累积发生率函数）和预测性能评估（AUC、C-index、Brier Score、MAEQ）。例如，在多生物标志物例子中，模型在5000个受试者、30399个观测上，使用并行计算，运行时间为7.38分钟。
  • 例子想说明什么：
  • 验证了软件的正确性和可用性（参数估计接近真实值）。
  • 展示了线性扫描算法的计算效率（多生物标志物例子在合理时间内完成）。
  • 演示了异质性WS变异性模型相对于同方差模型的改进（通过诊断图显示残差方差异质性）。
  • 展示了landmark框架的灵活性（支持时变关联结构）。

🔎 结论是否比证明窄

  • 线性扫描算法的O(n)复杂度：论文声称将复杂度从O(n²)降至O(n)。这个结论在时间独立协变量的假设下是严格成立的。如果生存子模型包含时间依赖协变量（如当前值关联），则风险集和不能简单地通过一次扫描计算，线性扫描算法不再适用。论文的landmark扩展正是为了在保留线性扫描优势的同时引入时变关联，但这是通过将关联固定在landmark时间来实现的，并非真正的时变协变量。
  • 正态近似的准确性：论文引用Murray and Philipson (2022) 的工作，声称正态近似在参数估计上表现合理，即使纵向随访次数较少。但这是一个经验观察，并非严格证明。对于某些极端情况（如随机效应后验分布严重非高斯），近似可能失效。
  • 标准误估计：论文使用profile似然方法，并通过线性扫描算法加速。但profile似然标准误的渐近性质（如一致性）依赖于半参数MLE的正则性条件，这些条件在论文中没有被显式验证或讨论。结论“计算可行”是成立的，但“统计上有效”依赖于未证明的假设。

四、开放问题

  1. 将线性扫描算法扩展到时间依赖协变量：当前线性扫描算法严格依赖于生存子模型中协变量的时间独立性。如何将其扩展到支持真正的时变协变量（如当前值关联），同时保持O(n)复杂度？这是一个开放的计算问题。扎根于：论文第5节承认共享随机效应参数化“does not directly accommodate several commonly used time-dependent latent association structures”，而landmark扩展只是部分解决方案。

  2. 将正态近似扩展到非高斯纵向结果：FastJM 目前只支持连续高斯纵向结果。如何将正态近似（或类似的高效近似方法）扩展到二元、序数或计数纵向结果，同时保持计算可扩展性？扎根于：论文第6节“Future development”提到将扩展到“non-Gaussian longitudinal responses”。

  3. 将landmark框架扩展到更一般的关联结构：当前landmark框架只支持两种关联结构（当前值和潜在过程当前值）。如何支持更复杂的关联，如斜率、面积或累积暴露？扎根于：论文第5节提到“accommodates different latent association structures through the specification of the association function Aig(s)”，但只实现了两种。

  4. 计算复杂度的更精细分析：论文声称线性扫描算法将复杂度从O(n²)降至O(n)。但这是针对最坏情况的分析。对于实际数据，常数因子和内存访问模式可能很重要。能否从tensor-contraction/einsum复杂度的角度，对联合模型EM算法的计算图进行更精细的分析，以识别进一步的优化机会？扎根于：论文第2.3.2节描述了线性扫描算法，但没有进行常数因子或内存复杂度的分析。这与研究者的武器库（高阶U统计量的树宽/张量收缩/einsum计算）直接相关。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论