某恶性肿瘤病人的生存情况分析-20150628

1、某恶性肿瘤病人的生存情况分析基于K-M、Cox和参数法三种模型课 程：生存分析 教 师： 学 校： SYSU 学 号：14212738、 14212775 2015年6月28日目 录一、背景概述2二、数据介绍2三、Kaplan-Meier模型3（一）K-M拟合结果3（二）不同治疗方法的生存率比较5四、Cox比例风险模型6（一）PHA检验61. 图形检验62. 统计量检验7（二）模型构建81. 单个变量的Cox模型82. 调整的Cox模型12五、参数模型15六、附录19一、背景概述生存分析是研究生存现象和响应时间数据及其统计规律的一门学科。该学科在生物学、医学、保险学、可靠性工程学、人口学、社会

2、学、经济学等方面都有重要应用。从定义上来讲，生存分析是指根据试验或调查得到的数据对生物或人的生存时间进行分析和推断，研究生存时间和结局与众多影响因素间关系及其程度大小的方法，因此也称为生存率分析或存活率分析。对于生存分析有一些基本的定义和性质。假设T是死亡时间，那么生存函数就定义为St=p(T>t)，定义域为t 0,)；而生存函数作为一个概率，St的范围为0,1之间。从上述定义中可以得到生存函数的如下几个性质：S0=1;limtSt=0;StaStbtatb;St=1-Ft=tfd.本文将利用患某恶性肿瘤病人的生存时间数据来分析不同治疗方法对生存率的影响以及其他例如年龄、性别等协变量因素

3、对生存率的影响，采用了K-M方法、Cox模型和参数模型三种方法。二、数据介绍该份数据共包含63条数据，8个变量，有关变量的符号表示及所代表的意义详见下表：表1：某恶性肿瘤的影响因素与赋值因素变量名赋值说明年龄X1（岁）性别X2女0，男1组织学类型X3低分化0，高分化1治疗方法X4传统疗法0，新型疗法1淋巴结转移X5否0，是1肿瘤浸润程度X6未突破浆膜层0，突破浆膜层1生存时间t（月）生存结局Y删失0，死亡1对数据的初步分析可以得到如下结论：Ø 数据中没有缺失值；Ø 时间范围为2,120；Ø 除去生存时间（t）和生存结局（Y）以外的6个协变量中，只有年龄（X1）为连

4、续性变量，其余5个协变量均为二分类变量；Ø 年龄（X1）的取值范围为19,74，均值为46.86，中位数为45；Ø 共有37个右删失值（Y=0），其中传统疗法15个，新型疗法22个。三、Kaplan-Meier模型（一）K-M拟合结果对数据进行K-M曲线拟合，得到如下所示图形。图1：K-M方法的整体生存曲线上图展示了该实验过程的整体生存曲线变化情况，横轴代表时间t，纵轴代表生存率S(t)。图中实现代表生存率的变化情况，虚线为95%的置信区间。由图可知，在40个月之前，生存率的降低比较明显，发生死亡事件的个体比较多；在40个月之后患者的情况变得较为稳定，除了有删失事件的发生，

5、几乎没有死亡事件发生。图2：K-M方法下不同治疗方法的生存曲线上图为所研究的对象不同治疗方法对患者的生存函数的影响曲线。图中，横坐标为时间（t），纵坐标为生存函数（S(t)），红色实线代表传统疗法的生存函数曲线，蓝色表示新型疗法的生存曲线；虚线为95%置信区间。从图中可以看到两种治疗方法的生存函数曲线存在着较大的差异，新型治疗方法取得了更好的治疗效果，下文中将进一步验证其是否存在显著差异。图3：两种治疗方法的K-M方法总结图由上图可以看到两种治疗方法随着时间的推进，死亡事件发生的情况统计。从上图可以得到以下结论：Ø 接受传统疗法（X4=0）的患者在2个月以后就开始发生了死亡事件，直到

