多结局生存分析模型与Cox模型的随机模拟比较

复旦大学公共卫生学院社会医学与卫生统计教研室　高峻　董　伟　赵耐青
上海市计划生育科学研究所　高尔生

【提　要】　目的　通过随机模拟评价不同多结局生存分析模型的特点。方法　利用随机模拟数据,比较多结局生存分析模型、将多个结局视为一种结局的单结局Cox模型及将各个结局分开拟合多个单结局Cox模型的回归系数的估计精度。结果　多结局生存分析模型回归系数的估计最准确,95%可信区间包含回归系数的百分比最高,95%可信区间不包含0的百分比也最高。结论　采用将各个结局分开单独拟合Cox模型和将多种结局视为单一结局拟合Cox模型分析多结局生存数据会导致回归系数估计不准及检验效能的降低。

【关键词】　多结局生存分析，夹心方差估计，风险区间

　　多结局现象(Multiple Failure Time Data)在各个领域都非常常见。所谓的多结局现象就是指研究对象可能出现的失效事件有多种,而不单单是经典生存分析中只有一种失效事件(即研究者所关心的结局不止一个)。研究中出现的多种结局类型可能是相同的,也可能是不同的。在多结局数据中,同一对象出现的多个不同结局所对应的生存时间之间存在着相关性〔1〕。
　　对于多结局的模型,偏似然估计对参数估计是渐近无偏的,但是由于多结局之间的相关性,导致参数估计的标准误偏小,因此采用”夹心方差估计”(sandwich variance estimator)〔2-4〕进行校正。经典的Cox模型回归系数的标准误是通过相应模型的方差协方差矩阵对角线各个数值的平方根而计算出。通过这种方法得到的标准误称为纯标准误(naive standard error)。但当同一个体的各生存时间出现相关性的时候该标准误不稳定〔5-11〕,此时计算出的标准误高估了样本中相关对象的个数,或者也可以认为“忽略了组内相关性,从而高估了组间的信息”〔12〕。这时上述方法估计出的标准误偏小。“调整夹心方差估计”(modified sandwich)就是为了校正这种组内的相关性,这种方法是Lin和Wei〔13〕基于Huber〔14〕和White〔15〕提出的夹心方差估计的基础上进一步提出的。目前所有的模型都采用“调整经验夹心方差估计”〔16-18〕(modified empirical sandwich variance estimator)来估计标准误,该方法在“调整夹心方差估计”的基础上进行了进一步的改进,使标准误的大小更接近实际情况。
　　由于Cox比例风险模型为半参数模型,适用范围比较广,所以大多数多结局生存分析模型都是多结局的Cox比例风险模型,并且众多的多结局比例风险模型都是基于以下两个最基本的模型:
　　　　基本模型1:hk(t;Xk)=h0(t)exp(β1x1k+β2x2k…+βpxpk) (1)
　　　　基本模型2: hk(t; Xk)=h0k(t)exp(β1x1k+β2x2k…+βpxpk) (2)
　　对于基本模型1〔19〕,假定不同的结局对应基础风险函数相同,而对于模型2〔20〕,假定不同结局对应基础风险函数不同。现在发展得比较成熟的多结局生存分析模型都是以这两个基本模型为基础,通过定义不同时间变量而产生各种不同的多结局模型,针对各种不同的实际情况,解决相应的问题。
时间变量有三种定义的形式(由结局开始观察的时间和该结局结束观察的时间所构成的一个时间区间称为风险区间〔21〕)。
图1　各个不同的生存时间变量定义方式
第一种时间变量的定义方式被称为“total time”。
　　各种结局之间为无序的情况,一般采用这种时间变量定义方式,称为无序生存分析模型。第二种定义方式被称为“counting process time”,第三种定义方式被称为“gap time”,当各种结局之间为有序的情况时,可采用“counting process time”和“gap time”定义生存时间变量,故称为有序生存分析模型。
　　由两个基本模型基础上根据不同的应用条件许多学者提出了各种不同的多结局生存分析模型,现有比较成熟的模型有PWP模型〔22〕(根据模型采用的风险区间的不同可以细分为PWP-CP模型和PWP-GT模型)、AG模型〔11〕、GTUR模型〔23〕、WLW模型〔24〕和LWA模型〔25〕。
　　上述六个多结局生存分析模型可以按照基础风险函数和风险区间两个因素进行如下划分(表1)。
表1　多结局生存分析模型的划分
【基础风险函数 生存时间定义 有序 Counting process time Gap time 无序 Total time 相同AG模型GT-UR模型LWA模型
不同PWP-CP模型PWP -GT模型WLW模型】
　　 多结局的生存分析模型在回归系数的估计和标准误的估计上都考虑了相关性的影响,因此估计值更接近研究因素实际的作用,但目前多结局生存分析在实际应用中没有得到很大的推广。为了进一步验证多结局生存分析的优越性,本文利用随机模拟产生多结局生存数据,采用三种不同的分析方法对模拟数据进行分析,通过比较结果之间的差别来说明多结局生存分析的重要性及不同模型应用的适用性。
方　　法
　　多结局模型的随机模拟中,产生随机数据的回归模型中性别和年龄的设定值分别为0.1和0.2,为了简化起见,模拟生存时间服从指数分布。
基础风险函数相同时,对应三次复发的风险函数为:
　　对应第一个结局:h(t,X)=h(t)exp(-14+0·1sex+0·2age)
　　对应第二个结局:h(t,X)=h(t)exp(-13+0·1sex+0·2age)
　　对应第三个结局:h(t,X)=h(t)exp(-12+0·1sex+0·2age)
基础风险函数不同时,对应三次复发的风险函数为:
　　对应第一个结局:h(t,X)=h01(t)exp(-14+0·1sex+0·2age)
　　对应第二个结局:h(t,X)=h01(t)exp(-13+0·1sex+0·2age)
　　对应第三个结局:h(t,X)=h01(t)exp(-12+0·1sex+0·2age)
　　本次研究共模拟了四种情况,即基础风险函数相同的有序多结局情况、基础风险函数相同的无序多结局情况、基础风险函数不同的有序多结局情况和基础风险函数不同的无序多结局情况。在四种情况中对应各种结局设定10%的失访率。样本含量500,重复抽样次数500次。
在比例风险模型中,生存时间服从指数分布是指固定协变量的情况下,因此其概率模型为:
其生存函数,S(t)=exp(-texp(βX))
其风险函数为exp(βX)
β=(β0,β1,…,βp),X(1,x1,x2,…,xp)τ
　　数理统计学上可以证明:设Z在(0,1)上均匀分布,则固定X,随机变量t=-ln(Z)exp(-βX)服从均数为exp(βX)的指数分布,并称上述模型为指数回归模型,属于比例风险模型。因此只要设定β和X,就可以用计算机产生在(0,1)上均匀分布,就可以得到服从指数分布随机模拟数据。
　　用统计软件产生服从标准正态分布的随机模拟数据X和Y,且相关系数为0·2(或0·3等),设对应标准正态分布的分布函数为F(X),则Z1=F(X),Z2=F( Y),则根据数理统计可知,Z1和Z2在(0,1)上呈均匀分布,其Spearman相关系数意义下的总体相关系数为0·2(或0·3等),然后用上述方法可以产生相关的生存时间随机模拟数据,并且服从指数分布模型。
2·分析方法
基础风险函数相同时,对应有序的多结局数据采用AG模型和GTUR模型拟合,对应无序的多结局数据采用LWA模型拟合。基础风险函数不同时,对应有序的多结局数据采用PWP-CP模型和PWP-GT模型拟合,对应无序的多结局数据采用WLW模型拟合。 此外,采用下列两种方法分析:
　　Cox模型1:将各种结局分开,认为三种结局之间没有任何联系,分别对应不同结局独立拟合经典Cox模型。
　　Cox模型2:将各种不同的结局视为一种,拟合单结局Cox模型。即每抽样一次,对应一个对象抽样三次。
结　　果
多结局生存分析模型(AG模型、GTUR模型、LWA模型)估计出的回归系数偏差很小,计算出回归系数95%的可信区间包含设定值的百分比除AG模型sex在96·20%外,其余都在95%左右,不包含0的百分比在40%以上。
将所有结局视为一种的单结局Cox模型估计出的回归系数偏差也很小,计算出回归系数95%的可信区间包含设定值的百分比都在95%左右,sex不包含0的百分比也在40%左右。
在将三种结局分割为对应三种结局的单结局Cox模型得到的结果中,回归系数的估计精度较多结局生存分析模型和将所有结局视为一种的单结局Cox模型相差不大,也较为精确。计算出回归系数95%可信区间包含回归系数估计值的百分比在95%左右,但sex不包含0的百分比在Cox模型1较低,在20%以下。
表2　有序多结局生存分析模型与单结局Cox模型的拟合效果之比(基础风险函数相同,性别对应回归系数为0·1,年龄为0·2,样本含量:500)
【模型变量|β-b|×103β∈95%CI(%) 0 95%CI(%) AG模型sex 4.0948 96.20 44.00age 0.5692 95.40 100.00 GTUR模型sex 3.9969 95.00 43.00age 0.0683 95.40 100.00 Cox模型1sex 3.0101 95.60 17.87age 0.4878 95.13 100.00 Cox模型2sex 3.9802 95.40 42.60age 0.0962 95.40 100.00】
　　 *性别对应的回归系数为0·1,年龄对应的回归系数为0·2
　　 *|β-b|表示回归系数估计值与设定值之差的绝对值的平均数
表3　无序多结局生存分析模型与单结局Cox模型的拟合效果之比(基础风险函数相同,性别对应回归系数为0.1,年龄为0.2,样本含量:500)
【模型变量|β-b|×103β∈95%CI(%) 0 95%CI(%) LWA模型sex 1.3478 95.20 45.40age 0.0049 96.00 100.00 Cox模型1sex 0.9513 95.67 18.33age 0.5271 95.67 100.00 Cox模型2sex 1.3198 95.40 45.20age 0.0472 96.20 100.00】
　　 *性别对应的回归系数为0·1,年龄对应的回归系数为0·2
　　 *|β-b|表示回归系数估计值与设定值之差的绝对值的平均数
表4　有序多结局生存分析模型与单结局Cox模型的拟合效果之比(基础风险函数不同,性别对应回归系数为0.1,年龄为0.2,样本含量:500)
【模型变量|β-b|×103β∈95%CI(%) 0 95%CI(%) PWP-CP模型sex 5.4837 95.00 44.20age 0.7935 95.00 100.00 PWP-GT模型sex 3.8173 95.40 43.20age 0.1074 95.00 100.00 Cox模型1sex 3.1511 95.73 17.93age 0.4315 95.20 100.00 Cox模型2sex 20.3448 95.40 23.80age 52.7882 0.00 100.00】
　　 *性别对应的回归系数为0·1,年龄对应的回归系数为0·2
　　 *|β-b|表示回归系数估计值与设定值之差的绝对值的平均数
　　 多结局生存分析模型(PWP-CT模型、PWP-GT模型)回归系数估计较精确,回归系数95%可信区间包含设定值的百分比在95%左右,回归系数95%可信区间不包含0的百分比在40%以上。
将所有结局视为一种的单结局Cox模型得到的结果中,回归系数的估计值偏差最大, 95%可信区间包含设定值的百分比Cox模型2 age最低为0, 95%可信区间sex不包含0的百分比在22%左右。
在将三种结局分割为对应三种结局的单结局Cox模型得到的结果中,回归系数的估计也较为精确,95%可信区间包含设定值的百分比在95%左右,95%可信区间不包含0的百分比最低,sex在20%左右。
表5　无序多结局生存分析模型与单结局Cox模型的拟合效果之比(基础风险函数不同,性别对应回归系数为0.1,年龄为0.2,样本含量:500)
【模型变量|β-b|×103β∈95%CI(%) 0 95%CI(%) WLW模型sex 1.117 95.00 45.00 age 0.0136 95.40 100.00 Cox模型1 sex 1.0164 95.60 18.40 age 0.3807 95.67 100.00 Cox模型2 sex 24.0655 96.00 22.00 age 54.8731 0.00 100.00】
　　 *性别对应的回归系数为0·1,年龄对应的回归系数为0·2
　　 *|β-b|表示回归系数估计值与设定值之差的绝对值的平均数
讨　　论
　　1·当基础风险函数相同且各研究因素对应不同结局的作用相同(影响因素对应不同结局回归系数设定值也相同)时,多结局生存分析模型较将所有结局视为一种结局拟合单结局Cox模型的优势在于考虑到各个不同结局之间的相关性。在本次模拟中各个时间变量之间的相关性较小,各时间之间的秩相关系数在0·2左右,因此两者的结果较为相近。
　　2·从模拟结果可以看出,将各个结局单独拟合Cox模型犯二类错误的可能性会增大,即检验效能降低,可能会遗漏那些实际上对失效事件的发生有影响的因素。而将所有结局视为一种拟合Cox模型犯二类错误的可能性比各种结局单独拟合Cox模型要低一些,但这种方法导致的是回归系数估计的较大偏差。
　　3·按照基础风险函数和风险区间两个因素可以将六种多结局生存分析模型分类。其中PWP-CP模型、PWP-GT模型、AG模型和GTUR模型为有序的多结局模型,WLW模型和LWA模型为无序的多结局模型。有序和无序多结局模型是针对同一对象发生的多种不同结局而言。有序就是指同一对象出现的不同结局之间存在一定的顺序,例如膀胱癌的第二次复发一定是在患者第一次复发以后出现。而在无序的情况下从研究开始时所有的对象在观察时间内发生中结局中的任意一种,例如糖尿病患者可能出现的不同的并发症(视网膜病变、神经病变等等)。故有序和无序的生存分析模型是对应两种完全不同的情况,不可以相互混用。基础风险函数的相同与否则必须根据实际的研究背景决定,不同的结局基础风险函数肯定不同,但相同结局的基础风险函数也可能不同,例如以孕妇出现流产为结局的研究中对象多次流产后导致的习惯性流产,此时发生相同结局(流产)的风险逐渐增大。
　　4·综上所述,对于多结局的生存数据应该根据实际情况采用合适的多结局生存分析模型来描述整个生存过程和寻找有关影响因素。

