[论文笔记] 因果模型：边缘结构模型MSM

阅读James M. Robins的文章Marginal Structural Models and Causal Inference in Epidemiology^[1]后的笔记

---

基本概念

  边缘结构模型（Marginal structural models，MSMs）是一种因果模型，用于调整影响随时间变化的处理/治疗方案（time-varying treatments）的时依性混杂因子（time-dependent confounding），使得我们能无偏估计因果效应（unbiased estimation of casual effects）^[2]。

  传统的估计时变性处理的因果效应的方法是对结果的概率建立回归方程，将过去的处理（past treatment/exposure）和过去的混杂因子（past confounder）作为方程的变量，通常使用分层法等（比如逻辑回归或者Cox比例风险回归模型）。这些方法，无论是否对混杂因子进行调整，都会存在偏倚（biased）。所以提出了MSM模型。

  时依性混杂因子的定义：⑴是时依性的协变量，能够用于预测目标事件（是目标事件的风险因子），且能预测处理；⑵过去的预测方案能预测此变量；⑶自身过去的状态能预测现在的状态。

MSM的基本思想

  记时间为$k$，将可观测的对结果有影响的风险因子记为$L_k$，将不可观测（即无法收集这个变量的值）的对结果有影响的风险因子记为$U_k$，将时变性处理记为$A_k$，将结果记为$Y$。MSM假设不存在影响处理的不可观测风险因子，即不可测假设（untestable assumption），因此MSM假设下的因果图应该如图(a)所示。

图(a) 左侧为多个时间构成的因果图，右侧为单时刻的因果图（未去混杂）

  MSM使用逆处理概率加权法（Inverse Probability of Treatment Weighting，IPTW）消去可观测的混杂因子对处理的影响。也就是说采用IPTW生成各对象组的伪分布，每组处理（treatment）的伪分布都相同，使得处理相对于混杂因子独立，反映在图上就是删除了从混杂因子到处理的弧，处理后的因果图如图(b)所示。

图(b) 左侧为多个时间构成的因果图，右侧为单时刻的因果图（已去混杂）

因果效应的估计函数

  取单个时刻（如图(b)右侧子图所示）为例，假设取值都是二元的，处理不受混杂因子影响。使用粗风险差/粗危险差/粗率差（crude risk difference，RD）、粗风险比/粗相对危险度/粗危险比/粗率比（crude risk ratio，RR）、粗优势比/粗比值比（crude odds ratio，OR）估计处理对结果的影响（因果效应）。三者的计算公式如下：
$$cRD=pr[Y=1|A_0=1]-pr[Y=1|A_0=0]\text{\hspace{1em}(1)}$$
$$cRR=\frac{pr[Y=1|A_0=1]}{pr[Y=1|A_0=0]}\text{\hspace{1em}(2)}$$
$$cOR=\frac{pr[Y=1|A_0=1]/pr[Y=1|A_0=0]}{pr[Y=0|A_0=1]/pr[Y=0|A_0=0]}\text{\hspace{1em}(3)}$$
  前者是从观测数据角度得到因果的。从因果分析角度出发，处理的因果对比形式与以上计算公式相同，但是涉及了反事实（counterfactual）变量的概念。将被处理后的结果记为$Y_{a_0=1}$，将未受到处理的结果记为$Y_{a_0=0}$，以上两者无法被同时观测到，这就是反事实的概念。则这个个体的因果效应为$Y_{a_0=1}-Y_{a_0=0}$。相应的，因果风险差（causal risk difference）、因果风险比（causal risk ratio）、因果优势比（causal odds ratio）的计算公式如下：
$$causalRD=pr[Y_{a_0=1}=1]-pr[Y_{a_0=0}=1]\text{\hspace{1em}(4)}$$
$$causalRR=\frac{pr[Y_{a_0=1}=1]}{pr[Y_{a_0=0}=1]}\text{\hspace{1em}(5)}$$
$$causalOR=\frac{pr[Y_{a_0=1}=1]/pr[Y_{a_0=0}=1]}{pr[Y_{a_0=1}=0]/pr[Y_{a_0=0}=0]}\text{\hspace{1em}(6)}$$
  当处理不受混杂因子影响时，以上六个公式两两相等。

单时刻建模以及MSM释义

  $causalRD$、$causalRR$和$causalOR$可以分别表示为自变量是处理$a_0$的三种模型，分别为线性（linear）、对数线性（log linear）和线性逻辑模型（linear logistic model），考虑单时刻情况，将(7-9)记为因果模型：