6、40个月之前，死亡事件的发生率都十分高；Ø 接受新型疗法（X4=1）的患者在12个月后才发生死亡时间，之后分别在40、66、120个月有一名患者死亡。Ø X4=0时，总体的生存率为S1t=0.405；95%的置信区间为0.274,0.599；Ø X4=1时，总体的生存率为S1t=0.793；95%的置信区间为0.617,1；Ø 结合前面的生存曲线图看，无论是哪种疗法在40个月之后的生存曲线都变得较为稳定。（二）不同治疗方法的生存率比较假设H0:S1=S2，即两种治疗方法下的患者生存率是一样的，则通过R中的survdiff( )函数可以对该假设做检验。该函

7、数中包含一个参数rho：表示分别给不同时间下的死亡个体赋予不同的权重，取值分别有0、1、0.5，分别对rho取三个不同参数的情况做假设检验，得到的结果如下所示：图4：rho=0、1、0.5时的两种治疗方法假设检验上述3图中，rho=0时，p=0.000264；rho=1时，p=0.000179；rho=0.5时，p=0.000206；三个检验结果的p<0.05，因此有充分的理由拒绝原假设，即说明两种不同的治疗方法是有显著差异的。四、Cox比例风险模型（一）PHA检验在构架Cox比例风险模型前，首先对变量是否满足PHA假定做检验，方法包括图形检验和统计量检验两种方式。1. 图形检验图形检验

8、的思想是假设协变量是符合PHA假定的，在这个假定下画出各个协变量的生存函数曲线并与K-M曲线进行比较，若两者比较接近，则说明PHA的假定是正确的；否则就是不正确的。在此需要说明的是“年龄”协变量X1为数值型变量，因此首先将其进行分段成（age>=45 和age<45）两段，最终得到的6个协变量的图如下所示： 图5：各个协变量在PHA假定下的Cox曲线与K-M曲线比较从上述六图中可以发现，每个协变量在PHA假定下的Cox曲线与K-M曲线都较为接近，差异比较大的是X2、X4和X6，由此可见，仅凭图形检验的方式并不能得出十分明朗的结论，因此进一步对每个协变量进行统计量检验。2. 统计量检

9、验对于每个协变量，假设H0:Xi满足PHA假设，i=1,6；即每个协变量可以表示Cox模型。对模型的检验结果如下：图6：6个协变量的卡方检验从上图中可以看到，X1-X5的p值都大于0.05，即接受原假设；而X6的的p=0.0204<0.05，拒绝原假设。由此说明X1-X5都可以表示成Cox模型，而X6则不可以，由此在模型构建时考虑X6与时间t的交互项影响。 （二）模型构建1. 单个变量的Cox模型数据中共包含6个协变量，首先分别对每个协变量构建cox比例风险模型的初步探索。对每个协变量的原假设为H0:i=0,i=1,6，即假设每个协变量的系数值都分别等于0；模型假设为：ht,xi=h0t

10、,expxi,i=1,6.对模型的显著结果有3种检验方法：Likelihood ratio test、Wald test和Score (logrank) test。下面分别展示了每个协变量的cox模型结果。（1）年龄（X1）图7：“年龄”的Cox系数显著性检验从3种检验结果来看，p=0.52左右，远大于0.05的置信水平，故不能拒绝系数为0 的原假设，也就是说“年龄（X1）”在Cox模型中不合适。（2）性别（X2）图8：“性别”的Cox系数显著性检验从3种检验结果看，p<0.05，因此模型是显著的；另一方面，对X2的系数检验结果p=0.00171也是显著的，由此可见，该Cox模型系数拟合

11、结果，模型为：ht,x=h0t,exp0.2289x,x=0,1其中x=0为女性，x=1为男性，由此说明当其他因素相同时，在任一时刻男性的即刻死亡风险是女性的4.369倍。（3）组织学类型（X3）图9：“组织学类型”的Cox系数显著性检验从上图的三个检验结果可以看到，p值约为0.08，略大于0.05；X3的系数显著性检验p=0.0875，同样略大于0.05，因此暂且判断其的显著性不明显，接受原假设。（4）治疗方法（X4）图10：“治疗方法”的Cox系数显著性检验从3种检验结果看，p<0.05，因此模型是显著的；另一方面，对X2的系数检验结果p=0.0012也是显著的，由此可见，该Cox模

