跳转至

Multistate approach for stochastic interventions on a time-to-event mediator in the presence of competing risks: A new R command within the CMAverse R package

作者: Ziqing Wang, Baoyi Shi, Cécile Proust-Lima, Hélène Jacqmin-Gadda, Linda Valeri
来源: Epidemiology
主题: 因果推断
相关性: 8/10
机构绿灯: Columbia University（US News 前 50，免分进入精读）
链接: https://doi.org/10.1097/ede.0000000000001791

一、领域脉络与小综述

这个方向是什么

本文的核心方向是因果中介分析（causal mediation analysis），具体聚焦于中介变量和结局变量均为时间-事件（time-to-event）且存在竞争风险（competing risks） 的设定。该子方向试图回答：在一个暴露（A）、一个时间-事件中介（M）和一个时间-事件结局（S）的因果路径中，如何定义和估计中介效应，尤其当暴露（如种族、社会经济地位）是“不可干预的社会建构”时，如何通过随机干预（stochastic intervention） 来量化“如果改变中介分布”会带来的结局差异。该领域的成熟度较低，方法虽在理论上有进展，但缺乏现成的、可复现的软件工具，这正是本文试图填补的缺口。

发展脉络（history）

  • 奠基工作：VanderWeele & Robinson (2014; 参考文献[7]) 系统讨论了种族等不可干预变量的因果解释，为“随机干预替代实际干预”提供了理论依据。Jackson & VanderWeele (2018; [6]) 进一步将分解分析（decomposition analysis） 引入健康差异研究，提出通过改变中介分布来量化差异来源。
  • 主要进展：Valeri 等人 (2023; [5]) 在本文引用的主要工作——将多状态模型（multistate modeling） 引入时间事件中介分析，在竞争风险存在的前提下，定义了 RD (残差异) 和 SD (转移分布效应) 两个因果估计量，并给出了其基于转移特定危险（transition-specific hazard）的表达。这是本文方法的理论前身。Shi 等人 (2021; [1]) 开发了 CMAverse R 包，提供统一的、可复现的因果中介分析函数集，但当时只支持连续/分类中介与连续/分类结局，不支持时间事件中介与竞争风险。
  • 当前frontier：学界开始关注更复杂的中介变量类型（如生存时间）和更现实的结局（如竞争风险）。已有多状态建模的软件工具（mstate 包，de Wreede et al., 2011; [8]）可供调用，但缺乏将因果参数（RD/SD）定义、估计与推断整合为一套现成命令的软件实现，导致方法虽在，但应用门槛很高。
  • 本文的位置：本文通过将 CMAverse 包与新实现的 cmest_multistate 命令结合，“封装”了 Valeri et al. (2023) 的理论方法，使其可直接由应用研究者使用。它本身不是方法论突破，而是一个工具性贡献：将一个已发表的、有效的统计方法转化为可复现的软件功能。

子线索聚类

被引文献大致落在两条线索上： 1. 理论方法线（Valeri 2023, VanderWeele 2014, Jackson 2018）：侧重于定义因果效应、提出识别假设、推导参数表达。这是“要做什么”和“理论上怎么做”的层面。 2. 软件实现线（Shi 2021, de Wreede 2011）：侧重于将方法封装为可执行的代码，提供可复现的管道。这是“能做多少”和“谁会用”的层面。CMAverse [1] 和 mstate [8] 是两个关键的底层依赖。

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

  • 核心问题 1：当中介和结局都是时间事件且存在竞争风险时，如何定义有意义的因果效应？
  • 答：已有理论提供了 RD 和 SD，通过随机干预来定义。
  • 核心问题 2：这些效应的识别假设是什么？如何验证？
  • 答：核心是“无未测量的中介-结局混杂”以及无暴露-中介/暴露-结局混杂。这些假设与标准中介分析相似但更严格。但本文并未展开讨论如何测试这些假设。
  • 核心问题 3：如何高效、稳健地估计和推断这些效应？
  • 答：使用多状态模型的转移特定半参 Cox 模型 + 非参 bootstrap。当前瓶颈在于：当时间点很多时，积分计算可能耗时；以及 bootstrap 的计算成本。