$$pr[Y_{a_0}=1]={\psi }_0+{\psi }_1{a}_0\text{\hspace{1em}(7)}$$
$$\log pr[Y_{a_0}=1]={\theta }_0+{\theta }_1{a}_0\text{\hspace{1em}(8)}$$
$$logitpr[Y_{a_0}=1]={\beta }_0+{\beta }_1{a}_0\text{\hspace{1em}(9)}$$
  分别将$a_0=1$和$a_0=0$代入(7-9)式，再代入到(4-6)式，可以得到$causalRD={\psi }_1$、$causalRR=e^{{\theta }_1}$、$causalOR=e^{{\beta }_1}$。将(7-9)记为饱和MSM模型（saturated MSMs）。

  同样的，从观测数据角度可以得到粗风险差、粗风险比和粗优势比的饱和线性、对数和线性逻辑模型，将(10-12)记为数据模型：
$$pr[Y=1|A_0=a_0]={\psi }_0^{{}^{\prime }}+{\psi }_1^{{}^{\prime }}{a}_0\text{\hspace{1em}(10)}$$
$$\log pr[Y=1|A_0=a_0]={\theta }_0^{{}^{\prime }}+{\theta }_1^{{}^{\prime }}{a}_0\text{\hspace{1em}(11)}$$
$$logitpr[Y=1|A_0=a_0]={\beta }_0^{{}^{\prime }}+{\beta }_1^{{}^{\prime }}{a}_0\text{\hspace{1em}(12)}$$
  可以得到$cRD={\psi }_1^{{}^{\prime }}$、$cRR=e^{{\theta }_1^{{}^{\prime }}}$、$cOR=e^{{\beta }_1^{{}^{\prime }}}$。(10-12)表示的模型是为观测到的数据关联关系构建的模型，因此仅当处理是非混淆的情况下，因果模型(7-9)的参数才与数据模型(10-12)的参数相同。

  对边缘结构模型MSM的解释如下：

• 边缘性（marginal）：对反事实随机变量$Y_{a_0=1}$和$Y_{a_0=0}$建立了边缘分布（(7-9)式左半边）而不是联合分布，也就是不对$Y_{a_0=1}$和$Y_{a_0=0}$之间的相关性建模。
• 结构化（structural）：对反事实变量的概率进行建模，在计量经济学和社会科学中常常把反事实变量的建模称为结构化。
• 饱和性（saturated）：$pr[Y_{a_0=1}=1]$和$pr[Y_{a_0=0}=1]$是两种未知的概率，模型也有两个未知的参数，因此模型未对两种未知的概率值进行限制。

逆处理概率加权法（IPTW）

  当处理存在混淆时（受混杂因子影响），因果模型的参数将与数据模型的参数不相等。MSM模型假设不存在不可观测的混杂因子（No unmeasured confounders），那么就可以采用加权分析方法去除混杂因子的影响。

  设$i$为种类（Subject）编号，每个种类的权重记为式(13)。
$$w_i=\frac{1}{pr[A_0=a_{0i}|L_0=l_{0i}]}\text{\hspace{1em}(13)}$$
  $l_{0i}$是第$i$类的$L_0$观测值。$w_i$的值可以通过式(14)使用回归方法进行估计，统计不同的$a_0$和$l_0$并回归得到参数。
$$logitpr[A_0=a_0|L_0=l_0]={\alpha }_0+{\alpha }_1{l}_0\text{\hspace{1em}(14)}$$
  当取$A_0=1$时，$w_i=1+\exp (-{\hat{a}}_0-{\hat{a}}_1{l}_{0i})$；当取$A_0=0$时，$w_i=1+\exp ({\hat{a}}_0+{\hat{a}}_1{l}_{0i})$。这种方法能消除混杂因子影响的原因在于，IPTW为每类复制$w_i$份，形成伪总体分布，这种情况下，⑴相同的$A_0$取值时，任意的$L_0$取值概率都相同，因此$A_0$是非混淆的；⑵因为$A_0$是非混淆的，因此在上述伪分布数据中，数据模型得到的因果效应将与因果模型得到的结果一致。