12、型拟合效果显著，模型为：ht,x=h0t,exp-1.769x,x=0,1其中x=0为传统疗法，x=1为新型疗法，由此说明当其他因素相同时，接受传统疗法的患者的即刻死亡风险是接受新型疗法患者的5.865倍。（5）淋巴结转移（X5）图11：“淋巴结是否转移”的Cox系数显著性检验从3种检验结果看，p<0.05，因此模型是显著的；另一方面，对X2的系数检验结果p=0.0369也是显著的，由此可见，该Cox模型拟合效果显著，模型为：ht,x=h0t,exp0.9240x,x=0,1其中x=0为淋巴结未转移，x=1为淋巴结已转移，由此说明当其他因素相同时，淋巴结未转移患者的即刻死亡风险是已转移患

13、者的2.519倍。（6）肿瘤浸润程度（X6）图12：“肿瘤浸润程度”的Cox系数显著性检验从上图的三个检验结果可以看到，p值约为0.56，远大于0.05；X6的系数显著性检验p=0.563，同样显著大于0.05，因此判断其的显著性不明显，不拒绝原假设。2. 调整的Cox模型在进过对每个变量的初步检验后，用“逐步回归”的方法对所有变量进行考虑，并包含交互项的作用。得到的结果如下图所示：图13：利用“逐步回归”方法进行Cox模型拟合过程图上图显示了“逐步回归”在筛选变量进入模型时的过程，每一次都将能最大幅度减小AIC值的变量加入到模型中去，最后一个模型的AIC值最小为182.31，模型中包含的变量

14、包括X4、X5、X2和X4:X5交互项。在之前的判断中认为X6与时间t可能存在交互项作用，因此在上述模型的基础上加入X6以及X6和时间t的交互作用，得到的结果如下所示：图14：Cox调整模型结果1上图中三种检验方法全都通过，然而X4和X2的系数检验不显著，因此进一步考虑删去X2和X4的情况。首先删去X2，得到的结果如下图：图15：Cox调整模型结果2删去X2后X4的系数检验仍然无法通过，因此进一步删去X4进行建模，得到如下结果：图16：Cox调整模型结果3此时模型的三种检验结果显著，每个变量的系数也十分显著，因此选取该方法为Cox模型的最终结果，在对上述3个模型的AIC值进行比较时，也得到了同

15、样的结果，模型3的AIC值最小。表2：三个调整模型的AIC值模型AIC值Em1182.2686Em2181.9716Em3180.1866因此最终的Cox模型为：hx=h0(x,)exp1.35x4+3.80x6-2.74x4x5-1.01x6logt五、参数模型常用参数生存分析模型：Exponential：St,x1,x2=exp-t exp-0-1x1-2x2Weibull：St,x1,x2=exp-t exp-0-1x1-2x2Log-logistic：St,x1,x2=exp1+t exp-0-1x1-2x2-1对参数模型的选取同样可以通AIC准则来实现，下面分别选取三种模型，筛选变量