⚠️ 作者的 framing

  • 作者的说法：作者把缺口 frame 成 “现有R包（CMAverse）无法处理时间事件中介 + 竞争风险”。他们的论文是“扩展 CMAverse”以填补这一空白，从而让健康差异研究者能直接使用这一方法。
  • 被淡化/回避的竞争路线：作者回避了讨论其他可能的建模选择（如基于Cox的逆概率加权、G-formula、或更灵活的加速失效时间模型）。他们直接选择半参 Cox + 多状态模型，这是一个特定但流行的技术选择。他们没有与其他更灵活的建模方式（如parametric Weibull 或 flexible additive hazards）做仿真对比，也没有提供理论戒律说明为何这个选择是“最佳”。
  • 什么明显该被引/存在但没出现在intro里？：mstate 包（[8]）被引了，但文献中大量关于时间事件因果中介分析的识别框架（如 VanderWeele 2011 在统计年鉴上的综述） 并未出现在参考文献中。另外，关于非参数bootstrap在时间事件分析中的置信区间覆盖率的理论讨论（如超高效但覆盖不足的 bootstrap 对 Cox 模型偏倚效应的影响）也未提及。这是一个值得研究者去查的问题——本文是否过于依赖bootstrap的简单应用，而忽略了更先进的推断方法（如稳健的 influence function-based 方差估计）？

张力

未见明显对立引用。主要被引工作之间（Valeri 2023, VanderWeele 2014, Jackson 2018）在理论视角上是一致的。唯一可能存在的张力是：mstate 包基于的统计理论（事件型数据）与因果推断的识别假设（如无未测量混杂）之间，存在一个强假设：即模型必须正确指定。本文和 Valeri 2023 承认这一点，但并未讨论模型错误指定的后果。这算不上“矛盾”，但可以被视为一个“保守的立场”。

二、最核心、最简单的例子 / 数学问题（先把符号 / 模型 / 可观测数据交代清楚）

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

  • 符号：
  • A：暴露变量（二值，0/1，如“社会处境不佳” vs “社会处境佳”）。这是研究者无法干预的“社会建构”。
  • M：时间事件中介变量（如“发生某种疾病事件的时间”）——是一个潜在时间。
  • S：时间事件结局变量（如“死亡时间”）——也是一个潜在时间。
  • C：基线协变量向量（如年龄、性别、收入）。
  • gₐ：随机干预分布，代表从暴露组A=a的中介时间分布中随机抽取一个时间来替代个体的实际中介时间。这里 g₀（暴露组=0的分布）和 g₁（暴露组=1的分布）是核心概念。
  • RDₛ (残差异)：P(S^{g₀} > s | A=1, C) - P(S^{g₀} > s | A=0, C)。在给定协变量C的条件下，暴露组(A=1)和未暴露组(A=0)在超过时间s的生存概率之差，假设中介时间M的分布被固定为“未暴露组”(g₀)的分布。
  • SDₛ (转移分布效应)：P(S^{g₁} > s | A=1, C) - P(S^{g₀} > s | A=1, C)。对于暴露组(A=1)的人来说，如果他们的中介时间分布从现有的变为未暴露组(g₀)的分布，其生存概率的改变。
  • δ：事件指示变量（0=删失，1=事件发生）。
  • λ₀₁、λ₀₂、λ₁₂：三个转移的转移特定风险函数。λ₀₁：从基线状态→中介事件的风险；λ₀₂：从基线状态→结局事件的风险；λ₁₂：从中介事件→结局事件的风险。图B展示了这三个转移。
  • λ₋₋₀：基准风险（baseline hazard），对应 Cox 模型中的常数项。

  • t：时间索引，为估计量定义在不同时间点s上。

  • 模型：

  • 数据生成机制假设满足三个独立的半参数Cox比例风险模型，分别建模三个转移（λ₀₁、λ₀₂、λ₁₂）：
    • λ₀₁(t|A,C₁,C₂) = λ₀₁₀ * exp(a₁A + a₂C₁ + a₃C₂)
    • λ₀₂(t|A,C₁,C₂) = λ₀₂₀ * exp(b₁A + b₂C₁ + b₃C₂)
    • λ₁₂(t|A, M, C₁, C₂) = λ₁₂₀ * exp(c₁A + c₂M + c₃C₁ + c₄C₂)

  • 关键是：暴露A和两个协变量C₁、C₂是“外生”的；中介事件时间M和结局事件时间S是依赖上述模型独立生成的，但存在一个竞争关系：如果中介事件时间先于结局事件时间（t_M < t_S），则个体先经历中介事件，然后结局时间重新生成（包含中介M作为协变量）；如果结局事件时间先于或等于中介事件时间，则个体直接经历结局事件（中介事件被“删失”）。删失时间U是独立于所有事件的均匀分布随机变量。

  • 可观测数据：

  • 每个个体的数据集样包括（A, C₁, C₂, t_M, δ_M, t_S, δ_S）。
  • 可观测的：暴露状态A、基线协变量C₁、C₂；事件时间 t_M 和 t_S；事件/删失指示变量 δ_M 和 δ_S。
  • 潜在/不可观测：在无随机干预情况下的生存概率 P(S^{gₐ} > s)。这些只能依靠多状态模型和假设来识别。