从二分类到多分类：多层处理与非饱和MSM模型

  当处理多层非二值（multilevel treatment）、是长为$N$的序列值（比如为服药剂量，0~15mg共$N=16$种取值，且剂量是线性变化的）时，相应的潜在结果（potential outcomes）也将有多种取值（比如16种）。这种情况下，为了构造饱和模型必须设置多个参数（比如16个），因此无法再使用饱和模型进行建模。
  为了克服这种问题，假设处理效果是线性变化的（即简化剂量反应关系，parsimonious dose-response relationship），那么可以将因果模型写作非饱和的式(15)。
$$logitpr[Y_{a_0}=1]={\beta }_0+{\beta }_1{a}_0\text{\hspace{1em}(15)}$$
  其中，${\beta }_1$是斜率参数。当剂量增加1时，$causalOR$增加$e^{{\beta }_1}$。对应的数据模型就可以写作式(16)。
$$logitpr[Y=1|A_0=a_0]={\beta }_0^{{}^{\prime }}+{\beta }_1^{{}^{\prime }}{a}_0\text{\hspace{1em}(16)}$$
  与前面的分析相同，当处理是非混淆的情况下，两者估计的参数是相同的。当处理仅受可观测混杂因子$L_0$影响时，可以使用IPTW调节种类分布达到去混杂的效果。对于序列变量，$w_i=1/pr[A_0=a_{0i}|L_0=l_{0i}]$可由式(17)进行估计。
$$pr[A_0=a_0|L_0=l_0]=\frac{\exp ({\alpha }_{0a_0}+{\alpha }_1{l}_0)}{1+{\Sigma }_{j=1}^N\exp ({\alpha }_{0j}+{\alpha }_1{l}_0)},a_0=1,\dots ,N;$$
$$pr[A_0=0|L_0=l_0]=\frac{1}{1+{\Sigma }_{j=1}^N\exp ({\alpha }_{0j}+{\alpha }_1{l}_0)},a_0=1,\dots ,N;\text{\hspace{1em}(17)}$$
  式(16)可以理解为类Softmax多元逻辑回归函数，在不处理时（剂量为0）分子设为常量1。

稳定权重（Stabilized Weights）

  当某些处理状态$A_0$与混杂因子$L_0$高度相关时，很有可能某些状态组合会缺乏观测数据。这会导致样本数量非常少，继而由$w_i$调整得到的伪总体分布中该部分占比非常大，会影响分析效果。稳定权重$s{w}_i$记为式(18)。
$$s{w}_i=\frac{pr[A_0=a_{0i}]}{pr[A_0=a_{0i}|L_0=l_{0i}]}\text{\hspace{1em}(18)}$$
  估计$s{w}_i$的值需要计算分子和分母两部分，分母可以采用式(16)计算，分子计算公式见式(19)。
$$pr[A_0=a_0]=\frac{\exp ({\alpha }_{0a_0}^{\ast })}{1+{\Sigma }_{j=1}^N\exp ({\alpha }_{0j}^{\ast })},a_0=1,\dots ,N;$$
$$pr[A_0=0]=\frac{1}{1+{\Sigma }_{j=1}^N\exp ({\alpha }_{0j}^{\ast })},a_0=1,\dots ,N;\text{\hspace{1em}(19)}$$
  ${\alpha }_{0a_0}^{\ast }$的星号表明当$A_0$是混淆的情况时此参数与${\alpha }_{0a_0}$不相同。这是因为在计算${\alpha }_{0a_0}$时是以不同的$L_0$作为条件的，也就是说分别对不同的子集进行计算，当存在混淆时，各子集的分布不相同，与总体的分布自然不同。

从分类到离散：连续处理下的稳定权重

  当处理变量是连续的情况下，$w_i$的方差趋于无穷，因此不能使用。假设$A_0$服从正态分布，则$s{w}_i$可写作式(20)，其中$f(\ast )$是概率密度函数。
$$s{w}_i=\frac{f(a_{0i)}}{f(a_{0i}|l_{0i})}\text{\hspace{1em}(20)}$$
  为了估计$f(a_{0i}|l_{0i})$，给定$L_0$，$A_0\sim N({\alpha }_0+{\alpha }_1{L}_0,{\sigma }^2)$，因此$f(a_{0i}|l_{0i})$可表示为式(21)；为了估计$f(a_{0i})$，给定$L_0$，$A_0\sim N({\alpha }_0^{\ast },{{\sigma }^{\ast }}^2)$，因此$f(a_{0i}|l_{0i})$可表示为式(22)。
