疾病预后研究的中介分析方法评价*

南京医科大学公共卫生学院生物统计学系(211166)

施倩雯 魏永越 李清雅 段巍巍 赵 杨 陈 峰△



疾病预后研究的中介分析方法评价*

南京医科大学公共卫生学院生物统计学系(211166)

施倩雯 魏永越 李清雅 段巍巍 赵 杨 陈 峰△

目的 对用于癌症预后研究的五种中介分析方法(VanderWeele法、Baron-Kenny法、Imai法、Sobel法和InverseWeight法)进行评价，为实例分析的方法选择提供依据。方法 基于模拟试验，产生不同参数设置下的模拟数据，并评价五种方法的第一类错误、检验效能和分析时间。结果 除Inverse Weight法在相关系数较大时第一类错误有所膨胀外，其余四种方法的第一类错误在不同参数情况下均在0.05附近。五种方法的检验效能趋势一致，均随着样本量、中介比、总效应的增大而增大，随着删失比的增大而减小。在样本量较小(N=100)且中介比不大于30%的情况下，Inverse Weight法的检验效能低于另四种方法。Inverse Weight法、Baron-Kenny法和Imai法的分析效率远低于VanderWeele法和Sobel法。结论 综合考虑一类错误控制、检验效能及分析效率，推荐VanderWeele法进行预后研究的中介分析。

生存资料 中介分析 VanderWeele法

近年来，因环境污染加剧、生活方式改变，以及医疗水平提高等多方面原因，越来越多的癌症病人被确诊，癌症的治疗和预后也愈发引人关注。根据2015年的数据显示，中国男性发病率最高的癌症前三位依次为肺癌、胃癌和肝癌，女性则为乳腺癌、肺癌和结直肠癌[1,2]。以肺癌为例，查阅美国国立癌症机构统计的数据可知，其五年相对生存率只有17.7%[3]。影响其预后的因素有很多，包括吸烟、某些SNP位点等[4-7]，但其作用机制尚不清楚。而中介分析作为一种能够探索暴露对结局的内在作用机制的因果分析方法，受到了研究者的广泛关注[8-9]。

中介分析雏形最早在1986年由Baron & Kenny提出(以下简称Baron-Kenny法)。1987年，Sobel基于Baron-Kenny法提出了中介效应的方差估计方法(以下简称Sobel法)，计算简单、原理易于理解，但是较为保守[10]。近年来，VanderWeele(VanderWeele法)、Imai(Imai法)和Tchetgen(Inverse Weight法)等人进一步拓宽了中介分析的理论基础及应用范畴[11-13]。然而，至今尚无报道对以上五种方法进行比较评价。

本研究将以预后的生存资料为切入点，建立一个暴露因素和一个中介因素的中介模型，通过统计模拟试验，从第一类错误及检验效能角度对上述五种方法进行统计学性质评价，为预后的中介分析方法选择提供建议。

原理与方法

1.Baron-Kenny法

Baron-Kenny法的原理是先建立中介因素(M)、暴露因素(A)和结局(Y)之间的模型，然后通过模型中的系数计算间接效应(indirect effect，IE)，涉及的模型包括Y模型、M模型与总模型(公式1-3)[14]。以加速失效时间(acceleratefailure time，AFT)模型建立生存函数模型(即Y模型)，各模型如下：

Y模型：

(1)

M模型：

(2)

总模型：

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

(3)

据此，间接效应的点估计为：

IE=eθ1×β2

(4)

直接效应的点估计为：

IE=eβ1

(5)

间接效应表示暴露通过中介因素介导对结局产生的作用大小。中介效应占总效应的比例(即，中介比)表示为：

(6)

Baron-Kenny法中，间接效应的可信区间通过Bootstrap得到(抽样1000次)。

2.Sobel法

模型假设和IE的估计方法与Baron-Kenny法一致，log(IE)的方差可表示为：

(7)

其中Sθ1和Sβ2分别表示θ1和β2的标准误[10]。

3.VanderWeele法

VanderWeele法的原理是对总效应进行分解[15]，虽然同样用到Y模型和M模型，但相较于Baron-Kenny法，VanderWeele法考虑了暴露与中介因素之间的交互作用，故Y模型可表示为：

