跳转至

Causal inference using multivariate generalized linear mixed-effects models

作者: Yizhen Xu, Ji Soo Kim, Laura K Hummers, Ami A Shah, Scott L Zeger
来源: Biometrics
主题: 因果推断
相关性: 9/10
链接: 期刊页 · arXiv

一、领域脉络与小综述

这个方向是什么： 这个子方向关注的是纵向观察性数据中的因果推断，核心难点在于如何处理时变混杂与未测量混杂的叠加问题。具体而言，在动态治疗方案下，研究者希望估计个体化或亚组层面的因果效应，但治疗分配机制未知、存在不可观测的个体异质性，且混杂因子随时间演化并受过去治疗影响。当前该领域已从经典的边际结构模型、g-估计等半参数方法，发展到利用混合效应模型、工具变量、敏感性分析等手段试图"部分识别"或"借力识别"未测量混杂的影响。成熟度方面，无未测量混杂假设下的纵向因果推断理论已相对成熟（MSM、SNM、g-computation），但引入未测量混杂后的识别与估计仍处于活跃探索期，尚无统一范式。

发展脉络： 1. 奠基工作：Robins 系列工作奠定了纵向因果推断的理论基石。Hernán & Robins (2010) 系统阐述了 g-formula、边际结构模型（MSM）与结构嵌套模型（SNM）在时变治疗下的识别与估计框架，核心假设是序列可忽略性——即给定过去历史后，治疗分配与潜在结果独立。这一假设排除了未测量混杂。本文引用此作为起点，明确指出"no unmeasured confounders"是现有主流方法的共同瓶颈。

  1. 主要进展——处理时变混杂：为处理可观测的时变混杂，MSM 与逆概率加权（IPTW）成为主流。Neugebauer et al. (2007) 提出历史限制边际结构模型（HRMSM），在降低维度的同时保持因果解释。Zhou, Elliott & Little (2019) 提出 PENCOMP，利用基于倾向性得分的样条插补，在 SUTVA、positivity、ignorability 下具有双稳健性质。这些工作都在"无未测量混杂"框架内打转。

  2. 主要进展——处理未测量混杂：

  3. 敏感性分析路线：Robins, Rotnitzky & Scharfstein (2000) 系统建立了未测量混杂的敏感性分析框架；Yang & Lok (2018) 针对粗粒化结构嵌套均值模型（SNMM），引入偏差函数量化未测量混杂的影响。本文作者明确引用此路线，但指出其属于"事后"敏感性分析，而非在主估计中直接建模。
  4. 工具变量路线：当存在有效工具变量时，可绕过未测量混杂。但本文指出，IV 往往不可得。

  5. 混合效应/随机效应路线：这是本文直接继承的线索。Shardell & Ferrucci (2018) 提出联合混合效应模型，将结果与治疗分配建模为共享随机效应，从而在"给定随机效应条件下序列可忽略"的假设下实现识别。本文引用此工作并指出其关键假设："no unmeasured confounders in the model... implies no treatment assignment heterogeneity"——即 Shardell & Ferrucci 的设定下，随机效应捕捉的是结果进程中的异质性，但治疗分配异质性被隐含排除了。

  6. 当前 Frontier 与本文位置：本文位于"随机效应路线"的前沿。作者明确要解决的问题是：当未测量混杂同时影响治疗分配与结果时，如何利用纵向数据的重复测量结构实现因果效应的部分识别与估计。相比 Shardell & Ferrucci (2018)，本文明确引入治疗分配异质性，将其建模为随机效应，并以此为基础提出新的识别条件——"条件于治疗分配异质性的序列可忽略性"。

子线索聚类： - 线索 A：边际结构模型与加权方法（MSM, IPTW, PENCOMP）。核心是构造伪人口，消除时变混杂影响。瓶颈：依赖无未测量混杂假设。 - 线索 B：结构嵌套模型与 g-估计（SNM, SNMM）。直接建模反事实结果与治疗历史的关系。瓶颈：同样依赖无未测量混杂假设；敏感性分析是事后补救。 - 线索 C：混合效应/随机效应模型（Shardell & Ferrucci, 本文）。利用纵向数据的重复测量结构，将未测量时间不变混杂参数化为随机效应，通过联合建模实现"条件独立性"。这是本文的主战场。 - 线索 D：固定效应与面板数据方法（Imai & Kim 2019, Allison et al. 2017）。利用个体固定效应消除时间不变混杂，但牺牲了动态因果结构。本文引用 Imai & Kim (2019) 的"瞬时处理效应"概念，但指出其无法处理长期效应与动态治疗方案。