$$f(a_{0i}|l_{0i})=\frac{1}{\sqrt{2\pi {\hat{\sigma }}^2}}{e}^{-\frac{[a_{0i}-({\hat{\alpha }}_0+{\hat{\alpha }}_1{l}_{0i}){]}^2}{2{\hat{\sigma }}^2}}\text{\hspace{1em}(21)}$$
$$f(a_{0i}|l_{0i})=\frac{1}{\sqrt{2\pi {{\hat{\sigma }}^{\ast }}^2}}{e}^{-\frac{[a_{0i}-{\hat{\alpha }}_0^{\ast }{]}^2}{2{{\hat{\sigma }}^{\ast }}^2}}\text{\hspace{1em}(22)}$$

多时刻建模：时依性处理（Time-dependent Treatments）

  记上标$\bar{A}$为时序处理序列（历史剂量），$\bar{A}=(A_0,\dots ,A_k)$。其他变量不再赘述。在最简单的情况下，即每个处理$a_i$都是二值的，那么$Y_{\bar{a}}$有$2^k$种可能取值。因此，在这里假设仍然服从简化剂量反应关系，则线性逻辑MSM因果模型可以写作式(23)。
$$logitpr[Y_{\bar{a}}=1]={\beta }_0+{\beta }_{1}cum(\bar{a})\text{\hspace{1em}(23)}$$
  式中$cum(\bar{a})=\Sigma a_k$是历史处理的累计剂量和。相应的，数据模型可以写作式(24)。
$$logitpr[Y=1|\bar{A}=\bar{a}]={\beta }_0^{{}^{\prime }}+{\beta }_1^{{}^{\prime }}cum(\bar{a})\text{\hspace{1em}(24)}$$
  多时刻建模传统方法的不足：使用未调整权重的逻辑回归模型将不可避免地引入偏倚。这是因为⑴$L_k$是后续处理的混杂因子，因此必须调整；⑵但同时$L_k$是由前面处理影响的，因此不能被标准回归方法调整。因此使用$L_k$计算$s{w}_i$用于处理$A_k$的去混杂，而不是将$L_k$加入回归方程中，这也是本文的目的所在。

多时刻建模：稳定权重的求解

  在仅有$L_k$的混淆影响下，可以使用$s{w}_i$得到无偏估计，见式(25)。
$$s{w}_i=\frac{{\prod }_{k=0}^{K}pr[A_k=a_{ki}|{\bar{A}}_{k-1}={\bar{a}}_{(k-1)i}]}{{\prod }_{k=0}^{K}pr[A_k=a_{ki}|{\bar{A}}_{k-1}={\bar{a}}_{(k-1)i},{\bar{L}}_k={\bar{l}}_{ki}]}\text{\hspace{1em}(25)}$$
  对$s{w}_i$的求解分为两部分，第一部分建立回归模型，第二部分进行参数估计并代入回$s{w}_i$计算公式。

  对$s{w}_i$的分母和分子建立回归模型，需要考虑处理和混杂因子的实际物理意义。例如若处理的概率与日期$k$、前两天的处理、今天和昨天的混杂因子、昨天的处理和今天的混杂因子相互作用、和基线混杂因子（baseline covariates）有关，那么模型可以写作式(26)（这里同样假设处理是二值的）。
$$\begin{gathered} logitpr[A_k=1|{\bar{A}}_{k-1}={\bar{a}}_{k-1},{\bar{L}}_k={\bar{l}}_k]={\alpha }_0+{\alpha }_{1}k+{\alpha }_2{a}_{k-1}+{\alpha }_3{a}_{k-2} \\ +{\alpha }_4{l}_k+{\alpha }_5{l}_{k-1}+{\alpha }_6{a}_{k-1}{l}_k+\alpha +7l_0 \end{gathered}$$
$$logitpr[A_k=1|{\bar{A}}_{k-1}={\bar{a}}_{k-1}]={\alpha }_0^{\ast }+{\alpha }_1^{\ast }k+{\alpha }_2^{\ast }{a}_{k-1}+{\alpha }_3^{\ast }{a}_{k-2}\text{\hspace{1em}(26)}$$
  对于每类$i$，可以由式(25)求得$pr$的最大似然估计值${\hat{p}}_{0i},\dots ,{\hat{p}}_{Ki}$和${\hat{p}}_{1i}^{\ast },\dots ,{\hat{p}}_{Ki}^{\ast }$。当$A_k=0$时，显然估计值为$1-{\hat{p}}_{0i}$，不再赘述。因此，将估计值代入式(24)，可得$s{w}_i$的计算式(27)。
