逆概率加權構建的邊際結構模型在中介效應分析中的應用*

朱鉦宏 羅家俊 陳 雯 王 瓊△

【提 要】 目的 基于傳統回歸模型的中介效應分析方法需要混雜因素滿足較強的假設,而對于一些復雜情景,如存在受到暴露影響的“中介-結果”混雜因素時,傳統的回歸方法可能不再適用。本文將介紹該情景下如何通過逆概率加權構建邊際結構模型進行中介效應分析。方法 詳細介紹存在受到暴露影響的“中介-結果”混雜因素時,傳統分析方法的局限性以及通過逆概率加權構建邊際結構模型進行替代分析的原理和方法。通過R模擬隨機數據集,分別采用傳統模型和逆概率加權構建邊際結構模型進行中介效應分析,并比較不同方法的評估結果。結果 存在受暴露影響的“中介-結果”混雜時,傳統模型效應的估計存在偏差。通過計算逆概率權重構建邊際結構模型進行中介效應分析,所估計的總效應、直接效應、間接效應與模擬的真實效應具有較好的一致性。結論 存在受暴露影響的“中介-結果”混雜的情況下,可以利用邊際結構模型有效地分析中介效應。

中介效應分析是探究風險因素和結局之間因果關系的重要方法之一[1],最初主要應用于心理學領域,成為許多經典和創新理論的基礎。過去十年間,中介效應分析在公共衛生領域也得到了廣泛的應用[2]。例如在流行病學研究中,研究者可以通過中介效應分析,評估危險因素和結局的總關聯中不同生物通路的貢獻,從而闡明危險因素影響健康的生物學機制[3]。近年來,隨著中介效應分析在多個領域的廣泛應用,其理論和方法也逐步得到發展和完善。本文將首先簡要介紹傳統的中介效應分析,進一步介紹當存在受到暴露影響的“中介-結果”混雜因素時,傳統的中介效應分析的局限,以及如何采用逆概率加權構建的邊際結構模型進行中介效應分析的思路和方法。最后,通過模擬隨機數據集對幾個方法的評估結果進行比較。

原理和方法

1.傳統的中介效應分析

中介效應分析中涉及的變量及關系如圖1(A)所示,自變量A如果通過某一變量M對因變量Y產生一定影響,則稱M為A和Y的中介變量;協變量C對自變量A和因變量Y都存在影響,但A并不會通過C而影響Y的發生,則稱C為影響A和Y之間關聯的混雜因素。

圖1 傳統中介效應分析的關聯模型

Baron和Kenny在1986年提出了經典的中介效應分析思路,即將總效應(TE)分解為直接效應(DE)和間接效應(IE)[4]。在這一思路的指導下,既往研究者主要基于圖1(B)所示的暴露、中介、結果和混雜因素間的關系,采用傳統回歸模型進行中介效應分析[5]。以暴露、中介和結果均為連續性變量為例:首先擬合中介M與暴露A的回歸模型(1),并控制協變量C;然后擬合結局Y與暴露A的回歸模型(2),并控制中介M和協變量C。其基本結構如下:

E(M|A,C)=α0+αaA+αcC

(1)

E(Y|A,M,C)=β0+βaA+βmM+βcC

(2)

由模型估計的系數可以得到直接效應以及間接效應的估計值:DE=βa;IE=βmαa。間接效應(IE)的標準誤差和95%置信區間推薦使用自舉法(bootstrap method)進行估計[6]。自舉法是在1個容量為n的原始樣本中重復抽取一系列容量也是n的隨機樣本,并保證每次抽樣中每一樣本觀察值被抽取的概率都是1/n,從而估計β的標準誤差并確定一定置信系數下β的置信區間[7]。此外,德爾塔法(delta method)也稱為誤差傳遞法,亦可用于計算間接效應的標準誤差和置信區間[8]。 需要注意的是,德爾塔法得到的標準誤差是一種近似解,其置信區間是對稱的;當使用德爾塔法得到的區間臨近統計顯著性邊界時,須謹慎考慮置信區間的意義。

2.傳統中介效應分析的局限性

