Causal inference using multivariate generalized linear mixed-effects models¶

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

一、领域脉络与小综述¶

这个方向是什么¶

本文所处理的根本问题是：在纵向观测研究中，如何利用贝叶斯方法对动态治疗方案（即随时间变化的治疗序列）进行亚组特异性因果效应预测。其核心挑战是未测时不变混杂——即存在某些不随时间变化的患者特征（如遗传倾向、基础健康状况），同时影响治疗分配和结果轨迹，且未被观测到。该方向当前处于“方法学应用扩展”阶段：已有大量识别策略（g-computation、MSM、SNM）和估计框架（贝叶斯、频率学派、双重稳健），但将未测时不变混杂通过随机效应纳入联合模型并用于动态治疗方案的亚组预测，仍是一个相对较新的尝试。

发展脉络¶

奠基工作：潜在结果框架与纵向因果推断的识别基础
- Rubin (1974)：定义了因果效应为同一组受试者在不同处理下的潜在结果比较，奠定了整个领域的语言基础。
- Hernán & Robins (2010); Richardson & Robins (2013)：发展了单世界干预图（SWIG），为纵向数据中序贯可忽略性等识别假设提供了图形化工具。本文用SWIG来展示其识别假设下的条件独立性。
- Robins, Rotnitzky, & Scharfstein (2000); Yang & Lok (2018)：提出了对未测混杂进行事后敏感性分析的方法。本文的定位是替代这种事后分析，转而将未测混杂的影响直接纳入模型。
主要进展：处理时变混杂的三大主流方法
- 边际结构模型 (MSM) 与逆概率加权 (IPTW)：如 Neugebauer et al. (2007) 的历史限制MSM。通过为每个时间点的治疗分配建模并加权，来模拟随机化。但权重不稳定时表现不佳。
- 结构嵌套模型 (SNM) 与 G-估计：如 Sitlani et al. (2012) 和 Qian et al. (2020) 的纵向结构混合模型。通过建模治疗“blip”效应来估计因果参数，对模型设定敏感。
- G-computation：通过模拟整个结果轨迹来估计因果效应。Schnitzer et al. (2019) 发展了纵向C-TMLE，实现了G-computation中的数据自适应模型选择。本文的核心方法正是贝叶斯G-computation。
当前 Frontier：处理未测混杂与个体异质性
- Shardell & Ferrucci (2018)：提出了联合混合效应模型，通过共享随机效应来连接结果和治疗选择模型，以处理未测混杂。本文引用其工作，并指出其假设“无未测混杂”等价于“无治疗分配异质性”。
- Imai & Kim (2019)：使用单位固定效应模型进行匹配，估计同期治疗效应。本文指出其方法排除了动态因果关系的可能性。
- Zhou et al. (2019) (PENCOMP)：提出了基于惩罚样条的倾向性预测方法，用于缺失潜在结果，具有双重稳健性。本文将其归为“结合倾向性得分和结果模型”的贝叶斯方法之一。
- 本文的位置：作者声称，他们的方法通过多变量广义线性混合效应模型 (MGLMM) 将未测时不变因素作为随机效应，同时纳入结果、时变混杂和治疗分配的联合分布中，从而部分识别了未测的永久性成分。这被视为对 Shardell & Ferrucci (2018) 的推广——后者假设未测混杂与结果异质性成比例，而本文允许它们通过随机效应更灵活地相关。

子线索聚类¶

贝叶斯因果推断：以 Li, Ding, & Mealli (2022) 的综述为代表，讨论了倾向性得分在贝叶斯框架下的角色、可识别性的定义、先验选择等。本文属于这一簇，强调贝叶斯方法在处理复杂模型和不确定性量化上的优势。
联合模型 / 共享参数模型：以 Shardell & Ferrucci (2018) 和 Achana et al. (2021) 为代表。核心思想是通过随机效应或共享参数，将结果模型和治疗/混杂模型联合起来，以处理未测混杂。本文的MGLMM属于这一簇。
动态治疗方案与精准医学：以 Kosorok & Laber (2019) 和 Rosen & Zeger (2019) 为代表。关注点在于如何利用观测数据为个体患者推荐最优的治疗序列。本文的应用场景（硬皮病患者的MMF用药）直接服务于这一目标。

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