$$s{w}_i=\frac{{\prod }_{k=0}^K({\hat{p}}_{ki}^{\ast }{)}^{a_{ki}}(1-{\hat{p}}_{ki}^{\ast }{)}^{1-a_{ki}}}{{\prod }_{k=0}^K({\hat{p}}_{ki}{)}^{a_{ki}}(1-{\hat{p}}_{ki}{)}^{1-a_{ki}}}\text{\hspace{1em}(27)}$$

多时刻建模：预处理协变量带来的效应修饰作用（Effect Modifier）

  效应修饰作用（effect modifier）反映的是处理与结果关系的强弱，通常情况下仅与结果有关，与处理无关，不会引起偏倚（bias），将其加入回归模型能提升各类（subject）的估计准确性。而混杂因子（confounder）反应的是处理与结果关系的有无，与处理和结果均存在关联，不考虑混杂会产生偏倚。直觉性的图解见图(c {} )^[3]，其中$A$表示效应修饰作用，$C_1$表示混杂因子。

图(c) 效应修饰作用与混杂因子

  假设$V$是预处理协变量$L_0$的子集，由于$V$已被IPTW调整过，因此仅当$V$确定对因果效应有很大影响时，才会考虑将其加入回归方程。将$V$加入到式(24)，一个数据模型的例子如式(28)所示，相应的$s{w}_i$概率调整如式(29)所示。
$$logitpr[Y=1|\bar{A}=\bar{a},V=v]={\beta }_0+{\beta }_{1}cum(\bar{a})+{\beta }_{2}v+{\beta }_{3}cum(\bar{a})v\text{\hspace{1em}(28)}$$
$$logitpr[A_k=1|{\bar{A}}_{k-1}={\bar{a}}_{k-1},V=v]={\alpha }_0^{\ast }+{\alpha }_1^{\ast }k+{\alpha }_2^{\ast }{a}_{k-1}+{\alpha }_3^{\ast }{a}_{k-2}+{\alpha }_4^{\ast }v\text{\hspace{1em}(29)}$$

失访情况下的因果效应分析

  失访表示失去随访（lose to follow-up），通常是由于观察对象死亡、结束或数据无法收集等导致的数据缺失情况。使用MSM模型能在失访的情况下进行因果效应分析。记$k$时刻失访为$C_k=1$，未失访为$C_k=0$，并假设失访后目标不会再次接受随访。将失访加入模型中的思想也比较简单，是在$A_{k-1}$和$L_k$间插入$C_k$，即$A_{k-1}⟶C_k⟶L_k$。其物理意义例如，当病人服用一定剂量的药物后，可能失访，导致结果不可知。

  接下来是如何把$C_k$加入MSM回归模型的问题。结果$Y$能观测必然的前提是$\bar{C}=(C_0,\dots ,C_{K+1})=0$，也就是说在整个随访周期中均未发生失访现象。因此，各个类别的权重$s{w}_i$就需要加入对随访丢失的考虑，随访丢失越严重，类别所占比例越小，权重调整就越大。失访情况下的权重可写作$s{w}_i\times s{w}_i^{†}$。$s{w}_i$仍然使用式(25)，$s{w}_i^{†}$是逆审查权重概率（inverse probability of censoring weighting），写作式(30)。
$$s{w}_i^{†}=\frac{{\prod }_{k=0}^{K+1}pr[C_k=0|{\bar{C}}_{k-1}=0,{\bar{A}}_{k-1}={\bar{a}}_{(k-1)i}]}{{\prod }_{k=0}^{K+1}pr[C_k=0|{\bar{C}}_{k-1}=0,{\bar{A}}_{k-1}={\bar{a}}_{(k-1)i},{\bar{L}}_k={\bar{l}}_{ki}]}\text{\hspace{1em}(30)}$$

MSM的局限性

  当混杂因子在某个状态下，所有观察对象的处理都是相同的，这时无法使用MSM模型进行计算。（因为有的概率值为0，连乘也为0。这表明MSM服从正值性Positivity假设）

---

参考

[1] James M. Robins, Miguel A. Hernan, Babette Brumback. Marginal Structural Models and Causal Inference in Epidemiology[J]. Epidemiology, 2000, 11(5): 550-60.
[2] Thomas Fuchs. Slide 1 of MSM[DB/OL]. 2006. (Accessed in October 29, 2020)
[3] Tyler J. VanderWeele. Confounding and Effect Modification: Distribution and Measure[J]. Epidemiologic Methods, 2012, 1(1): 55-82.