这个方向在追问的核心问题： 1. 识别问题：在存在未测量时间不变混杂时，因果效应是否可识别？需要何种额外假设？如何量化识别的不确定性？ 2. 异质性问题：如何区分"治疗效应异质性"与"治疗分配偏好异质性"？两者如何与未测量混杂纠缠？ 3. 动态治疗方案问题：如何定义与估计动态治疗方案下的因果效应？如何处理治疗历史与时变混杂的复杂交互？ 4. 计算问题：在高维纵向设定下，贝叶斯后验推断的计算可行性如何？

⚠️ 作者的 framing： 作者将缺口 frame 为：现有方法要么假设无未测量混杂（MSM、SNM），要么依赖工具变量，要么在敏感性分析中事后处理未测量混杂；而混合效应路线虽有潜力，但先前工作（Shardell & Ferrucci 2018）隐含排除了治疗分配异质性。因此，本文的贡献被定位为：首次在混合效应框架下显式建模治疗分配异质性，并提出"条件于治疗分配异质性的序列可忽略性"假设，从而实现因果识别。

被淡化或回避的竞争路线： - 双稳健方法（如 augmented IPTW）：在无未测量混杂假设下具有稳健性，但无法处理未测量混杂。Intro 未深入讨论。 - 近端因果推断：利用代理变量处理未测量混杂，是近年热点。Intro 完全未提及，可能是作者刻意回避或认为不适用。 - 负控制方法（Negative Control）：另一条处理未测量混杂的路线，同样未提及。 - 主分层方法（Principal Stratification）：处理治疗后变量问题，与本文设定有交集但侧重点不同。

明显该被引但未出现的： - Proximal Causal Inference 系列工作（Tchetgen Tchetgen et al.）：利用负控制结果/暴露处理未测量混杂，与本文目标高度重合，但技术路线迥异。研究者应核查为何未被引用——是作者不了解，还是认为假设不兼容？ - M-estimation 与半参数效率理论：本文采用贝叶斯参数化方法，未讨论半参数效率界或无偏估计量。对于熟悉半参数理论的研究者，这是一个明显的缺口。

张力： 未见明显对立引用。但存在一个隐含张力：Shardell & Ferrucci (2018) 假设随机效应捕捉的是结果进程异质性，而本文假设随机效应捕捉的是治疗分配异质性——两者是否兼容？如果两者同时存在（结果异质性 + 分配异质性），模型如何设定？本文是否过度简化？这是研究者应追问的点。

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

第一步：符号、模型、可观测数据

符号定义： - \(i = 1, \ldots, n\)：个体下标，\(n\) 为样本量。 - \(t = 0, 1, \ldots, T\)：时间下标，\(T\) 为最终时间点。 - \(A_{it}\)：个体 \(i\) 在时间 \(t\) 的治疗变量（二值或连续）。 - \(\bar{A}_i = (A_{i0}, A_{i1}, \ldots, A_{iT})\)：治疗历史向量。 - \(Y_{it}\)：个体 \(i\) 在时间 \(t\) 的结果变量（连续或离散）。 - \(\bar{Y}_i = (Y_{i0}, Y_{i1}, \ldots, Y_{iT})\)：结果历史向量。 - \(L_{it}\)：个体 \(i\) 在时间 \(t\) 的时变混杂变量（向量）。 - \(\bar{L}_i = (L_{i0}, L_{i1}, \ldots, L_{iT})\)：时变混杂历史。 - \(b_i\)：个体随机效应（标量或向量），代表未测量的时间不变混杂因素。 - \(g\)：动态治疗方案，是一个从历史 \((\bar{L}_i, \bar{Y}_i)\) 到治疗决策的确定性规则。 - \(Y_i(g)\)：在治疗方案 \(g\) 下的潜在结果。 - \(\theta\)：因果效应参数（如亚组平均处理效应）。