识别问题：在纵向观测研究中，如何放松“无未测混杂”这一强假设？本文的答案是：通过将未测混杂建模为随机效应，并假设条件于治疗分配异质性的序贯可忽略性。
估计问题：如何高效、稳定地估计复杂动态治疗方案下的因果效应？本文的答案是：贝叶斯G-computation，通过MCMC从后验分布中抽样，自然处理了模型不确定性。
个性化问题：如何估计亚组特异性的因果效应，而非仅仅平均效应？本文的答案是：通过MGLMM中的随机效应和协变量交互，直接输出每个亚组的后验预测分布。

⚠️ 作者的 framing¶

作者把缺口 frame 成什么：作者声称，现有方法（如Shardell & Ferrucci, 2018）在处理未测混杂时，隐含地假设了“无治疗分配异质性”。本文通过MGLMM，允许未测时不变因素对治疗分配和结果轨迹有不同的影响，从而更灵活地处理未测混杂。作者将这一推广定位为“显然的下一步”。
哪些竞争路线被他淡化或回避了：
- 双重稳健估计：作者在引言中提到了“后验逆概率加权或双重稳健估计”（Li et al., 2022），但本文并未采用或比较双重稳健方法。作者选择了一个纯模型依赖的贝叶斯G-computation路线，回避了双重稳健性带来的理论复杂性。
- 工具变量法：作者在引言中提到了“未混杂的工具变量”，但本文并未使用。对于处理未测混杂，IV是另一种强大的策略，但本文完全回避了这一路线。
什么明显该被引 / 该存在、却没出现在 intro 里？
- 深度学习方法：近年来，用循环神经网络（RNN）或Transformer进行纵向因果推断的工作（如 Bica et al., 2020 “Estimating counterfactual treatment outcomes over time through adversarially balanced representations”）没有被引用。这可能是由于本文是纯贝叶斯参数模型，与深度学习范式差异较大，但作为处理复杂纵向数据的竞争方法，其缺失值得注意。
- 更近期的贝叶斯非参数方法：例如用高斯过程（GP）建模结果轨迹的贝叶斯因果推断工作没有被引用。这可能是由于本文的MGLMM框架是参数化的，但GP作为处理个体异质性的另一种强大工具，其缺失是一个信号。

张力¶

未见明显对立引用。所有被引工作基本在同一个“处理未测混杂”的共识下，以不同的技术路线（联合模型、MSM、SNM、敏感性分析）进行探索。本文与Shardell & Ferrucci (2018) 的关系是“推广”而非“对立”。

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

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

符号：
- \(i = 1, \dots, N\)：个体索引。
- \(t = 1, \dots, T_i\)：时间点索引（每个个体观测次数可能不同）。
- \(Y_{it}\)：个体 \(i\) 在时间 \(t\) 的结果变量（如肺功能指标FVC）。这是可观测的。
- \(A_{it}\)：个体 \(i\) 在时间 \(t\) 的治疗分配（如是否使用MMF）。这是可观测的。
- \(L_{it}\)：个体 \(i\) 在时间 \(t\) 的时变混杂（如疾病活动度评分）。这是可观测的。
- \(\bar{A}_{i,t} = (A_{i1}, \dots, A_{it})\)：个体 \(i\) 截至时间 \(t\) 的治疗历史。
- \(\bar{L}_{i,t} = (L_{i1}, \dots, L_{it})\)：个体 \(i\) 截至时间 \(t\) 的混杂历史。
- \(\bar{Y}_{i,t-1} = (Y_{i1}, \dots, Y_{i,t-1})\)：个体 \(i\) 截至时间 \(t-1\) 的结果历史。
- \(b_i\)：个体 \(i\) 的随机效应向量。这是潜在/不可观测的，代表未测时不变因素。
- \(Y_{it}(\bar{a})\)：在治疗序列 \(\bar{a}\) 下的潜在结果。这是不可观测的（除了在观测到的治疗序列下）。
- estimand：对于给定的动态治疗方案 \(q\)（定义了一组治疗分配规则），我们关心的是亚组特异性的潜在结果均值，例如 \(E[Y_{i,T}(\bar{a}^q) | X_i = x]\)，其中 \(X_i\) 是基线协变量。
模型：
- 作者假设一个多变量广义线性混合效应模型 (MGLMM)，它联合建模了结果、时变混杂和治疗分配。核心是，这三个过程都依赖于同一个个体随机效应 \(b_i\)。
- 具体来说，对于每个时间点 \(t\)，模型可以写成：
  - 结果模型：\(g_Y(E[Y_{it} | \bar{Y}_{i,t-1}, \bar{L}_{i,t}, \bar{A}_{i,t}, b_i]) = \text{一些线性预测项}\)
  - 混杂模型：\(g_L(E[L_{it} | \bar{Y}_{i,t-1}, \bar{L}_{i,t-1}, \bar{A}_{i,t-1}, b_i]) = \text{一些线性预测项}\)
  - 治疗模型：\(g_A(E[A_{it} | \bar{Y}_{i,t-1}, \bar{L}_{i,t}, \bar{A}_{i,t-1}, b_i]) = \text{一些线性预测项}\)