(8)

VanderWeele法的M模型表示为：

(9)

通过公式(8)和公式(9)的推导可知间接效应：IEAFT=e(θ1β2+θ1β3a)(a-a*)，直接效应：

4.Inverse Weight法

中介因素M对暴露A的回归模型为：

(10)

调整过权重的总模型：

(11)

5.Imai法

Imai法基于反事实原理，将平均总效应(averagetotaleffect)分成平均因果中介效应(ACMEs,averagecausalmediationeffects)和平均直接效应(ADE,averagedirecteffect)两部分[15,17-18]。

平均因果中介效应：

(12)

平均直接效应：

(13)

平均总效应(averagetotaleffect)：

(14)

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

模拟试验

1.模拟数据的产生

基于如图1的模型，按照以下步骤产生模拟数据：

(1) 指定总效应(HR)，根据中介效应比例(M%)把总效应(HR)分解成直接效应和间接效应两部分，分别可表示为：HR×(1-M%)和HR×M%，因DE=eβ1故可求得β1；

(2) 指定样本量(N)、A与M间的相关系数r，可产生二元标准正态分布随机数A和M；此时r即为公式(2)中的θ1值，因IE=eθ1×β2故可得到β2；

(3) 把变量A、M、系数β1和β2带入指数回归模型的风险函数得到参数λ，进而基于指数分布随机产生生存时间；

(4) 根据样本量(N)和删失比例(cen%)随机产生截尾，1代表死亡，0代表删失。

参考若干流行病学研究，总效应(HR)设为1.5、2.0、2.5、3.0[19-24]。由于生存资料的删失比越高，模型预测越不准确[25]，故保险起见，模拟试验时删失比取值不超过50%。理论上中介效应比例(M%)越大，中介分析的检验效能越高，而在实际分析中更希望能检出中介效应小的中介因素，不放过任何可能存在的因果关联，故中介比(M%)设置在50%及以下。考虑到预后研究的实际情况，模拟试验中样本量取100到1000。参数设置如表1。

表1 模拟试验参数设置情况

2.模拟试验流程

(1)根据参数设置产生模拟数据；

(2)用五种方法对产生的同一模拟数据分别进行分析，提取相关结果，包括：间接效应的点估计、区间估计(95%CI)、假设检验P值及检验结果(是否拒绝H0)；

(3)在不同参数情况下重复(1)～(2)步骤1000次，对1000次的结果进行汇总，得到五种方法的第一类错误和检验效能。

结 果

1.一类错误

当总效应不存在(HR=1)的时候，IE的期望为0。如图1为该情况下不同方法的第一类错误。横坐标为暴露A与中介M的相关系数r，纵坐标为第一类错误。当相关系数较低情况下，五种方法的第一类错误均接近于0，且随着相关系数的增大而逐渐增大。样本量较小情况下更为明显，样本量较大情况下，该趋势有所缓和。相关系数较大的情况下，IW法的第一类错误略有膨胀。删失比对第一类错误的影响不大。

图1 总效应不存在时暴露-中介的相关系数与第一类错误的关系

同理，当暴露与中介因素不存在相关性的时候(r=0)，IE的期望亦为0。图2为该情况下不同方法的第一类错误。横坐标为总效应，纵坐标为第一类错误。如图所示，在样本量较小(N=100)的情况下，五种方法的第一类错误随着总效应的增大逐渐增大；在样本量较大(N=500,1000)的情况下，该趋势不再明显，五种方法的第一类错误均保持在0.05附近。删失比对五种方法的第一类错误影响不明显。

2.检验效能

如图3所示(r=0.3，HR=1.5)，随着样本量的增加五种方法的检验效能均逐渐增大；IW法的检验效能明显低于其余四种方法。同样，如图4所示(r=0.3，N=100)，随着总效应的增加，检验效能亦逐渐增大；IW法亦明显低于其余四种方法。

当考察相关系数r不同水平对间接效应检验效能的影响时，发现随着相关系数的增大，检验效能先增大后减小(如图5)。当中介比为10%时，由于检验效能过低，无法反映出方法间的差异；而当中介比为30%及50%时，r在0.5以上时，检验效能逐渐下降，其中，中介比为30%、r在0.7以上时，检验效能较之峰值下降超过一半。

