南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)系(211166)

施倩雯 魏永越 李清雅 段巍巍 趙 楊 陳 峰△

?

疾病預(yù)后研究的中介分析方法評價*

南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)系(211166)

施倩雯 魏永越 李清雅 段巍巍 趙 楊 陳 峰△

目的 對用于癌癥預(yù)后研究的五種中介分析方法(VanderWeele法、Baron-Kenny法、Imai法、Sobel法和InverseWeight法)進(jìn)行評價，為實例分析的方法選擇提供依據(jù)。方法 基于模擬試驗，產(chǎn)生不同參數(shù)設(shè)置下的模擬數(shù)據(jù)，并評價五種方法的第一類錯誤、檢驗效能和分析時間。結(jié)果 除Inverse Weight法在相關(guān)系數(shù)較大時第一類錯誤有所膨脹外，其余四種方法的第一類錯誤在不同參數(shù)情況下均在0.05附近。五種方法的檢驗效能趨勢一致，均隨著樣本量、中介比、總效應(yīng)的增大而增大，隨著刪失比的增大而減小。在樣本量較小(N=100)且中介比不大于30%的情況下，Inverse Weight法的檢驗效能低于另四種方法。Inverse Weight法、Baron-Kenny法和Imai法的分析效率遠(yuǎn)低于VanderWeele法和Sobel法。結(jié)論 綜合考慮一類錯誤控制、檢驗效能及分析效率，推薦VanderWeele法進(jìn)行預(yù)后研究的中介分析。

生存資料 中介分析 VanderWeele法

近年來，因環(huán)境污染加劇、生活方式改變，以及醫(yī)療水平提高等多方面原因，越來越多的癌癥病人被確診，癌癥的治療和預(yù)后也愈發(fā)引人關(guān)注。根據(jù)2015年的數(shù)據(jù)顯示，中國男性發(fā)病率最高的癌癥前三位依次為肺癌、胃癌和肝癌，女性則為乳腺癌、肺癌和結(jié)直腸癌[1,2]。以肺癌為例，查閱美國國立癌癥機(jī)構(gòu)統(tǒng)計的數(shù)據(jù)可知，其五年相對生存率只有17.7%[3]。影響其預(yù)后的因素有很多，包括吸煙、某些SNP位點等[4-7]，但其作用機(jī)制尚不清楚。而中介分析作為一種能夠探索暴露對結(jié)局的內(nèi)在作用機(jī)制的因果分析方法，受到了研究者的廣泛關(guān)注[8-9]。

中介分析雛形最早在1986年由Baron & Kenny提出(以下簡稱Baron-Kenny法)。1987年，Sobel基于Baron-Kenny法提出了中介效應(yīng)的方差估計方法(以下簡稱Sobel法)，計算簡單、原理易于理解，但是較為保守[10]。近年來，VanderWeele(VanderWeele法)、Imai(Imai法)和Tchetgen(Inverse Weight法)等人進(jìn)一步拓寬了中介分析的理論基礎(chǔ)及應(yīng)用范疇[11-13]。然而，至今尚無報道對以上五種方法進(jìn)行比較評價。

本研究將以預(yù)后的生存資料為切入點，建立一個暴露因素和一個中介因素的中介模型，通過統(tǒng)計模擬試驗，從第一類錯誤及檢驗效能角度對上述五種方法進(jìn)行統(tǒng)計學(xué)性質(zhì)評價，為預(yù)后的中介分析方法選擇提供建議。

原理與方法

1.Baron-Kenny法

Baron-Kenny法的原理是先建立中介因素(M)、暴露因素(A)和結(jié)局(Y)之間的模型，然后通過模型中的系數(shù)計算間接效應(yīng)(indirect effect，IE)，涉及的模型包括Y模型、M模型與總模型(公式1-3)[14]。以加速失效時間(acceleratefailure time，AFT)模型建立生存函數(shù)模型(即Y模型)，各模型如下：

Y模型：

(1)

M模型：

(2)

總模型：

log(OS)=α0+α1A+α2c+γε

(3)

據(jù)此，間接效應(yīng)的點估計為：

IE=eθ1×β2

(4)

直接效應(yīng)的點估計為：

IE=eβ1

(5)

間接效應(yīng)表示暴露通過中介因素介導(dǎo)對結(jié)局產(chǎn)生的作用大小。中介效應(yīng)占總效應(yīng)的比例(即，中介比)表示為：

(6)