16、使得模型AIC值达到最小。1. Exponential模型：先建立只与时间有关的初始模型，然后通过stepAIC函数来选取变量，得到AIC值最小的Exponential模型。exp=survreg(S1,dist="exponential",data=data)Scope=list(upper=(X1+X2+X3+X4+X5+X6)2,lower=1)stepAIC(exp,Scope,direction="both")最终得到的模型为exp=survreg(formula = S X4 + X5 + X2 + X6 + X4:X5, data = da

17、ta, dist = "exponential")该模型AIC值为289.8446。 结果显示如下图：图17：指数模型拟合结果从结果可以看出X2、X4、X5、X6、X4与X5的交互项都表现显著，其中模型p值为3.1×10-7。最终模型可写成下列表达式：St,x1,x2=exp-t exp-5.613-1.116X2+0.574X4+1.861X5+0.659X6-3.194X4X5 2. Weibull模型：同样先建立只与时间有关的初始模型，然后通过stepAIC函数来选取变量，得到AIC值最小的Weibull模型。wei=survreg(S1,dist=&quo

18、t;w",data=data)Scope=list(upper=(X1+X2+X3+X4+X5+X6)2,lower=1)stepAIC(wei,Scope,direction="both")wei=survreg(S X4 + X5 + X2 + X4:X5,dist="w",data=data)最终得到的模型为：survreg(formula = S X4 + X5 + X2 + X4:X5, data = data, dist = "w")该模型AIC值为288.8021。 结果显示如下图：图18：威布尔模型拟合结果从

19、结果可以看出X2、X4、X5、X4与X5的交互项都表现显著，相比Exponential模型少了X6变量，其中模型p值为2.4×10-7。最终模型可写成下列表达式：St,x1,x2=exp-t1.31 exp-5.640-1.420X2+0.697X4+2.1121X5-3.988X4X51.31 3. Log-logistic模型：同样先建立只与时间有关的初始模型，然后通过stepAIC函数来选取变量，得到AIC值最小的Log-logistic模型。loglog=survreg(S1,dist="loglogistic",data=data)Scope=list(

20、upper=(X1+X2+X3+X4+X5+X6)2,lower=1)stepAIC(loglog,Scope,direction="both")loglog=survreg(S X4 + X5 + X2 + X4:X5,dist="loglogistic",data=data)最终得到的模型为survreg(formula = S X4 + X5 + X2 + X4:X5, data = data, dist = " loglogistic ")该模型AIC值为283.4984。 结果显示如下图：图19：log-logistic模型

21、拟合结果从结果可以看出X2、X4、X5、X4与X5的交互项都表现显著，与Weibull模型一样，其中模型p值为7.4×10-7。最终模型可写成下列表达式：St,x1,x2=exp1+t0.906 exp-4.8995-1.2766X2+0.5003X4+2.0706X5-3.5709X4X50.906-1 从上面三个模型可以看到比较显著的变量包括X2、X4、X5、X6、X4与X5的交互项，这与前面非参模型得到的结论一致，即性别、治疗方法、是否有淋巴结转移、肿瘤浸润程度对患者的生存率有显著影响。三者的AIC值以Log-logistic模型最小，为283.4984。由于无法利用AIC值来

22、对参数模型和cox模型进行比较，故此处只限参数模型的比较。六、附录library(survival)data=read.csv("g:/数据.csv",head=T)summary(data)table(data$X4,data$Y) #删失数据在不同治疗方法中的个数#K-M曲线#S=Surv(data$t,data$Y)fit0=survfit(S1,data=data)plot(fit0)fit=survfit(Surv(t,Y)X4,data=data) #构建K-M模型plot(fit,xlab="t",ylab=expression(hat(S

23、)(t), = F,col=c(2,4)#plot(fit,xlab="t",ylab=expression(hat(S)(t), = F,lty=1:2,col=c(2,4)legend("bottomleft",c("传统疗法","新型疗法"),lty=1,col=c(2,4)summary(fit) #得到两个类别的K-M模型结果，有死亡时间，死亡个体数，生存概率，标准差，CIsurvdiff(Surv(t,Y)X4,data=data,rho=1) #做H0：S1(t)=S2

24、(t)检验,rho=1在时间值小的部分加更多权重survdiff(Surv(t,Y)X4,data=data,rho=0) #做H0：S1(t)=S2(t)检验,rho=0在时间值大的部分加更多权重survdiff(Surv(t,Y)X4,data=data,rho=0.5) #做H0：S1(t)=S2(t)检验,rho=0.5在时间值小的部分加更多权重#cox模型#图形检验fit=list()fit_X1=survfit(Surv(t,Y)X1,data=data) #构建K-M模型fit1=survfit(Surv(t,Y)(X1<45),data=data) fit2=survfi