隨著中介效應分析在多個領域的廣泛應用,研究者們在基于反事實理論的因果推斷框架下進一步發展了中介效應分析[9-10],將總效應分解為自然直接效應(NDE)和自然間接效應(NIE)[11],其中的“自然”強調了它是通過反事實進行定義的。在此基礎上,統計學家們提出中介效應分析結果用于進行因果解釋時,混雜因素必須滿足的假設條件:(1)沒有未測量的混雜因素(包括“暴露-結果”混雜、“暴露-中介”混雜和“中介-結果”混雜);(2)不存在受到暴露影響的“中介-結果”混雜因素。

當假設(1)不能滿足,即存在未測量的混雜因素U,如圖2(A)所示。此時,對未測量的混雜因素進行敏感性分析,可以觀察到暴露與結果之間的關聯是源于暴露對結果的真實影響,還是僅僅歸因于混雜因素。以“暴露-結果”混雜為例,敏感性分析的基本思想是指定未測量的混雜因素和暴露之間(即U-A),以及未測量的混雜因素和結果之間(即U-Y)的關聯系數來構建變量U。在此基礎上構建模型,控制協變量和未測量的混雜因素U以后獲得“校正”后的效應估計。這些自定義的關聯系數本身沒有實際的意義,但研究者可以通過在合理的范圍內取一系列不同系數獲得對應的一系列“校正”效應,從而獲得因果效應的合理估計范圍[10]。

圖2 不滿足中介效應分析中混雜因素假設的兩種情景

當假設(2)不能滿足,即存在受到暴露影響的“中介-結果”混雜因素L,如圖2(B)所示。此時情況變得更為復雜,因為L既是一個“中介-結果”混雜因素,同時又處于從暴露到結果的因果路徑上,傳統的回歸方法此時受到了挑戰。假設將L作為上述介紹模型(2)中的協變量進行控制,可能會阻斷通過L的因果路徑,即 A-L-Y,對A-Y的總效應估計會產生偏差。另一種情況,如果不在回歸模型中調整L,那么對M-Y的效應估計又會存在偏差,因為L是M-Y的混雜因素。因此,無論是否在回歸模型中控制L,都會得到有偏的估計,傳統回歸方法不再適用。此時,通過逆概率加權構建的邊際結構模型可以較好地解決這個問題,也是本文將主要介紹的方法。

3.逆概率加權構建的邊際結構模型

采用逆概率加權構建邊際結構模型,是基于反事實理論的因果推斷框架下一種中介分析的新方法?；驹硎抢媚娓怕始訖鄻嫿ㄒ粋€新的虛擬人群。在新的虛擬人群中暴露和中介相互獨立,且不再與其他混雜存在關聯,所以此時的模型無需調整混雜C和L,且模型的因變量是反事實結局的邊緣均值,因而也被稱為邊際結構模型[12]。

當暴露A為二分類變量時,每個人暴露A的逆概率權重被定義為1/P(A|C),每個人受到暴露的條件概率P(A|C)也被稱為傾向性評分,因此逆概率權重也可以認為是傾向性評分的倒數。當暴露為連續性變量時,我們需要計算穩定的逆概率權重P(A)/P(A|C)。穩定權重的分子常用原始樣本中暴露的概率P(A),亦被稱為穩定因子,它可以有效縮小逆概率權重的取值范圍。穩定權重的均值是1,這是因為此時虛擬人群和原樣本的人數一樣,因此需要檢驗穩定權重的均值是否為 1,如果偏離 1,則表明估計權重的模型設定可能有誤,需要檢驗暴露變量的分布是否滿足正態分布。

當存在受到暴露影響的“中介-結果”混雜因素L時,采用逆概率加權構建邊際結構模型可分別將M看作獨立中介因素,也可將M和L看作聯合中介。

(1)M作為獨立中介構建邊際結構模型

第一步:計算暴露的逆概率權重(inverse-probability-of-treatment weight,IPTW)。分子為暴露A的估計邊際概率,分母為協變量C條件下每個樣本的實際暴露A的條件概率。IPTW的作用是移除C-A之間的關聯,即使得暴露A不再受到混雜因素C的影響。

特別地,當暴露A為二分類變量時:

第二步:計算中介的逆概率權重(inverse-probability-of-mediator weight,IPMW),分子為中介M在暴露A為條件下的估計概率,分母為中介M在以暴露A,協變量C和受暴露影響的“中介-結果”混雜L為條件下的估計概率。IPMW的作用是移除C-M、A-M以及L-M之間的關聯,即使得中介變量M既不受到暴露A影響,也不受到混雜C和L影響。然后通過IPMW乘以IPTW獲得總權重W。

同理,當中介M為二分類變量時:

第三步:構建逆概率加權的邊際結構模型。擬合IPTW加權的邊際結構模型(3);擬合W加權的邊際結構模型(4)。虛擬人群的變量間關系如圖3(A)所示,由于IPTW和IPMW消除了混雜的影響,此時,構建邊際結構模型可直接估計自然直接效應和自然間接效應:NDE=βa;NIE=βmαa。在邊際結構模型中,反事實變量Ma表示暴露A取值為a時,觀察到的中介變量M;反事實變量Ya,m表示暴露A取值為a以及中介M取值為m時,觀察到的結局變量Y。

圖3 逆概率加權構建的虛擬人群中變量的關聯模型

E(Ma)=α0+αaa

(3)

E(Ya,m) =β0+βaa+βmm

(4)

(2)M和L作為聯合中介構建邊際結構模型

將M和L看作一個整體,即聯合中介(M,L),那么假設(2)將再次得到滿足[13]。在這樣的處理思路下,也可以通過逆概率加權構建的邊際結構模型來進行分析。

把M和L看作聯合中介因素時,需要在第二步的基礎上重新擬合一個IPMW,這個新的IPMW由兩部分組成,分別對應兩個中介變量M和L:①分子為中介M在暴露A、中介L為條件下的估計概率,分母為中介M在以暴露A、中介L、協變量C為條件下的估計概率;②分子為中介L在暴露A為條件下的估計概率,分母為中介L在以暴露A、協變量C為條件下的估計概率。

此時,通過逆概率加權,構建了一個暴露A不受到混雜因素C影響,中介變量M和L既不受到暴露A影響也不受到混雜因素C影響的虛擬人群,如圖3(B)所示。

E(Ya) =β0+βaa

(5)

E(Ya,m,l) =β′0+β′aa+β′mm+β′ll

(6)

然后使用IPTW加權模型(5)可估計暴露影響結果的總效應:TE=NDE+NIE=βa;W加權的邊際結構模型(6)可估計暴露影響結果的自然直接效應:NDE=β′a;采用差分法[2]從模型總效應中減去估計的自然直接效應,從而得到自然間接效應:NIE=βa-β′a。

數據模擬和結果

在R中,本文基于正態分布產生一個10000行的隨機數據集,模擬了存在受暴露影響的“中介-結果”混雜因素時的數據結構。具體模擬參數設置如下:

1.混雜C滿足均值為零,標準差為1的正態分布;

2.暴露A為二分類變量,與C的關系滿足:

Log Odds=0.1-0.5C;

3.受暴露影響的“中介-結果”混雜因素L滿足均值為1+0.2A,標準差為1的正態分布;

4.中介M滿足均值為1+0.3A+0.5L+0.7C,標準差為1的正態分布;

5.結局Y滿足均值為1+0.4A+0.5M+0.5L+0.7C,標準差為1的正態分布。

圖4反映了模擬數據集中各變量之間的關聯,本文通過簡單的計算可得到:當M作為獨立中介時,總效應為0.7,直接效應為0.5,間接效應為0.2;當M和L作為聯合中介時,總效應為0.7,直接效應為0.4,間接效應為0.3。

圖4 模擬隨機數據集的變量間關聯

基于該隨機數據集,本文采用了以下分析策略:

1.通過傳統的回歸模型進行中介效應分析,將L作為混雜進行控制;

2.通過傳統的回歸模型進行中介效應分析,不對L進行控制;

3.通過逆概率加權構建的邊際結構模型進行中介效應分析,將M作為獨立中介;

4.通過逆概率加權構建的邊際結構模型進行中介效應分析,將M和L作為聯合中介。