模型（数据生成机制）： 作者假设一个多元广义线性混合效应模型（MGLMM），联合建模治疗分配、时变混杂与结果进程：

  1. 随机效应分布：\(b_i \sim N(0, D)\)，其中 \(D\) 为协方差矩阵。\(b_i\) 捕捉未测量的时间不变混杂。

  2. 治疗分配模型（给定随机效应与时变混杂）：

    \[A_{it} \mid b_i, \bar{A}_{i,t-1}, \bar{L}_{it}, \bar{Y}_{it} \sim \text{Bernoulli}(\text{logit}^{-1}(\alpha_0 + \alpha_1 b_i + \alpha_2 L_{it} + \alpha_3 Y_{i,t-1} + \ldots))\]
    关键点：\(b_i\) 进入治疗分配模型，意味着治疗分配存在异质性，且与未测量混杂相关。

  3. 结果进程模型（给定随机效应、治疗与时变混杂）：

    \[Y_{it} \mid b_i, \bar{A}_{it}, \bar{L}_{it}, \bar{Y}_{i,t-1} \sim \text{Normal}(\beta_0 + \beta_1 b_i + \beta_2 A_{it} + \beta_3 L_{it} + \ldots, \sigma^2)\]
    关键点：\(b_i\) 同时进入结果模型，意味着未测量混杂同时影响治疗与结果。

  4. 时变混杂进程模型（可选，若需完整联合分布）：

    \[L_{it} \mid b_i, \bar{A}_{i,t-1}, \bar{L}_{i,t-1}, \bar{Y}_{i,t-1} \sim \ldots\]

可观测数据： 研究者实际观测到的是 \(\{ (A_{it}, Y_{it}, L_{it}) : i = 1, \ldots, n; t = 0, \ldots, T \}\)，即每个个体在多个时间点的治疗、结果与时变混杂。不可观测的是 \(b_i\)（随机效应/未测量混杂）与潜在结果 \(Y_i(g)\)。

核心识别假设： 作者提出的"条件于治疗分配异质性的序列可忽略性"：

\[Y_i(g) \perp\!\!\!\perp A_{it} \mid b_i, \bar{A}_{i,t-1}, \bar{L}_{it}, \bar{Y}_{i,t-1}\]
直觉：在控制了随机效应 \(b_i\)（代表治疗分配偏好/异质性）后，治疗分配与潜在结果独立。这相当于"平衡了潜在的治疗偏好"。

第二步：最小内核

最简特例：\(T=1\)，单个时间点，二值治疗，连续结果，\(b_i\) 为标量

在这个特例下，模型简化为： - \(b_i \sim N(0, \sigma_b^2)\) - \(A_i \mid b_i \sim \text{Bernoulli}(\text{logit}^{-1}(\alpha_0 + \alpha_1 b_i))\) - \(Y_i(0), Y_i(1) \mid b_i \sim N(\beta_0 + \beta_1 b_i + \beta_2 A_i, \sigma^2)\)

要估计的因果效应：平均处理效应 \(\tau = E[Y_i(1) - Y_i(0)]\)。

识别问题： 在标准设定下，若无未测量混杂，\(\tau\) 可通过调整可观测混杂识别。但这里，\(b_i\) 同时影响 \(A_i\) 与 \(Y_i\)，且不可观测。因此，\(A_i\) 与 \(Y_i\) 之间存在混杂路径 \(A_i \leftarrow b_i \rightarrow Y_i\)，标准方法失效。

本文的解决方案： 1. 利用重复测量结构：假设每个个体有多个时间点的观测（在 \(T=1\) 特例下不适用，但推广到 \(T \geq 2\) 即可）。在 \(T \geq 2\) 时，通过联合建模 \(A_{i0}, A_{i1}, \ldots\) 与 \(Y_{i0}, Y_{i1}, \ldots\)，随机效应 \(b_i\) 可被"借力"识别——因为 \(b_i\) 解释了个体内相关性。 2. 贝叶斯 g-computation： - 步骤 1：基于 MGLMM，写出联合似然函数 \(p(\{A_{it}, Y_{it}, L_{it}\} \mid b_i, \theta)\)。 - 步骤 2：对 \(b_i\) 积分，得到边际似然 \(p(\{A_{it}, Y_{it}, L_{it}\} \mid \theta) = \int p(\{A_{it}, Y_{it}, L_{it}\} \mid b_i, \theta) p(b_i) db_i\)。 - 步骤 3：通过 MCMC 采样获得参数 \(\theta\) 的后验分布 \(p(\theta \mid \text{data})\)。 - 步骤 4：对于给定的治疗方案 \(g\)，计算潜在结果的预测分布：