25、t(Surv(t,Y)X2,data=data) #构建K-M模型fit3=survfit(Surv(t,Y)X3,data=data) #构建K-M模型fit4=survfit(Surv(t,Y)X4,data=data) #构建K-M模型fit5=survfit(Surv(t,Y)X5,data=data) #构建K-M模型fit6=survfit(Surv(t,Y)X6,data=data) #构建K-M模型phm=list()phm1=coxph(Surv(t,Y)X1,data=data) #年龄phm2=coxph(Surv(t,Y)X2,data=data) #性别phm3=co

26、xph(Surv(t,Y)X3,data=data) #组织性类型phm4=coxph(Surv(t,Y)X4,data=data) #治疗方法phm5=coxph(Surv(t,Y)X5,data=data) #淋巴结转移phm6=coxph(Surv(t,Y)X6,data=data) #肿瘤浸润程度pdf("e:/v4.pdf") #画年龄X1的图形plot(fit1,xlab="t",ylab=expression(hat(S)(t), = F,col=c(2,4)d.phm=coxph.detail(phm1)times=c(0

27、,d.phm$t)h0=c(0,d.phm$hazard)S0=exp(-cumsum(h0)beta=c(phm1$coef)x1=mean(data$X1data$X1>=45)-mean(data$X1)Sx1=S0exp(t(beta)%*%x1)x2=mean(data$X1data$X1<45)-mean(data$X1)Sx2=S0exp(t(beta)%*%x2)lines(times,Sx1,col=2,type='s',lty=2)lines(times,Sx2,col=4,type='s',lty=2)legend("

28、bottomleft",c("age>=45","age<45"),lty=1,col=c(2,4),cex=1) #画X2-X6的图形d.phm=list()for(i in 2:6) plot(fiti,xlab="t",ylab=expression(hat(S)(t), = F,col=c(2,4) d.phmi=coxph.detail(phmi) times=c(0,d.phmi$t) h0=c(0,d.phmi$hazard) S0=exp(-cumsum(h0) beta=c(ph

29、mi$coef) x1=c(0)-mean(data,i) Sx1=S0exp(t(beta)%*%x1) x2=c(1)-mean(data,i) Sx2=S0exp(t(beta)%*%x2) lines(times,Sx1,col=2,type='s',lty=2) lines(times,Sx2,col=4,type='s',lty=2) legend("bottomleft",c(paste(names(data)i,"=0"),paste(names(data)i,"=1"),lty=1,c

30、ol=c(2,4),cex=1)dev.off()#图形检验完毕#pha检验phm.data=coxph(SX1+X2+X3+X4+X5+X6,data=data)#初步coxpha.data=cox.zph(phm.data,transform='rank',global=F)#pha检验pha.data#模型构建phm=list()phm1=coxph(Surv(t,Y)X1,data=data) #年龄phm2=coxph(Surv(t,Y)X2,data=data) #性别phm3=coxph(Surv(t,Y)X3,data=data) #组织性类型phm4=coxp

31、h(Surv(t,Y)X4,data=data) #治疗方法phm5=coxph(Surv(t,Y)X5,data=data) #淋巴结转移phm6=coxph(Surv(t,Y)X6,data=data) #肿瘤浸润程度summary(phm1)summary(phm2)summary(phm3)summary(phm4)summary(phm5)summary(phm6)#筛选变量library(MASS)attach(data)Scope=list(upper=(X1+X2+X3+X4+X5+X6)2,lower=1)phm_0=coxph(Surv(t,Y)1)phm_f=stepAIC(phm_0,Scope,direction="both")detach(data)fitstep1=coxph(SX4 + X5 + X2 + X4:X5,data=data)summary(fitstep1)#extending modelem1=coxph(SX4 * X5 + X2 +X6+X6:log(t),data=data)#summ