Baron-Kenny法中，間接效應(yīng)的可信區(qū)間通過Bootstrap得到(抽樣1000次)。

2.Sobel法

模型假設(shè)和IE的估計方法與Baron-Kenny法一致，log(IE)的方差可表示為：

(7)

其中Sθ1和Sβ2分別表示θ1和β2的標(biāo)準(zhǔn)誤[10]。

3.VanderWeele法

VanderWeele法的原理是對總效應(yīng)進(jìn)行分解[15]，雖然同樣用到Y(jié)模型和M模型，但相較于Baron-Kenny法，VanderWeele法考慮了暴露與中介因素之間的交互作用，故Y模型可表示為：

(8)

VanderWeele法的M模型表示為：

(9)

通過公式(8)和公式(9)的推導(dǎo)可知間接效應(yīng)：IEAFT=e(θ1β2+θ1β3a)(a-a*)，直接效應(yīng)：

4.Inverse Weight法

中介因素M對暴露A的回歸模型為：

(10)

調(diào)整過權(quán)重的總模型：

(11)

5.Imai法

Imai法基于反事實原理，將平均總效應(yīng)(averagetotaleffect)分成平均因果中介效應(yīng)(ACMEs,averagecausalmediationeffects)和平均直接效應(yīng)(ADE,averagedirecteffect)兩部分[15,17-18]。

平均因果中介效應(yīng)：

(12)

平均直接效應(yīng)：

(13)

平均總效應(yīng)(averagetotaleffect)：

(14)

其中，a代表暴露水平，Mi(1)表示暴露存在時中介因素的水平，Mi(0)表示暴露不存在時中介因素的水平，Mi(a)表示暴露處于a水平時中介因素的水平。

模擬試驗

1.模擬數(shù)據(jù)的產(chǎn)生

基于如圖1的模型，按照以下步驟產(chǎn)生模擬數(shù)據(jù)：

(1) 指定總效應(yīng)(HR)，根據(jù)中介效應(yīng)比例(M%)把總效應(yīng)(HR)分解成直接效應(yīng)和間接效應(yīng)兩部分，分別可表示為：HR×(1-M%)和HR×M%，因DE=eβ1故可求得β1；

(2) 指定樣本量(N)、A與M間的相關(guān)系數(shù)r，可產(chǎn)生二元標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)A和M；此時r即為公式(2)中的θ1值，因IE=eθ1×β2故可得到β2；

(3) 把變量A、M、系數(shù)β1和β2帶入指數(shù)回歸模型的風(fēng)險函數(shù)得到參數(shù)λ，進(jìn)而基于指數(shù)分布隨機(jī)產(chǎn)生生存時間；

(4) 根據(jù)樣本量(N)和刪失比例(cen%)隨機(jī)產(chǎn)生截尾，1代表死亡，0代表刪失。

參考若干流行病學(xué)研究，總效應(yīng)(HR)設(shè)為1.5、2.0、2.5、3.0[19-24]。由于生存資料的刪失比越高，模型預(yù)測越不準(zhǔn)確[25]，故保險起見，模擬試驗時刪失比取值不超過50%。理論上中介效應(yīng)比例(M%)越大，中介分析的檢驗效能越高，而在實際分析中更希望能檢出中介效應(yīng)小的中介因素，不放過任何可能存在的因果關(guān)聯(lián)，故中介比(M%)設(shè)置在50%及以下?？紤]到預(yù)后研究的實際情況，模擬試驗中樣本量取100到1000。參數(shù)設(shè)置如表1。

表1 模擬試驗參數(shù)設(shè)置情況

2.模擬試驗流程

(1)根據(jù)參數(shù)設(shè)置產(chǎn)生模擬數(shù)據(jù)；

(2)用五種方法對產(chǎn)生的同一模擬數(shù)據(jù)分別進(jìn)行分析，提取相關(guān)結(jié)果，包括：間接效應(yīng)的點估計、區(qū)間估計(95%CI)、假設(shè)檢驗P值及檢驗結(jié)果(是否拒絕H0)；

(3)在不同參數(shù)情況下重復(fù)(1)～(2)步驟1000次，對1000次的結(jié)果進(jìn)行匯總，得到五種方法的第一類錯誤和檢驗效能。

結(jié) 果

1.一類錯誤