- 其中 \(g_Y, g_L, g_A\) 是连接函数（如恒等、logit）。关键假设是，给定 \(b_i\) 和过去的历史，当前的结果、混杂和治疗是条件独立的。\(b_i\) 捕捉了所有未测时不变混杂。
可观测数据：
- 研究者实际能观测到的是：\(\{Y_{it}, A_{it}, L_{it}, X_i\}_{i=1, t=1}^{N, T_i}\)。
- 想要但观测不到的是：
  1. 潜在结果 \(Y_{it}(\bar{a})\) 在非观测的治疗序列下。
  2. 随机效应 \(b_i\)。
  3. 未测时不变混杂（即 \(b_i\) 所代表的所有因素）。

第二步：讲最小内核¶

本文的最小内核可以简化为一个两时间点、二值治疗、连续结果的特例。

最简特例设定：
- \(T = 2\)：只有两个时间点 \(t=1, 2\)。
- \(A_{it} \in \{0, 1\}\)：治疗是二值的。
- \(Y_{it} \in \mathbb{R}\)：结果是连续的，且 \(g_Y\) 是恒等连接（线性模型）。
- 没有时变混杂 \(L_{it}\)（或假设 \(L_{it}\) 被 \(b_i\) 完全捕捉）。
- 我们关心的是动态治疗方案 \(\bar{a} = (a_1, a_2)\)，例如“始终治疗” \((1,1)\) 和“从不治疗” \((0,0)\)。
在这个特例下，模型退化为：
- 结果模型：
  \[Y_{i1} = \beta_0 + \beta_1 A_{i1} + b_i + \epsilon_{i1}\]
  
  \[Y_{i2} = \beta_0 + \beta_1 A_{i1} + \beta_2 A_{i2} + \beta_3 Y_{i1} + b_i + \epsilon_{i2}\]
  其中 \(b_i \sim N(0, \sigma_b^2)\) 是个体随机效应，\(\epsilon_{it} \sim N(0, \sigma_\epsilon^2)\) 是独立误差。
- 治疗模型（为了简化，假设治疗只依赖于 \(b_i\)，不依赖于过去结果）：
  \[\text{logit}(P(A_{i1}=1 | b_i)) = \gamma_0 + \gamma_1 b_i\]
  
  \[\text{logit}(P(A_{i2}=1 | A_{i1}, Y_{i1}, b_i)) = \delta_0 + \delta_1 A_{i1} + \delta_2 Y_{i1} + \delta_3 b_i\]
核心思路：
1. 识别假设：作者声称，通过将 \(b_i\) 同时放入结果和治疗模型，我们条件于 \(b_i\) 后，治疗分配是序贯可忽略的。即：
  \[Y_{it}(a_1, a_2) \perp A_{it} | \bar{Y}_{i,t-1}, \bar{A}_{i,t-1}, b_i\]
  这等价于说，所有未测时不变混杂都通过 \(b_i\) 被捕捉了。
2. 贝叶斯 G-computation：
  - 第一步（拟合模型）：用MCMC从观测数据中拟合上述MGLMM，得到所有参数（\(\beta, \gamma, \delta, \sigma_b^2, \sigma_\epsilon^2\)）和随机效应 \(b_i\) 的后验分布。
  - 第二步（模拟潜在结果）：对于每个MCMC后验样本，对于每个个体 \(i\)，我们强制设定一个治疗序列 \(\bar{a} = (a_1, a_2)\)。然后，利用拟合好的模型，从后验预测分布中模拟出在该治疗序列下的潜在结果轨迹：
    - 从后验中抽取 \(b_i^{(s)}\)。
    - 设定 \(A_{i1} = a_1\)。
    - 从 \(P(Y_{i1} | A_{i1}=a_1, b_i^{(s)})\) 中抽取 \(Y_{i1}^{(s)}\)。
    - 设定 \(A_{i2} = a_2\)。
    - 从 \(P(Y_{i2} | A_{i1}=a_1, A_{i2}=a_2, Y_{i1}^{(s)}, b_i^{(s)})\) 中抽取 \(Y_{i2}^{(s)}\)。
  - 第三步（计算因果效应）：对于每个后验样本 \(s\)，计算亚组（例如，基线协变量 \(X_i\) 相同的个体）的平均潜在结果 \(\bar{Y}_T^{(s)}(\bar{a})\)。所有后验样本的分布就构成了该亚组在治疗方案 \(\bar{a}\) 下的后验预测分布。比较不同治疗方案的后验分布，即可得到因果效应。