不同建模方式下的估計結果對比如表1所示,當存在受到暴露影響的“中介-結果”混雜因素時,如果仍然使用傳統回歸方法,若控制變量L,將低估直接效應和總效應,而不控制變量L,將低估直接效應并高估間接效應。采用逆概率加權構建邊際結果模型,不管是以M作為唯一中介,還是將M和L作為聯合中介,模擬真實值均處于估計值的95%置信區間內,表明逆概率加權構建的邊際結構模型有效地估計了中介效應。

表1 不同建模方式下的估計結果對比

討 論

逆概率加權構建的邊際結構模型在中介效應分析中已經得到了成熟的應用,如Kaisla Komulainen等人將其用于兒童社會心理環境與成年期心臟健康的分析中[14];Yongfu Yu等人采用該方法分析了早產和小于胎齡兒在產婦教育與嬰兒死亡率關聯中的中介作用[15]。本文系統介紹了當存在受暴露影響的“中介-結果”混雜時,如何通過逆概率加權技術構建邊際結構模型,進而實現中介效應分析。并通過在R中模擬隨機數據集,詳細介紹了計算逆概率權重構建邊際結構模型的步驟和代碼。本文發現,當存在受暴露影響的“中介-結果”混雜時,采用傳統回歸模型估計效應存在偏差,而利用逆概率權重構建邊際結構模型估計的總效應、直接效應、間接效應與模擬的真實效應具有較好的一致性。

但需要指出的是,逆概率加權構建的邊際結構模型本身也存在一定的局限性。如果暴露和中介是連續性變量,此時逆概率加權通常不夠穩定。對連續性變量構建逆概率權重首先需要對連續變量的分布進行假設。當變量是二分變量的時候,“逆概率”中的“概率”就是預測變量取值為1(或0)時的概率;當變量是連續變量的時候,此時的“概率”應該是變量不同觀測值下的概率密度。要得到準確的概率密度,就需要先驗地假設變量分布。通常情況下,多數研究者會假設連續變量服從正態分布。但是在實踐中,很難保證所有變量都服從正態分布。此外,如果暴露或者中介是連續性變量,逆概率權重的分母可能在樣本之間變化很大,這種可變性將導致少數樣本的權重值非常大,進而會為虛擬人群貢獻大量的自身副本,從而主導加權分析。

解決連續性變量的逆概率權重有以下幾種可能的途徑:(1)對逆概率權重本身做一些處理,如剪切[16]或截斷[17]。剪切的原理是將極端傾向性評分或權重的個體直接刪除,截斷是將預先設定范圍之外的傾向性評分或權重,統一賦為范圍臨界點的值。在數據分析時,可將穩定權重與截斷法結合使用,對比取不同截斷點對應的效應值,從而獲得穩健的估計。(2)使用更加穩健的權重,如協變量平衡的傾向性評分[18]或者雙重穩健估計的逆概率權重[19]。(3)如果研究者對直接效應更感興趣,可以使用結構均值模型來估計連續暴露的受控直接效應[10]。此外,即使不存在受暴露影響的“中介-結果”混雜,逆概率加權構建的邊際結構模型依舊可以使用,作為傳統回歸方法的替代方法。然而,當不存在受暴露影響“中介-結果”混雜時,傳統回歸方法更有效(即有較小的標準誤差)。因此,除非受暴露影響的“中介-結果”混雜確實存在,否則傾向于使用傳統的基于回歸的方法。

中介效應分析在流行病學研究中已經得到了快速的發展,但其所依賴的強假設使其在研究具體問題時遇到了許多挑戰。與此同時,統計學家們不斷地開發出了新的中介分析方法來適應不同的復雜情景,除本文介紹的基于回歸的經典方法和逆概率加權構建的邊際結構模型之外,還有VanderWeele提出的基于加權的方法[13],Tchetgen Tchetgen提出的逆比值比加權方法[20],Vansteelandt等人提出的自然效應模型[21]和Robins提出的g-formula方法[22]等。事實上,并沒有一種完美的中介效應分析方法可以徹底解決所有問題,根據研究目的和感興趣的中介變量,研究者可以靈活地選擇適合的假設框架與分析方法,識別出感興趣的效應并實現對因果機制路徑的探索。

附錄：R代碼

#建立隨機數據集

n<- 10^4

C<- rnorm(n,mean=0,sd=1)