當(dāng)總效應(yīng)不存在(HR=1)的時候，IE的期望為0。如圖1為該情況下不同方法的第一類錯誤。橫坐標(biāo)為暴露A與中介M的相關(guān)系數(shù)r，縱坐標(biāo)為第一類錯誤。當(dāng)相關(guān)系數(shù)較低情況下，五種方法的第一類錯誤均接近于0，且隨著相關(guān)系數(shù)的增大而逐漸增大。樣本量較小情況下更為明顯，樣本量較大情況下，該趨勢有所緩和。相關(guān)系數(shù)較大的情況下，IW法的第一類錯誤略有膨脹。刪失比對第一類錯誤的影響不大。

圖1 總效應(yīng)不存在時暴露-中介的相關(guān)系數(shù)與第一類錯誤的關(guān)系

同理，當(dāng)暴露與中介因素不存在相關(guān)性的時候(r=0)，IE的期望亦為0。圖2為該情況下不同方法的第一類錯誤。橫坐標(biāo)為總效應(yīng)，縱坐標(biāo)為第一類錯誤。如圖所示，在樣本量較小(N=100)的情況下，五種方法的第一類錯誤隨著總效應(yīng)的增大逐漸增大；在樣本量較大(N=500,1000)的情況下，該趨勢不再明顯，五種方法的第一類錯誤均保持在0.05附近。刪失比對五種方法的第一類錯誤影響不明顯。

2.檢驗效能

如圖3所示(r=0.3，HR=1.5)，隨著樣本量的增加五種方法的檢驗效能均逐漸增大；IW法的檢驗效能明顯低于其余四種方法。同樣，如圖4所示(r=0.3，N=100)，隨著總效應(yīng)的增加，檢驗效能亦逐漸增大；IW法亦明顯低于其余四種方法。

當(dāng)考察相關(guān)系數(shù)r不同水平對間接效應(yīng)檢驗效能的影響時，發(fā)現(xiàn)隨著相關(guān)系數(shù)的增大，檢驗效能先增大后減小(如圖5)。當(dāng)中介比為10%時，由于檢驗效能過低，無法反映出方法間的差異；而當(dāng)中介比為30%及50%時，r在0.5以上時，檢驗效能逐漸下降，其中，中介比為30%、r在0.7以上時，檢驗效能較之峰值下降超過一半。

圖2 暴露-中介關(guān)系為零時總效應(yīng)與第一類錯誤的關(guān)系

圖3 r=0.3，HR=1.5時檢驗效能與樣本量的關(guān)系

其余參數(shù)情況下，趨勢類似。