Comparison Between Multivariate Generalization of the Cox Proportional Hazards Models and Cox Proportional Hazards Model
by Simulation　Gao Jun,Dong Wei, Gao Ersheng,et al.School of Public Health,Fudan University(200032),Shanghai
【Abstract】　Objective　Evaluate the characteristics of multivariate generalization of the Cox proportional hazards models by simulation.
Methods　Compare estimate biases of three strategies of analysis (fitting multivariate generalization of the Cox proportional hazards models,fitting Cox proportional hazards model by assuming all the outcomes are same,fitting Cox proportional hazards model by dividing different outcomes) by simulated data.Results　Estimators of multivariate generalization of the Cox proportional hazards models are the most closed to it’s real values,the proportion of including estimator and the proportion of not-including zero in 95%CI of according estimator is the highest in three analysis strategies.
Conclusion　Fitting Cox proportional hazards model by assuming all the outcomes are same and fitting Cox proportional hazards model by dividing different outcomes will conduce great biases in estimating coefficients and lower the efficiency.
【Key words】　Multivariate survival analysis;Sandwich variance estimator;Risk interval
参　考　文　献
1·Hougaard Philip.Analysis of Multivariate Survival Data.New York,Springer-Verlag,2000.
2·Crowley JJ.A comparison of survival life table estimates.Master’s thesis,University of Washington.Seattle,1970.
3·Lawless JF.Statistical Models and Methods for Life Time Data.Wiley,New York,1982.
4·Cox D R.and Oakes D.Analysis of Survival Data.London,Chapman and Hall,1984.
5·Mann N R,Schafer R.E and Singpurwalla N.D.Method for Statistical Analysis of Reliability and Life data,Wiley,New York,1974
6·Nelson W B.Hazard plotting for incomplete failure data. Journal of Quality Technology,1969,1(3):27-52.
7·Kalbfleisch JD and Prentice RL.The statistical analysis of failure time data.Wiley,New York,1980.
8·Andersen PK, Borgan O, Gill RD and Keiding N. Statistical Models based on Counting Processes.London,Spring-Verlag,1993.
9·Cox D R.Regression models and Life Tables.Journal of Royal Statistical Society(B),1972,34:187-220.
10·Cox DR.Partial Likelihood.Biometrika,1975,62:269-276.
11·Anderson PK and Gill RD.Cox regression model for counting processes:A large sample Study.The Annals of Statistics, 1982, 10: 1100-1120.
12·Aalen OO.Nonparametric inference for a family of counting process.The Annals of Statistics,1978,6:701-726.
13·Wei LJ,Lin DY and Weissfeld L.Regression analysis of multivariate incomplete failure time data by modeling marginal distribution.Journal of America Statistical Association,1989,84:1065-1073.
14·Lin DY and Ying Z.Semi-parametric analysis of the additive risk model.Biometric,1994,81:61-71.
15·Nelson WB.Graphical analysis of system repair data.Journal of Quality Technology 1988,1:27-52.
16·Kelly Patrick Jand Lynette LY.Lim.Survival Analysis for Recurrent Event Data:An Application to Childhood Infectious Diseases.Statistics in Medicine,1999,19:13-33.
17·Therneau Terry M.and Patricia M Grambsch.Modeling Survival Data:Extending the Cox Model.New York,Springer-Verlag,2000.
18·Finkelstein Dianne M,David A Schoenfeld and Eva Stamenovic.Analysis of Multiple Failure Time Data from an Aids Clinical Trial.Statistics in Medicine,1997,16:951-961.
19·Gross A J andClark V A.Survival Distribution:Reliability Applications in the Biomedical Sciences.Wiley,New York,1975.
20·Mann N R,Schafer R E and Singpurwalla N D.Method for Statistical Analysis of Reliability and Life data.New York,Wiley,1974.
21·Kalbfleisch JD andPretience R L.Marginal Likelihoods based on Cox regression and Life Model.Biometrika,1973,60:267-278.
22·Nelson WB.Hazard plotting for incomplete failure data.Journal of Quality Technology,1969,1:7-52.
23·Johenson S.The product limit estimator as maximum likelihood estimator.Scandinavian Journal of Statistics,1983,5:165-174.
24·Pepe MS and Fleming TR. Weighted Kaplan-meier statistics: Large sample and optimality considerations.Journals of the Royal Statistical Society,Series B,1991,53:341-352.
25·Lawless JF.Statistical Models and Methods for Life Time Data.New York,Weily,1982.