第二步：讲最小内核（去掉一般性假设后的特例）

最简特例：无竞争风险 + 中介事件总发生

设想一个更极端的场景，允许简化 RD 和 SD 的积分计算：假设没有“直接”从基线状态到结局状态的转移（即 λ₀₂ 在所有情况下都设为零）。那么所有结局都必须经由中介事件后才发生（t_M < t_S 必然成立）。那此时： - 可观测的路径就只剩一条：基线→中介→结局。 - 在多状态模型下，这退化为一个串联两阶段模型。 - RD 变为：在未暴露组的中介分布 g₀ 下，暴露组 vs 未暴露组的生存概率差。因为未暴露组的中介时间分布是固定的，暴露组的生存概率完全由 λ₁₂ 中暴露A的条件效应（c₁）和中介时间M（c₂）来决定。这是一个“直接效应”。 - SD 变为：对于暴露组，将其中介分布从实际的 g₁ 切换到未暴露组的 g₀，生存概率的增量。这恰恰是“中介分布改变”带来的间接效应。 - 核心思路：在这最简特例下，论文的核心不是引入新的模型，而是提出一个计算框架：将 RD 和 SD 用多状态模型的转移特定危险和累积危险函数来表达成一个积分结果，进而通过 Cox 模型的输出（hazard ratios, baseline hazards 估计）和 bootstrap 来直接计算和推断。去掉竞争风险这个“壳”，这篇论文的核心数学操作就是对一个三阶段串行过程多状态模型的乃积分点估计与 bootstrap 推断。

三、这篇论文做了什么

三句话

① 开发了一个新的R命令 cmest_multistate，专门用于处理中介变量和结局均为时间事件且存在竞争风险的因果中介分析。② 核心方法基于多状态模型，先通过转移特定半参数Cox模型拟合数据，然后用公式将 RD 和 SD 两个随机干预因果估计量表达为基线风险与hazard ratio的积分函数，最终通过非参数bootstrap获得置信区间。③ 通过模拟数据证明：估计偏差低且覆盖率达到理论水平；该命令现已集成到CMAverseR包中，可直接供流行病学研究者使用。

关键设定与假设

  • 设定：暴露A二值，中介M是时间事件变量（可能有删失），结局S是时间事件变量（有竞争风险），基线协变量C。存在三个转移状态：0→1（暴露到中介事件）、0→2（暴露到结局事件）、1→2（中介事件到结局事件）。
  • 假设：核心假设包括：
  • 无未测量的暴露-中介/暴露-结局混杂：即给定C，A不对M或S产生未观测混淆。
  • 无未测量的中介-结局混杂：在给定C和A的条件下，M和S之间的因果路径未受未观测的U混淆。
  • 模型正确性：三个转移的Cox比例风险模型是同一时空的正确定义。
  • 与已有文献相比：相比Valeri et al. (2023)【5】，本文无新增理论假设——本质上是对该理论的封装实现。相比CMAverse【1】，其放宽了对中介和结局变量类型的限制（从连续/分类扩展到时间事件）。

主要结果

  • 理论结果：无。本文是应用/方法型论文，不包含新的统计定理或效率界。
  • 数值结果：基于一个n=1000的模拟数据。在给定时间点 s = {0.1, 0.5, 1, 2, 3, 4, 5, 6} 处，分别计算并比较了 cmest_multistate 给出的RD和SD的点估计与理论真值。结果显示（如图C所示）：
  • 低偏倚：1000次bootstrap抽样的均值几乎等同于理论真值。
  • 成功覆盖：95% bootstrap置信区间在所有时间点都成功包围了理论真值。
  • 无与任何baseline方法的对比。