根據(jù)間接效應(yīng)公式IE=eθ1×β2=er×β2可知，隨著相關(guān)系數(shù)r的增大，間接效應(yīng)理論上呈指數(shù)增加(在log尺度上呈線性增加)，其檢驗效能也應(yīng)升高。但模擬結(jié)果與此推測相違背。結(jié)合模型(1)，提出假設(shè)：隨著A和M相關(guān)系數(shù)的增大，模型(1)中A及M的系數(shù)估計的隨機(jī)擾動增加，導(dǎo)致系數(shù)β2的變異變大，從而間接效應(yīng)的變異變大，最終導(dǎo)致在較高相關(guān)系數(shù)的情況下檢驗效能減小。據(jù)此，設(shè)置模擬試驗參數(shù)(β1=β2=0.5，samplesize=100，cen%=0.2，r=0.01(0.01)0.99，基于VanderWeele法模擬5000次，發(fā)現(xiàn)：隨著相關(guān)系數(shù)的增大，β2的點估計波動變大，其方差呈現(xiàn)指數(shù)級膨脹；log(IE)的估計值呈線性上升，其方差亦呈指數(shù)級膨脹(圖6)。該模擬試驗結(jié)果驗證了前文的假設(shè)。由此得出，當(dāng)A與M相關(guān)性較高時(如，r>0.5)，五種中介分析方法皆存在檢驗效能較低的缺陷。

圖4 N=100，r=0.3時檢驗效能與總效應(yīng)的關(guān)系

圖5 N=100，HR=2時檢驗效能與相關(guān)系數(shù)的關(guān)系

3.分析效率

理論上，由于Baron-Kenny法、IW法和Imai法皆采用Bootstrap進(jìn)行區(qū)間估計，故分析效率較低。以樣本量為1000的生存數(shù)據(jù)為例，用VanderWeele法分析僅需1.1073×10-2秒，用Sobel法需1.2705×10-2秒；如Bootstrap次數(shù)設(shè)置為1000次，則用Baron-Kenny法、IW法和Imai法分析分別所需時間為：12.2784秒、21.3041秒和9.2768秒，分析效率較低。故當(dāng)樣本量較大且維度較高的情況下，不推薦使用Baron-Kenny法、IW法和Imai法。

圖6 相關(guān)系數(shù)與中介分析檢驗效能關(guān)系探討結(jié)果

討 論

除IW法，其余四種方法均能較好地控制第一類錯誤。五種方法的檢驗效能趨勢一致，均隨著樣本量、中介比、總效應(yīng)的增大而增大，隨著刪失比的增大而減小。總體上來看，IW法較其它方法檢驗效能較低，Baron-Kenny法和Imai法的分析效率遠(yuǎn)低于VanderWeele法和Sobel法，而Sobel法由于要求β2×θ1服從正態(tài)分布故在實際應(yīng)用中存在較大的局限性[26]。

本研究亦提示，當(dāng)X-M相關(guān)性較高(如r>0.5)情況下，各種方法皆存在檢驗效能下降的問題。殘差法及主成分分析等方法可為后續(xù)研究提供思路[27-28]。在實例分析中，本研究以簡單中介模型為基礎(chǔ)進(jìn)行模型試驗，而在實際分析中，可能存在中介因素同樣是生存數(shù)據(jù)的情況，往往還存在“多個暴露”、“多個中介因素”，或者“暴露與中介因素存在交互作用”等情況，在這些情況下該如何選擇中介分析方法，仍需進(jìn)一步探討。

Baron-Kenny法和Sobel法雖然模型簡單、原理便于理解，但這兩種方法存在一定的局限性——不能處理暴露與中介因素存在交互作用以及暴露非線性的情況，而VanderWeele法則能較好地處理這些問題，故應(yīng)用范圍更廣[13]。IW法對間接效應(yīng)的點估計是通過總效應(yīng)減去直接效應(yīng)間接得到的，同樣不適用于“暴露與中介因素存在交互作用”的情況[12]。故當(dāng)“暴露與中介因素存在交互作用”的情況下，優(yōu)先考慮VanderWeele法和Imai法，此時這兩種方法的統(tǒng)計性質(zhì)將是下一步的研究重點。

綜上所述，VanderWeele法能較好地控制第一類錯誤，檢驗效能亦可，且分析效率具備一定的優(yōu)勢，故推薦作為生存數(shù)據(jù)的中介分析方法，尤其是維度比較高的情況下。當(dāng)數(shù)據(jù)維度較低的情況下，Imai法亦是個不錯的選擇。

[1]Chen WQ,Zheng RS,Zhang SW,et al.Cancer statistics in China,2015.CA Cancer J Clin,2016,66(2):115-132.

[2]曾倩,崔芳芳,宇傳華，等.中國癌癥發(fā)病、死亡現(xiàn)狀與趨勢分析.中國衛(wèi)生統(tǒng)計.2016,33(2):321-323.

[3]SEER Stat Fact Sheets:Lung and Bronchus Cancer.Available from:http://seer.cancer.gov/statfacts/html/lungb.html.

[4]Jiang LP,Zhu ZT,He CY.Effects of CYP3A5 genetic polymorphism and smoking on the prognosis of non-small-cell lung cancer.Onco Targets Ther,2016,9:1461-1469.

[5]Dong L,Wei LX,Xu BH,et al.Association of GWAS-identified lung cancer susceptibility loci with survival length in patients with small-cell lung cancer treated with platinum-based chemotherapy.PLoS One,2014,9(11):e113574.

[6]Meguid RA,Hooker CM,Harris J,et al.Long-term survival outcomes by smoking status in surgical and nonsurgical patients with non-small cell lung cancer:comparing never smokers and current smokers.Chest.2010,138(3):500-509.

[7]Sasaki H,Tatemaysu T,Okuda K,et al.PD-1 gene promoter polymorphisms correlate with a poor prognosis in non-small cell lung cancer.Mol Clin Oncol,2014,2(6):1035-1042.

[8]Greenland S,Robins J.Invited commentary:ecologic studies-biases,misconceptions,and counterexamples.Am J Epidemiol.1994,139(8):747-760.

[9]Jo B.Causal inference in randomized experiments with mediational processes.Psychol Methods.2008,13(4):314-336.

[10]Sobel.Direct and indirect effects in linear structural equation models.Sociological Methods Research.1987,16(1):155-176.

[11]Imai K,Keele L,Tingley D.A general approach to causal mediation analysis.Psychol Methods,2010,15(4):309-334.

[12]Tchetgen Tchetgen EJ.Inverse odds ratio-weighted estimation for causal mediation analysis.Stat Med.2013,32(26):4567-4580.

[13]VanderWeele TJ.Explanation in causal inference:methods for mediation and interaction.2015,New York:Oxford University Press.706 pages.

[14]Baron RM,Kenny DA.The moderator-mediator variable distinction in social psychological research:conceptual,strategic,and statistical considerations.J Pers Soc Psychol.1986,51(6):1173-1182.

[15]VanderWeele TJ.A unification of mediation and interaction:a 4-way decomposition.Epidemiology.2014,25(5):749-761.

[16]Tchetgen Tchetgen EJ.On causal mediation analysis with a survival outcome.Int J Biostat.2011,7(1):Article 33.

[17]VanderWeele TJ.A three-way decomposition of a total effect into direct,indirect,and interactive effects.Epidemiology.2013,24(2):224-232.

[18]Bohnke JR.Explanation in causal inference:Methods for mediation and interaction.Q J Exp Psychol(Hove).2016,69(6):1243-1244.

[19]Osaki Y,Okamote M,Kaetsu A,et al.Retrospective cohort study of smoking and lung cancer incidence in rural prefecture,Japan.Environ Health Prev Med.2007,12(4):178-182.

[20]Wang A,Kubo J,Luo J,et al.Active and passive smoking in relation to lung cancer incidence in the Women′s Health Initiative Observational Study prospective cohort.Ann Oncol,2015,26(1):221-230.

[21]Takahashi K,Saito S,Kamamura Y,et al.Prognostic value of CD4+ lymphocytes in pleural cavity of patients with non-small cell lung cancer.Thorax,2001,56(8):639-642.

[22]Taguchi O,Gabazza EC,Yasui H,et al.Prognostic significance of plasma D-dimer levels in patients with lung cancer.Thorax.1997,52(6):563-565.

[23]Huovinen E,Kaprio J,Vesterinen E,et al.Mortality of adults with asthma:a prospective cohort study.Thorax,1997,52(1):49-54.

[24]Qiu F,Yang L,Lu X,et al.The MKK7 p.Glu116Lys Rare Variant Serves as a Predictor for Lung Cancer Risk and Prognosis in Chinese.Plos Genetics,2016,12(3).

[25]Watt DC,Aitchison TC,MacKie RM,et al.Survival analysis:the importance of censored observations.Melanoma Res,1996,6(5):379-385.

[26]Mackinnon DP,Lockwood CM,Williams J.Confidence Limits for the Indirect Effect:Distribution of the Product and Resampling Methods.Multivariate Behav Res,2004,39(1):99.

[27]Schroeder MA.Diagnosing and dealing with multicollinearity.West J Nurs Res,1990,12(2):175-187.

[28]Vatcheva KP,Lee M,McCormick JB,et al.Multicollinearity in Regression Analyses Conducted in Epidemiologic Studies.Epidemiology(Sunnyvale),2016,6(2).

(責(zé)任編輯：郭海強(qiáng))

Evaluation of Mediation Analysis Methods in Survival Data

Shi Qianwen,Wei Yongyue,Li Qingya,et al

(DepartmentofBiostatistics,SchoolofPublicHealth,NanjingMedicalUniversity(211166),Nanjing)

Objective To evaluate the five different mediation analysis methods for survival data including VanderWeele method,Baron and Kenny approach,Imai method,Sobel method and Inverse Weight method.Methods Statistical simulationswere performed under different parameter settings to evaluate the methods in type I error,power and analysis efficiency.Results Type I errors of all the five methods are close to the expected level(0.05)except the Inverse Weight method.Increased power was observed according to the increase of sample size,total effect,and proportion of mediation effect,while which was decreased according to the increased censoringrate.However,Inverse Weight method was inferior to the other four methods at the situation of small sample size(such as 100)and low mediation proportion(such as 30%).The analysis efficiency of Inverse Weight method,Baron-Kenny approach and Imai method were far slower than the other two.Conclusion Considering the type I error,power and test efficiency,we recommend VanderWeele method to analysis survival data.

Survival data;Mediation analysis;VanderWeele method

國家自然科學(xué)基金(81402764，81473070、81373102，81530088);江蘇省自然科學(xué)基金(BK20140907)

△通信作者：陳峰，E-mail：fengchen@njmu.edu.cn