\[p(Y_i(g) \mid \text{data}) = \int p(Y_i(g) \mid b_i, \theta) p(b_i \mid \text{data}) p(\theta \mid \text{data}) db_i d\theta\]
- 步骤 5：计算因果效应的后验分布，如 \(\tau(g) = E[Y_i(g) \mid \text{data}] - E[Y_i(g') \mid \text{data}]\)。

为什么这个最小内核能抓住核心： - 识别的数学本质：在 \(T \geq 2\) 的纵向设定下，随机效应 \(b_i\) 通过"个体内相关性"被部分识别。即使 \(b_i\) 不可观测，但多个时间点的数据提供了关于 \(b_i\) 分布的信息。这类似于随机效应模型中的"借力识别"（shrinkage estimation）。 - 因果推断的本质：在给定 \(b_i\) 的条件下，序列可忽略性成立，因此可以通过 g-computation 计算潜在结果。关键在于 \(b_i\) 不是直接观测的，而是通过贝叶斯后验推断"估计"出来的。 - 证明路线的雏形：一般情形的证明只是这个最小内核的"加壳"——将 \(T=1\) 推广到 \(T \geq 2\)，将单个治疗推广到动态治疗方案，将单个结果推广到纵向结果进程，将标量 \(b_i\) 推广到向量随机效应。

三、这篇论文做了什么

三句话： 1. 研究了纵向观察性数据中存在未测量时间不变混杂时，动态治疗方案下因果效应的识别与估计问题。 2. 核心方法是多元广义线性混合效应模型（MGLMM）联合建模治疗分配、时变混杂与结果进程，配合贝叶斯 g-computation 算法计算因果效应的后验分布。 3. 主要结论是：在"条件于治疗分配异质性的序列可忽略性"假设下，通过随机效应捕捉未测量混杂，可以实现因果效应的部分识别，模拟与真实数据验证了方法的有限样本性能。

关键设定与假设：

  1. 设定：
  2. 纵向观察性研究，\(n\) 个个体，\(T+1\) 个时间点。
  3. 每个个体在每个时间点观测到治疗 \(A_{it}\)、结果 \(Y_{it}\)、时变混杂 \(L_{it}\)。
  4. 存在未测量的时间不变混杂 \(b_i\)（随机效应）。

  5. 目标是估计动态治疗方案 \(g\) 下的因果效应，如亚组平均处理效应。

  6. 核心假设：

  7. 假设 1：条件序列可忽略性。

    \[Y_i(g) \perp\!\!\!\perp A_{it} \mid b_i, \bar{A}_{i,t-1}, \bar{L}_{it}, \bar{Y}_{i,t-1}, \quad \forall t, g\]
    统计含义：在控制了随机效应 \(b_i\)、过去治疗历史、当前时变混杂与过去结果后，当前治疗分配与潜在结果独立。这是一个强假设，因为它要求 \(b_i\) 捕捉了所有未测量混杂。 与已有文献的关系：相比标准序列可忽略性（假设无未测量混杂），本文的假设更弱——允许未测量混杂存在，但要求其被 \(b_i\) 捕捉。相比 Shardell & Ferrucci (2018)，本文明确引入治疗分配异质性，而非仅结果进程异质性。

  8. 假设 2：随机效应的正确设定。 \(b_i\) 服从特定分布（如正态分布），且进入治疗分配与结果模型的形式正确。 统计含义：这是一个参数化假设，若误设会导致偏差。本文未提供稳健性分析或敏感性分析来检验此假设的影响。

  9. 假设 3：MGLMM 模型的正确设定。 治疗分配模型、结果模型、时变混杂模型的函数形式正确（如线性、logit 连接函数等）。 统计含义：同样是参数化假设，误设风险存在。

  10. 假设 4：Positivity（正定性）。 在给定 \(b_i\) 与历史的条件下，每个治疗水平的概率非零。 统计含义：这是因果推断的标准假设，但在实践中可能被违反（如某些患者因禁忌症绝不接受某种治疗）。

主要结果：

  1. 定理 1：识别性结果（非正式陈述）。 在假设 1-4 下，动态治疗方案 \(g\) 下的潜在结果分布可识别，且可通过 g-formula 表达：

    \[p(Y_i(g) \mid \text{data}) = \int \prod_{t=0}^{T} p(Y_{it} \mid b_i, \bar{A}_{i,t-1}, A_{it}=g_t, \bar{L}_{it}, \bar{Y}_{i,t-1}, \theta) \times p(b_i \mid \text{data}) p(\theta \mid \text{data}) db_i d\theta\]
    直觉：这是 Robins g-formula 的贝叶斯版本，在给定 \(b_i\) 的条件下，通过积分消除历史依赖，得到潜在结果的预测分布。

  2. 算法：贝叶斯 g-computation。

  3. 输入：观测数据 \(\{ (A_{it}, Y_{it}, L_{it}) \}\)，MGLMM 模型设定，先验分布 \(p(\theta)\)。
  4. 步骤：
    1. MCMC 采样获得参数后验 \(p(\theta \mid \text{data})\) 与随机效应后验 \(p(b_i \mid \text{data})\)。
    2. 对于每个 MCMC 迭代 \(s\)，从 \(p(b_i \mid \theta^{(s)}, \text{data})\) 中抽取 \(b_i^{(s)}\)。
    3. 根据治疗方案 \(g\)，模拟潜在结果轨迹：
      \[Y_{it}^{(s)}(g) \sim p(Y_{it} \mid b_i^{(s)}, \bar{A}_{i,t-1}^{(s)}(g), A_{it}=g_t, \bar{L}_{it}^{(s)}(g), \bar{Y}_{i,t-1}^{(s)}(g), \theta^{(s)})\]
    4. 计算因果效应，如 \(\tau^{(s)}(g) = \frac{1}{n} \sum_{i=1}^{n} Y_{iT}^{(s)}(g)\)。

  5. 输出：因果效应的后验样本 \(\{ \tau^{(s)}(g) \}_{s=1}^{S}\)，用于计算后验均值、可信区间等。

  6. 模拟研究结果：

  7. 设定：生成纵向数据，包含未测量混杂 \(b_i\)，比较本文方法与标准方法（忽略未测量混杂的 MGLMM、MSM）。

  8. 结果：当存在未测量混杂时，标准方法有偏；本文方法在正确设定模型下偏差较小。当模型误设时（如 \(b_i\) 分布误设），偏差增大但仍优于完全忽略未测量混杂的方法。

  9. 真实数据应用：

  10. 数据：Johns Hopkins Scleroderma Center 的硬皮病患者纵向数据，\(n \approx 1000\)，\(T \approx 5\) 年。
  11. 治疗：霉酚酸酯（MMF）的持续使用。
  12. 结果：肺功能指标（FVC）。
  13. 时变混杂：疾病活动度、并发症等。
  14. 目标：估计不同亚组（如疾病严重程度分层）中 MMF 持续使用对肺功能的影响。
  15. 发现：MMF 持续使用对早期患者的肺功能保护效果更显著，与临床认知一致。

证明路线与技术技巧：

  1. 整体路线：
  2. 第一步：建立 MGLMM 联合模型，将治疗分配、结果进程、时变混杂建模为共享随机效应 \(b_i\) 的函数。
  3. 第二步：提出条件序列可忽略性假设，将因果识别问题转化为"给定 \(b_i\) 下的标准 g-computation"。
  4. 第三步：利用贝叶斯框架，通过 MCMC 获得参数与随机效应的后验分布。

  5. 第四步：通过 g-computation 计算潜在结果的预测分布，进而得到因果效应的后验分布。

  6. 关键跳跃点：

  7. 跳跃点 1：如何从观测数据识别 \(b_i\)？ 难点：\(b_i\) 是潜在变量，不可直接观测。 解决：利用纵向数据的重复测量结构。在 \(T \geq 2\) 时，个体内的相关性提供了关于 \(b_i\) 的信息。通过联合建模多个时间点的数据，\(b_i\) 的后验分布可被"借力"识别。 技术工具：贝叶斯分层模型 + MCMC 后验采样。

  8. 跳跃点 2：如何处理动态治疗方案下的反事实轨迹？ 难点：在治疗方案 \(g\) 下，时变混杂 \(L_{it}\) 的轨迹也是反事实的——因为 \(L_{it}\) 受过去治疗影响。 解决：采用 g-computation 的标准策略——递归模拟反事实轨迹。在每个时间点，根据 \(g\) 设定治疗，根据模型模拟 \(L_{it}\) 与 \(Y_{it}\)，然后进入下一时间点。 技术工具：Monte Carlo 模拟 + 递归反事实轨迹生成。

  9. 技术技巧点名：

  10. 贝叶斯分层模型：用于联合建模治疗、结果、时变混杂，共享随机效应 \(b_i\)。
  11. MCMC（Markov Chain Monte Carlo）：用于后验推断。具体采样方案未详细说明，可能使用 Gibbs 采样或 Hamiltonian Monte Carlo。
  12. g-computation：Robins 的经典方法，用于计算潜在结果分布。本文将其嵌入贝叶斯框架。
  13. 单世界干预图：用于可视化因果假设与识别条件。

真实例子与应用： - 数据场景：硬皮病患者的纵向观察性研究，关注免疫抑制剂 MMF 的持续使用对肺功能的影响。 - 方法应用： 1. 构建 MGLMM 模型，将 MMF 使用（二值）、肺功能（连续）、疾病活动度（时变混杂）联合建模，引入个体随机效应 \(b_i\)。 2. 设定动态治疗方案：\(g_1\) = 持续使用 MMF，\(g_0\) = 从不使用 MMF。 3. 通过贝叶斯 g-computation 计算每个亚组的平均处理效应后验分布。 - 结果：早期疾病患者中，MMF 持续使用对肺功能有显著保护效果；晚期患者效果不明显。这与临床认知一致——早期干预更有效。 - 说明什么：验证了方法的实用性，展示了如何从观察性数据中提取因果信号，同时控制未测量混杂。

🔎 结论是否比证明窄： - 本文的核心识别结果（定理 1）依赖于条件序列可忽略性假设与MGLMM 模型的正确设定。这两个假设都是不可检验的（untestable）。 - 作者在讨论部分承认：若 \(b_i\) 未能捕捉所有未测量混杂，或模型误设，估计会有偏。但未提供敏感性分析来量化这种偏差。 - 作者 claim 本文方法可"部分识别"未测量混杂的影响，但严格来说，这是在参数化模型下的"点识别"，而非集合识别。若模型误设，识别可能完全失效。 - 一个明显的缺口：本文未与 Proximal Causal Inference 或负控制方法进行比较，而这些方法同样处理未测量混杂问题。研究者应追问：在何种设定下，本文的随机效应方法优于或劣于 Proximal 方法？

四、开放问题

  1. 模型误设的敏感性分析：本文假设 MGLMM 模型正确设定，特别是随机效应 \(b_i\) 的分布与进入模型的形式。若 \(b_i\) 误设（如应为多峰分布而非正态），偏差有多大？如何构建敏感性分析框架？扎根点：Section 5 讨论部分提到 "model misspecification" 但未深入。

  2. 半参数拓展：本文采用参数化贝叶斯方法，能否拓展到半参数框架？例如，保持 \(b_i\) 的参数化假设，但对结果模型采用非参数或机器学习方法？这涉及半参数效率理论与 debiased ML。扎根点：Section 1 提到 "Bayesian g-computation"，但未讨论半参数方法。

  3. 与 Proximal Causal Inference 的比较：Proximal 方法利用代理变量处理未测量混杂，与本文的随机效应路线目标相同但技术迥异。两者在何种设定下等价或互补？扎根点：Intro 未引用 Proximal 相关文献，这是一个明显的缺口。

  4. 计算可扩展性：MCMC 在高维设定下（如 \(n\) 很大或 \(T\) 很大）的计算成本如何？能否用变分推断或其他近似方法加速？扎根点：Section 4 提到使用 MCMC，但未讨论计算成本或近似方法。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论