閱讀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.