图2 暴露-中介关系为零时总效应与第一类错误的关系

图3 r=0.3，HR=1.5时检验效能与样本量的关系

其余参数情况下，趋势类似。

根据间接效应公式IE=eθ1×β2=er×β2可知，随着相关系数r的增大，间接效应理论上呈指数增加(在log尺度上呈线性增加)，其检验效能也应升高。但模拟结果与此推测相违背。结合模型(1)，提出假设：随着A和M相关系数的增大，模型(1)中A及M的系数估计的随机扰动增加，导致系数β2的变异变大，从而间接效应的变异变大，最终导致在较高相关系数的情况下检验效能减小。据此，设置模拟试验参数(β1=β2=0.5，samplesize=100，cen%=0.2，r=0.01(0.01)0.99，基于VanderWeele法模拟5000次，发现：随着相关系数的增大，β2的点估计波动变大，其方差呈现指数级膨胀；log(IE)的估计值呈线性上升，其方差亦呈指数级膨胀(图6)。该模拟试验结果验证了前文的假设。由此得出，当A与M相关性较高时(如，r>0.5)，五种中介分析方法皆存在检验效能较低的缺陷。

图4 N=100，r=0.3时检验效能与总效应的关系

图5 N=100，HR=2时检验效能与相关系数的关系

3.分析效率

理论上，由于Baron-Kenny法、IW法和Imai法皆采用Bootstrap进行区间估计，故分析效率较低。以样本量为1000的生存数据为例，用VanderWeele法分析仅需1.1073×10-2秒，用Sobel法需1.2705×10-2秒；如Bootstrap次数设置为1000次，则用Baron-Kenny法、IW法和Imai法分析分别所需时间为：12.2784秒、21.3041秒和9.2768秒，分析效率较低。故当样本量较大且维度较高的情况下，不推荐使用Baron-Kenny法、IW法和Imai法。

图6 相关系数与中介分析检验效能关系探讨结果

讨 论

除IW法，其余四种方法均能较好地控制第一类错误。五种方法的检验效能趋势一致，均随着样本量、中介比、总效应的增大而增大，随着删失比的增大而减小。总体上来看，IW法较其它方法检验效能较低，Baron-Kenny法和Imai法的分析效率远低于VanderWeele法和Sobel法，而Sobel法由于要求β2×θ1服从正态分布故在实际应用中存在较大的局限性[26]。

本研究亦提示，当X-M相关性较高(如r>0.5)情况下，各种方法皆存在检验效能下降的问题。残差法及主成分分析等方法可为后续研究提供思路[27-28]。在实例分析中，本研究以简单中介模型为基础进行模型试验，而在实际分析中，可能存在中介因素同样是生存数据的情况，往往还存在“多个暴露”、“多个中介因素”，或者“暴露与中介因素存在交互作用”等情况，在这些情况下该如何选择中介分析方法，仍需进一步探讨。

Baron-Kenny法和Sobel法虽然模型简单、原理便于理解，但这两种方法存在一定的局限性——不能处理暴露与中介因素存在交互作用以及暴露非线性的情况，而VanderWeele法则能较好地处理这些问题，故应用范围更广[13]。IW法对间接效应的点估计是通过总效应减去直接效应间接得到的，同样不适用于“暴露与中介因素存在交互作用”的情况[12]。故当“暴露与中介因素存在交互作用”的情况下，优先考虑VanderWeele法和Imai法，此时这两种方法的统计性质将是下一步的研究重点。

综上所述，VanderWeele法能较好地控制第一类错误，检验效能亦可，且分析效率具备一定的优势，故推荐作为生存数据的中介分析方法，尤其是维度比较高的情况下。当数据维度较低的情况下，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]曾倩,崔芳芳,宇传华，等.中国癌症发病、死亡现状与趋势分析.中国卫生统计.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).

(责任编辑：郭海强)

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

国家自然科学基金(81402764，81473070、81373102，81530088);江苏省自然科学基金(BK20140907)

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