🔎 结论是否比证明窄

  • 是。作者的结论“我们相信新命令将有助于健康差异研究”（最后一句），是基于一个特定的、理想化的模拟数据集。模拟中，λ₀₁₀等参数是常数（Weibull情形），且假设模型形式完全正确。在实际应用中，当Cox比例风险假设不成立或删失机制更复杂时，该结论的普遍性并未在理论或实证上被证明。作者自身也在最后一段承认：“在更强的无未测量混杂假设下……cmest_multistate可用于估计自然直接/间接效应。” 这种条件性的表述说明，论文的结论实际上比其声称的“直接工具价值”要窄——只有在核心识别假设成立、模型指定正确、数据规模足够大的情况下，该命令才保证有效。

证明路线与技术技巧（应用型论文：重点拆方法设计与实证）

  • 方法设计的整体路线（区别于定理证明，这里是建模与计算路线）：
  • 模型阶段：用 mstate 包拟合三个转移的半参数Cox模型，获得回归系数和基准风险函数的非参数估计。
  • 因果估计量表达：将RD和SD用转移特定危险函数和累积危险函数表达为积分。核心公式见Valeri et al. (2023)【5】。
  • 计算阶段：在给定的时间点s和协变量C的条件下，利用步骤1的估计值，对生存概率 P(S^{gₐ} > s) 进行数值积分计算。
  • 推断阶段：采用非参数 bootstrap（1000次重复抽样），分别对每一步都重复进行计算，从而得到点估计的分布并构造置信区间。
  • 关键细节：
  • 积分：估计量的计算本质上是数值积分，需要处理多个时间点上的累积危险函数。这是计算上最“吃劲”的部分，但作者未提及具体积分算法（如是否使用简单梯形法则还是Gauss-Legendre积分）。
  • Bootstrap：使用了并行计算来加速bootstrap（文中提到的“parallel-computed”）。这是保证软件实用性的重要工程细节。
  • 核心技巧：非参数 bootstrap + 并行计算。这是论文的核心推断技巧。亮点在于：避免了复杂的方差公式推导（如delta method或influence function），选择了一种“黑箱”但计算上可扩展的方法。

真实例子

  • 模拟数据：本文没有使用任何真实数据集，但通过一个精心设计的模拟实验来验证方法的有效性。
  • 数据来源：基于Cox模型人工合成数据（n=1000）。
  • 怎么用：用 cmest_multistate 命令计算特定时间点的RD和SD。
  • 结果说明：验证了命令的可复现性和精度。没有真实数据案例意味着无法评估该方法在存在模型错误指定或复杂删失模式下的稳健性——这是纯理论模拟的天然局限。
  • 结论：本文是纯应用/无真实数据例子。它明确说：“我们在模拟数据集上验证了cmest_multistate的有效性……”，而非真实数据。

四、开放问题（点到为止）

  • 开放问题 1：时间可变混杂（time-varying confounders）。本文和Valeri 2023【5】均在最后一句提到“在未来的工作中，我们计划允许时间可变混杂”。扎根点：原文最后一句 “In future works, we aim to allow for time-varying confounders...” 这需要更复杂的G-methods或结构化嵌套模型（如仿生II型模型）。
  • 开放问题 2：更强的识别假设验证与松弛。作者承认“在更强的无未测量混杂假设下……”才能估计自然直接/间接效应。扎根点：原文 “Under stronger no unmeasured confounding assumptions...”。目前该命令只支持RD/SD，不提供对“无混杂”假设的检验（如基于负对照变量的方法或敏感性分析）。研究者可探索如何将该命令与敏感性分析框架整合。
  • 开放问题 3：模型错误指定的影响。模拟中，λ₀₁₀ 等为常数（Weibull的简单情形）。当真实数据中的基准风险更复杂（如分片常数、非单调、或违反比例假设）时，该方法的稳健性未知。扎根点：模拟描述中只给出了“常数基准风险”的生成假设。没有讨论在模型错误指定下的性能。
  • 开放问题 4 (力学计算)：积分的数值精度与效率。本文未讨论数值积分的具体实现。当时间区间很细或维数较高时，积分计算的误差会如何影响置信区间的覆盖？扎根点：文章未描述其积分算法。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论