A<- rbinom(n,size=1,prob=plogis(0.1 - 0.5*C))

L<- rnorm(n,mean=1 + 0.2*A,sd=1)

M<- rnorm(n,mean=1 + 0.3*A + 0.5*L + 0.7*C,sd=1)

Y<- rnorm(n,mean=1 + 0.4*A + 0.5*M + 0.5*L + 0.7*C,sd=1)

dat<- data.frame(id=1 :n,C,A,M,L,Y)

head(dat)

#第一步:計算IPTW(二分類變量算法)

mod<- glm(A ～ C,family=binomial(link="logit"),data=dat)

dat$ps<- predict(mod,dat,type="response")

summary(dat$ps)

dat$iptw<- ifelse(dat$A==1,1/dat$ps,1/(1 - dat$ps))

#第二步:計算IPMW(連續性變量算法)

model_num<- lm(M ～ as.factor(A),data=dat)

num<- dnorm(dat$M,predict(model_num),sd(model_num$residuals))

model_den<- lm(M ～ as.factor(A)+ L + C,data=dat)

den<- dnorm(dat$M,predict(model_den),sd(model_den$residuals))

dat$ipmw<- num/den

summary(dat$ipmw)

#計算總權重W

dat$w<- dat$iptw*dat$ipmw

#第三步:建模(3)-(4)

mod3<- lm(M ～ as.factor(A),weights=iptw,data=dat)

mod4<- lm(Y ～ as.factor(A)+ M,weights=w,data=dat)

#自然直接效應

confint(mod4,"as.factor(A)1",level=0.95)

#自然間接效應

NIE<- coef(summary(mod3))[2,1]*coef(summary(mod4))[3,1]

#使用自舉法估計間接效應的標準誤差和95%置信區間

library(boot)

fc<- function(dat,i){

newdat<- dat[i,]

mod3<- lm(M ～ as.factor(A),weights=iptw,data=newdat)

mod4<- lm(Y ～ as.factor(A)+ M,weights=w,data=newdat)

out<- coef(summary(mod3))[2,1]*coef(summary(mod4))[3,1]

return(out)

}

set.seed(123)

bootm=boot(dat,fc,R=1000)

boot.ci(boot.out=bootm,type=c("basic"))

#使用德爾塔法估計間接效應的標準誤差和95%置信區間

library(msm)

estmean<- c(coef(summary(mod3))[2,1],coef(summary(mod4))[3,1])

estvar<- matrix(c(vcov(mod3)[2,2],0,0,vcov(mod4)[3,3]),2,2)

beta_se<- deltamethod(～ x1 *x2,estmean,estvar)

beta_lwr<- NIE - qnorm(0.975)*beta_se

beta_upr<- NIE + qnorm(0.975)*beta_se

round(c(beta_lwr,beta_upr),2)

#計算IPMW(聯合中介算法)

model_num1<- lm(M ～ as.factor(A)+ L,data=dat)

num1<- dnorm(dat$M,predict(model_num1),sd(model_num1$residuals))

model_den1<- lm(M ～ as.factor(A)+ L + C,data=dat)

den1<- dnorm(dat$M,predict(model_den1),sd(model_den1$residuals))

model_num2<- lm(L ～ as.factor(A),data=dat)

num2<- dnorm(dat$M,predict(model_num2),sd(model_num2$residuals))

model_den2<- lm(L ～ as.factor(A)+ C,data=dat)

den2<- dnorm(dat$M,predict(model_den2),sd(model_den2$residuals))

dat$ipmw<-(num1/den1)*(num2/den2)

summary(dat$ipmw)

#計算總權重W

dat$w<- dat$iptw*dat$ipmw

#建模(5)-(6)

mod5<- lm(Y ～ as.factor(A),weights=iptw,data=dat)

mod6<- lm(Y ～ as.factor(A)+ M + L,weights=w,data=dat)

#總效應

TE<- coef(summary(mod5))[2,1]

confint(mod5,"as.factor(A)1",level=0.95)

#自然直接效應

NDE<- coef(summary(mod6))[2,1]

confint(mod6,"as.factor(A)1",level=0.95)

#自然間接效應

NIE<- TE - NDE