逆概率加權(quán)構(gòu)建的邊際結(jié)構(gòu)模型在中介效應(yīng)分析中的應(yīng)用*

朱鉦宏 羅家俊 陳 雯 王 瓊△

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

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

原理和方法

1.傳統(tǒng)的中介效應(yīng)分析

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

圖1 傳統(tǒng)中介效應(yīng)分析的關(guān)聯(lián)模型

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

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

(1)

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

(2)

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

2.傳統(tǒng)中介效應(yīng)分析的局限性

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

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

圖2 不滿足中介效應(yīng)分析中混雜因素假設(shè)的兩種情景

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

3.逆概率加權(quán)構(gòu)建的邊際結(jié)構(gòu)模型

采用逆概率加權(quán)構(gòu)建邊際結(jié)構(gòu)模型,是基于反事實(shí)理論的因果推斷框架下一種中介分析的新方法。基本原理是利用逆概率加權(quán)構(gòu)建一個新的虛擬人群。在新的虛擬人群中暴露和中介相互獨(dú)立,且不再與其他混雜存在關(guān)聯(lián),所以此時的模型無需調(diào)整混雜C和L,且模型的因變量是反事實(shí)結(jié)局的邊緣均值,因而也被稱為邊際結(jié)構(gòu)模型[12]。

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

當(dāng)存在受到暴露影響的“中介-結(jié)果”混雜因素L時,采用逆概率加權(quán)構(gòu)建邊際結(jié)構(gòu)模型可分別將M看作獨(dú)立中介因素,也可將M和L看作聯(lián)合中介。

(1)M作為獨(dú)立中介構(gòu)建邊際結(jié)構(gòu)模型

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

特別地,當(dāng)暴露A為二分類變量時:

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

同理,當(dāng)中介M為二分類變量時:

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

圖3 逆概率加權(quán)構(gòu)建的虛擬人群中變量的關(guān)聯(lián)模型

E(Ma)=α0+αaa

(3)

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

(4)

(2)M和L作為聯(lián)合中介構(gòu)建邊際結(jié)構(gòu)模型

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

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

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

E(Ya) =β0+βaa

(5)

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

(6)

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

數(shù)據(jù)模擬和結(jié)果

在R中,本文基于正態(tài)分布產(chǎn)生一個10000行的隨機(jī)數(shù)據(jù)集,模擬了存在受暴露影響的“中介-結(jié)果”混雜因素時的數(shù)據(jù)結(jié)構(gòu)。具體模擬參數(shù)設(shè)置如下:

1.混雜C滿足均值為零,標(biāo)準(zhǔn)差為1的正態(tài)分布;

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

Log Odds=0.1-0.5C;

3.受暴露影響的“中介-結(jié)果”混雜因素L滿足均值為1+0.2A,標(biāo)準(zhǔn)差為1的正態(tài)分布;

4.中介M滿足均值為1+0.3A+0.5L+0.7C,標(biāo)準(zhǔn)差為1的正態(tài)分布;

5.結(jié)局Y滿足均值為1+0.4A+0.5M+0.5L+0.7C,標(biāo)準(zhǔn)差為1的正態(tài)分布。

圖4反映了模擬數(shù)據(jù)集中各變量之間的關(guān)聯(lián),本文通過簡單的計算可得到:當(dāng)M作為獨(dú)立中介時,總效應(yīng)為0.7,直接效應(yīng)為0.5,間接效應(yīng)為0.2;當(dāng)M和L作為聯(lián)合中介時,總效應(yīng)為0.7,直接效應(yīng)為0.4,間接效應(yīng)為0.3。

圖4 模擬隨機(jī)數(shù)據(jù)集的變量間關(guān)聯(lián)

基于該隨機(jī)數(shù)據(jù)集,本文采用了以下分析策略:

1.通過傳統(tǒng)的回歸模型進(jìn)行中介效應(yīng)分析,將L作為混雜進(jìn)行控制;

2.通過傳統(tǒng)的回歸模型進(jìn)行中介效應(yīng)分析,不對L進(jìn)行控制;

3.通過逆概率加權(quán)構(gòu)建的邊際結(jié)構(gòu)模型進(jìn)行中介效應(yīng)分析,將M作為獨(dú)立中介;

4.通過逆概率加權(quán)構(gòu)建的邊際結(jié)構(gòu)模型進(jìn)行中介效應(yīng)分析,將M和L作為聯(lián)合中介。

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

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

討 論

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

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

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

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

附錄：R代碼

#建立隨機(jī)數(shù)據(jù)集

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(連續(xù)性變量算法)

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)

#計算總權(quán)重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)

#自然直接效應(yīng)

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

#自然間接效應(yīng)

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

#使用自舉法估計間接效應(yīng)的標(biāo)準(zhǔn)誤差和95%置信區(qū)間

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"))

#使用德爾塔法估計間接效應(yīng)的標(biāo)準(zhǔn)誤差和95%置信區(qū)間

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(聯(lián)合中介算法)

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)

#計算總權(quán)重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)

#總效應(yīng)

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

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

#自然直接效應(yīng)

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

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

#自然間接效應(yīng)

NIE<- TE - NDE