为什么这个特例抓住了核心：
- 它清晰地展示了随机效应 \(b_i\) 如何作为“未测混杂的代理”，通过联合建模来满足识别假设。
- 它展示了贝叶斯G-computation的“预测-模拟”流程：先拟合模型，再强制干预、模拟轨迹。
- 它去掉了多变量、多类型数据、复杂协方差结构等一般性设定，让读者一眼看清“用随机效应处理未测混杂”这一核心想法。

三、这篇论文做了什么¶

三句话¶

研究了什么问题：在纵向观测研究中，如何估计动态治疗方案的亚组特异性因果效应，同时处理未测时不变混杂。
核心工具 / 方法：提出了一个多变量广义线性混合效应模型 (MGLMM) 来联合建模结果、时变混杂和治疗分配，并开发了一个贝叶斯 G-computation 算法来计算亚组特异性干预效益的后验分布。
主要结论：通过模拟研究和真实数据（硬皮病患者持续使用MMF）分析，作者表明该方法能够有效处理由未测时不变因素导致的治疗分配异质性，并提供合理的亚组特异性因果效应估计。

关键设定与假设¶

设定：纵向观测数据，每个个体有多个时间点的观测。数据包括结果 \(Y_{it}\)、治疗 \(A_{it}\)、时变混杂 \(L_{it}\) 和基线协变量 \(X_i\)。
假设：
1. SUTVA (Stable Unit Treatment Value Assumption)：个体间无交互，且治疗水平唯一。这是标准假设。
2. Positivity (Positivity)：对于任何可能的治疗历史和混杂历史，每个个体接受任何治疗的概率都大于0。这是G-computation的标准要求。
3. 序贯可忽略性 (Sequential Ignorability)：这是本文的核心假设。作者将其表述为条件于治疗分配异质性的序贯可忽略性。具体来说，给定过去的历史和个体随机效应 \(b_i\)，当前的治疗分配与未来的潜在结果独立。即：
  \[Y_{it}(\bar{a}) \perp A_{it} | \bar{Y}_{i,t-1}, \bar{L}_{i,t}, \bar{A}_{i,t-1}, b_i\]
  这比标准的“无未测混杂”假设更强，因为它要求所有未测混杂都被 \(b_i\) 捕捉。但作者认为，通过MGLMM对 \(b_i\) 的灵活建模，这个假设比Shardell & Ferrucci (2018) 中“无治疗分配异质性”的假设更弱、更合理。
4. 模型正确设定 (Correct Model Specification)：MGLMM的联合分布被假定为正确。这是贝叶斯参数方法的固有弱点，也是与双重稳健方法的主要区别。

主要结果¶

理论结果：本文没有提供新的渐近理论或效率界。其理论贡献在于识别：作者证明了在MGLMM框架下，条件于随机效应的序贯可忽略性假设足以识别动态治疗方案的因果效应。这个证明是通过SWIG和条件独立性论证完成的，而非严格的数学定理。
模拟研究：
- 设计：模拟了两种场景：一种是数据由MGLMM生成（符合模型假设），另一种是数据由更复杂的机制生成（模型误设定）。
- 核心量化结论：在模型正确设定时，本文提出的贝叶斯G-computation方法能够无偏地估计亚组特异性因果效应，且后验区间覆盖率达到名义水平。在模型误设定时，估计结果出现偏差，但作者声称其偏差程度“可接受”。
- 与 baseline 对比：模拟中与“朴素”的线性混合模型（忽略治疗分配机制）进行了对比。结果显示，朴素方法在有未测混杂时存在严重偏差，而本文方法能有效校正。
真实数据应用：
- 数据：来自约翰霍普金斯硬皮病中心研究数据库的观测数据。关注的是霉酚酸盐 (MMF) 对硬皮病患者肺功能（FVC）的疗效。
- 怎么用：将患者按基线特征（如年龄、疾病亚型）分为亚组。对于每个亚组，用贝叶斯G-computation估计“持续使用MMF” vs. “从未使用MMF”两种动态治疗方案下的FVC轨迹。
- 结果：结果显示，MMF对某些亚组（如弥漫性皮肤型硬皮病患者）的肺功能下降有显著延缓作用，而对其他亚组效果不明显。这展示了方法的个性化预测能力。
- 这个例子想说明什么：验证了方法在真实世界复杂数据中的可行性，并展示了其核心卖点——能够发现异质性治疗效应，为精准医学决策提供依据。

证明路线与技术技巧¶

整体路线：
1. 模型构建：定义MGLMM，明确结果、混杂、治疗三个子模型的连接函数和随机效应结构。
2. 识别论证：利用SWIG，在MGLMM的假设下，展示条件于随机效应 \(b_i\) 后，治疗分配与潜在结果独立。这是G-computation有效性的理论基础。
3. 贝叶斯推断：为所有参数和随机效应设定先验分布。利用MCMC（如Gibbs采样或Hamiltonian Monte Carlo）从后验分布中抽样。
4. G-computation：对于每个MCMC后验样本，对每个个体，在指定的动态治疗方案下，从后验预测分布中模拟整个结果轨迹。
5. 汇总：对每个亚组，汇总所有后验样本的模拟结果，得到亚组特异性因果效应的后验分布。
关键跳跃点：
- 从“无未测混杂”到“条件于随机效应的序贯可忽略性”：这是本文的核心概念跳跃。作者没有证明这个假设比标准假设更弱，而是通过论证“未测时不变混杂”可以被随机效应捕捉，来为这个假设的合理性辩护。这个跳跃是概念性的，而非数学性的。
- 从联合模型到G-computation：如何从拟合好的MGLMM中“读出”因果效应？作者没有推导出闭式解，而是采用了模拟的方法。这个跳跃是计算性的，依赖于MCMC的抽样能力。
技术技巧点名：
- MCMC (Markov Chain Monte Carlo)：用于从复杂的后验分布中抽样。这是贝叶斯推断的标准工具。
- SWIG (Single-World Intervention Graph)：用于图形化地展示条件独立性假设，辅助识别论证。这是一个图形化工具，而非计算工具。
- 后验预测分布 (Posterior Predictive Distribution)：用于模拟在指定干预下的潜在结果。这是贝叶斯预测的核心概念。

🔎 结论是否比证明窄¶

是。作者在引言和摘要中声称方法能处理“未测时不变混杂”，但证明（识别论证）完全依赖于“条件于随机效应的序贯可忽略性”这一假设。这个假设本质上等价于“所有未测时不变混杂都被随机效应捕捉”。如果存在时变的未测混杂，或者未测混杂与随机效应的关系是非线性的，这个假设就会失效。因此，结论的适用范围被证明所依赖的假设严格限制。
作者在模拟研究中测试了模型误设定的情况，但并未在理论上证明方法的稳健性。因此，“能处理未测混杂”这一结论，在严格意义上，只适用于“未测混杂是时不变的且能被线性随机效应捕捉”的情形。作者在文中（如讨论部分）可能承认了这一点，但引言中的表述可能给读者更宽泛的印象。

四、开放问题¶

放松线性随机效应假设：本文假设未测混杂通过线性随机效应影响结果和治疗。如何将模型扩展到非线性随机效应（如高斯过程）或非参数随机效应分布？这扎根于本文对MGLMM的线性假设。
处理时变未测混杂：本文只处理了时不变的未测混杂。如果存在随时间变化的未测混杂（如患者情绪波动），模型将失效。如何扩展框架以处理时变未测混杂？这扎根于本文对“未测时不变因素”的强调。
与双重稳健方法的结合：本文的贝叶斯G-computation完全依赖于模型正确设定。如何将双重稳健性（如TMLE）的思想融入贝叶斯框架，以提供对模型误设定的保护？这扎根于本文回避了双重稳健方法这一事实。
计算效率与可扩展性：MCMC在高维随机效应和大规模数据下计算成本高昂。如何开发更高效的变分贝叶斯或近似贝叶斯计算方法，以扩展本文方法到更大规模的数据集？这扎根于本文对MCMC的依赖。

Maintained by 陈星宇 · Homepage · Source on GitHub