基于原发性胆汁性肝硬化数据的生存分析

1、基于原发性胆汁性肝硬化数据的生存分析一、 数据背景1、数据名称：原发性胆汁性肝硬化的数据2、数据来源：/datasets/pbc3、数据简介：该数据以pbc的名字包含在程序包survival中，数据集有20个变量，共有418个患者的观测,每个患者只有一次观测。其数据是来自于梅奥诊所所审理的从1974年到1984年这十年间的原发性胆汁性肝硬化的病例。总共有424个PBC患者，按照资格标准随机采用安慰剂对照药物D-青霉胺进行。312个个案在数据集中参加随机试验，而附加的112个个案没有参加临床试验，但有基本的测量记录。6个个案失去了确诊后没有跟进，所以数

2、据上附加的是106个个案以及312个随机参与者，共418组观测值。4、变量介绍：变量名描述性质id个人代码分类time登记到换肝或死亡天数数量status状态（活着，肝移植，死亡)分类trt用药（D-青霉胺，安慰剂）分类age年龄数量sex性别分类ascites腹水（No,Yes)分类hepatomegaly肝肿大（No,Yes)分类spiders蜘蛛纹（No,Yes)分类edema水肿（无水肿，未用利尿剂水肿，用利尿剂仍水肿）分类bili血清胆红素（mg/dl)数量chol血清胆固醇（mg/dl)数量albumin白蛋白（mg/dl)数量copper尿铜（ug/day)数量alk.phos碱

3、性磷酸酶（U/liter)数量ast血清谷氨酸谷草转氨酶（U/liter)数量trig甘油三酯（mg/dl)数量platelet血小板（每升千个）数量protime凝血酶原（秒）数量stage疾病组织学分期（1,2,3,4期）分类5、分析思路：如果对后面106个观测值的缺失值整列进行弥补然后对418个观测值进行分析意义不大，因此只将前312个观测值单独提取出来进行非参数、参数、Cox比例风险模型的建立。本文的分析包括：选择对生存函数有显著性影响的因子进行建模；通过图形、统计检验的方法比较各个因子变量的不同取值下的生存曲线是否存在差异，差异是否显著；建立的模型是否满足前提假定等。6、数据预处理：

4、由于前312个观测值存在缺失的情况，因此首先进行弥补。之后由于原始数据中status变量分为活着，肝移植和死亡三个结果，而生存分析要求结果为两分类互斥事件，因此将肝移植和死亡的变量值合并。最终status变量的取值为0=活着，1=换肝或死亡，并在312个观测值中将分类变量status、trt、ascites、hepato、spiders、edema、stage转换为因子。经过处理后的数据概况如下：数据集一共包括19个变量，其中11个数值型变量：time、age、bili、chol、albumin、copper、alk.phos、ast、trig、platelet、protime;8个分类变量:

5、 status、trt、sex、ascites、hepato、spiders、edema、stage；无缺失值。2、 非参数模型1、 不含变量的Kaplan-Meier估计得到的95%的置信区间的拟合图如下所示：2、 检验D-青霉胺和安慰剂组生存函数的差异（1）图形检验上图为用Kaplan-Meier 比较D-青霉胺和安慰剂两种不同的药对生存时间的差异，红色的曲线代表安慰剂组，黑色的曲线代表D-青霉胺组，两种不同水平下的曲线大致是重合的，只有小部分存在分歧。因此可以初步判断两组不同的用药对生存函数没有影响。（2） log-rank test（3） Tarone-Ware test（4） Wil

6、coxon test三种检验方法的原假设均为，而得到的p值均远远大于0.05，因此可以认为两种不用的用药对生存函数是没有显著性影响的。3、 检验不同性别的生存函数的差异（1）图形检验（2）log-rank test（3） Tarone-Ware test（4） Wilcoxon test分析：从生存曲线可以看出，在起初的一小段时间里男性和女性的生存曲线重合，没有太大的差异，但随着时间的推移，不同的性别的生存曲线逐渐显示出了较大的差异，从而说明了性别对于生存函数是有影响的，从上面的三种检验得到的p值均为0.02左右也可以得到相应的结论。而且男性患病的风险比女性更高。3、 Cox 比例风险模型1、

7、基于不同用药组（trt)建立Cox比例风险模型从回归系数以及三种检验得到的p值可以看出，不同用药组对于生存函数的影响不显著。2、基于不同性别(sex)建立Cox比例风险模型从图和分析结果均可以看出，不同性别的生存曲线有显著性的差异，女性的风险是男性的0.617倍，男性患pbc疾病的风险更大。之后对sex做PHA的检验，得到的p值是大于0.05的，因此满足风险成比例的假定。对于不同性别组，将建立的KM曲线和Cox比例风险模型的曲线进行比较：(说明：下图中的红色和绿色曲线分别代表的是用Cox比例风险模型对男性和女性的生存状况的拟合，黑色和蓝色分别代表的是用KaplanMeier估计的男性和女性的生

8、存状况）从曲线中可以看出用Cox比例风险模型拟合的女性生存曲线与KM曲线基本重合，男性生存曲线大部分是重合的，小部分存在差异。但从整体来看用Cox比例风险模型拟合的曲线基本与KM曲线重合，因此说明Cox比例风险模型拟合效果较好，也说明了sex变量是满足风险成比例的假定的。3、 基于水肿程度（edema：无水肿，未用利尿剂水肿，用利尿剂仍水肿）建立Cox比例风险模型（1）直接建立Cox比例风险模型从上面的分析结果可以看出，不同水肿程度对生存函数有显著性的影响，其中edema1（用利尿剂仍水肿）的风险是最高的。从医学的角度来看，用利尿剂仍水肿的个体在一定程度上就反映了该个体的生存状况最差，因此其风

9、险理应最高。下图是不同水肿程度拟合的Cox比例风险曲线：检验edema变量是否满足风险比例假定：检验得到的p值均小于0.05，因此认为edema变量是不满足风险比例假定的。（2） 建立edema与时间的交互项模型将heaviside function设置为，建立edema变量与gt的交互项，得到的模型为：得到的交互项的p值为0.00978，小于0.05，因此认为该项是显著的。同时也说明了将heaviside function设置为较为合理。4、基于不同疾病组织学分期(stage)建立Cox比例风险模型从上面的分析结果以及下面的图可以看出，不同疾病组织学分期对生存函数有显著性的影响，相对于其他三

10、个stage=4的风险是最高的。但是从医学的角度来看，疾病组织学分期本身就是根据病情的严重程度来划分的，因此属于stage=4的个体本身已经处于最差的情况当然其风险越高。 从另外一个角度来看，由于有些指标本身就可以反映个体的生存状况，因此指标本身代表的值就可以反映其风险高低程度。而且在观测个体患病时记录的很多变量其实本身就存在高度的相关性。检验stage变量是否满足风险比例假定：检验结果得到的p值都是大于0.05的，因此stage变量是满足风险成比例的。5、利用AIC逐步选择变量：（1）通过10次得到最终的模型，逐步选择的结果如下：逐步选择次数变量AIC1无1477.72bili1394.13

11、3bili、stage1353.824bili、stage、copper1332.95bili、 stage、 copper、 albumin1318.626bili、stage、copper、albumin、protime1312.977bili、stage、copper、albumin、protime、edema1310.758bili、stage、copper、albumin 、protime、edema、 alk.phos1309.399bili、stage、copper、albumin 、protime 、edema 、 alk.phos、 ast1308.210bili、stage

12、、copper、 albumin、 protime、edema、alk.phos、ast、sex1307.67最终得到的模型为：从回归结果的p值可以看出，并不是所有的变量都是显著的，比较显著的有bili、copper、albumin，其余的都不太显著。PHA检验结果：结果显示并不是所有的变量都满足风险成比例的假设，其中bili、protime、edema0.5均不满足。（2）模型改进由于edema0.5既不满足风险成比例，在回归中的系数也不是特别显著，因此直接去掉这个变量，对于bili变量，可以将该变量分段，之后分层建立Cox比例风险模型。分段是按照个变量的均值进行分段，bili的均值为3.2

13、56，因此若，则，若,则。得到的模型为：PHA检验结果为：结果显示protime变量仍不满足风险成比例的假定，去除该变量的拟合结果如下：PHA检验结果如下：所有变量均满足风险成比例的假定。四、参数模型1、没有任何变量的Weibull参数模型：得到的生存函数和风险函数的图如下所示（已将时间转化为年）：2、将拟合的Weibull参数曲线和K-M曲线进行对比发现在开始阶段Weilbull曲线更低，之后的时间里两条曲线基本重合。3、 引入sex变量的Weibull参数曲线：得到的sex1的系数为0.，其对应的p值为0.03，视为显著。得到的kappa值为0.。因此对于女性，其参数模型为：对于男性，其参

14、数模型为：下图为拟合的Weibull参数曲线4、 参数模型的选择选入的变量有stage、copper、albumin、sex，这些变量是在Cox比例风险模型中极为显著的，参数模型的形式选择了4种，最终得到每种模型的AIC值如下：其中最小的是log-logistic模型，输出的模型如下所示：五、代码library(survival)data=pbcw=data1:312,2:20w$statusw$status!=0=1sum(complete.cases(w)library(missForest)w=missForest(w)$ximpwrite.table(w,c:usersthinkDes

15、ktopw.txt)w=read.table(c:usersthinkDesktopw.txt,header=T)w$sex=factor(w$sex);w$trt=factor(w$trt);w$ascites=factor(w$ascites);w$hepato=factor(w$hepato);w$spiders=factor(w$spiders);w$edema=factor(w$edema);w$stage=factor(w$stage)summary(w)fit1=survfit(Surv(time,status)1,data=w)summary(fit1)plot(fit1, m

16、ain=Kaplan-Meier estimate with 95% confidence bounds,xlab=t, ylab=expression(hat(S)(t)fit2=survfit(Surv(time,status)trt,data=w)plot(fit2,=FALSE,mark.time=TRUE,col=c(black,red),lty=1:2,ylab=expression(hat(S)(t),xlab=time,main=Survival Functions)legend(topright,box.lwd=1,cex=0.6,c(D-penicillam

17、ine,placebo),lty=1:2,col=c(black,red)survdiff(Surv(time,status)trt, data=w,rho=0)survdiff(Surv(time,status)trt, data=w,rho=0.5)survdiff(Surv(time,status)trt, data=w,rho=1)fit3=survfit(Surv(time,status)sex,data=w)plot(fit,=FALSE,mark.time=TRUE,col=c(blue,red),lty=1:2,ylab=expression(hat(S)(t)

18、,xlab=time,main=Survival Functions)legend(topright,box.lwd=1,cex=0.6,c(men,female),lty=1:2,col=c(blue,red)survdiff(Surv(time,status)sex, data=w,rho=0)survdiff(Surv(time,status)sex, data=w,rho=0.5)survdiff(Surv(time,status)sex, data=w,rho=1)fit4=coxph(Surv(time,status)trt,data=w)summary(fit4)fit5=cox

19、ph(Surv(time,status)sex,data=w)summary(fit5)d.fit5=coxph.detail(fit5)times=c(0,d.fit5$time)h0=c(0,d.fit5$hazard)s0=exp(-cumsum(h0)beta=c(fit5$coef)xf=c(1)-mean(as.numeric(w$sex)sf=s0exp(t(beta)%*%xf)xm=c(0)-mean(as.numeric(w$sex)sm=s0exp(t(beta)%*%xm)plot(times,sf,type=s,xlab=t,ylab=expression(hat(S

20、)(t),ylim=0:1)lines(times,sm,col=2,type=s)legend(topright,box.lwd=1,cex=0.6,c(female,man),lty=1,col=c(black,red)cox.zph(fit5,transform=rank,global=F)fit6=coxph(Surv(time/365.25,status)sex,data=w)g-levels(w$sex)plot(survfit(fit6,newdata=data.frame(sex=g1),=F,col=2,xlab=t,ylab=expression(hat(S

21、)(t)lines(survfit(fit6, newdata=data.frame(sex=g2),=F,col=3)km=survfit(Surv(w$time/365.25, w$status)w$sex)lines(km,col=c(black,blue),lty=1)legend(topright,box.lwd=1,cex=0.55,c(Cox-man,Cox-female,KM-man,KM-female),lty=1,col=c(red,green,black,blue)fit7=coxph(Surv(time,status)edema,data=w)summa

22、ry(fit7)g-levels(w$edema)plot(survfit(fit7,newdata=data.frame(edema=g1),=F,col=2,xlab=t,ylab=expression(hat(S)(t)lines(survfit(fit7,newdata=data.frame(edema=g2),=F,col=3)lines(survfit(fit7,newdata=data.frame(edema=g3),=F,col=4)legend(topright,c(edema0,edema0.5,edema1),col=c(2

23、:4),lty=1,cex=0.6)cox.zph(fit7,transform=rank,global=F)w2=survSplit(w,cut=c(2000),end=time,event=status,start=start)w2$gt=(w2$start=2000)+0fit8- coxph(Surv(start,time,status)as.numeric(edema)+gt:as.numeric(edema),data=w2)summary(fit8)fit9=coxph(Surv(time,status)stage,data=w)summary(fit9)g-levels(w$s

24、tage)plot(survfit(fit9, newdata=data.frame(stage=g1),=F,col=2,xlab=t,ylab=expression(hat(S)(t)lines(survfit(fit9, newdata=data.frame(stage=g2),=F,col=3)lines(survfit(fit9,newdata=data.frame(stage=g3),=F,col=4)lines(survfit(fit9,newdata=data.frame(stage=g4),=F,col=5)le

25、gend(bottomleft,c(stage1,stage2,stage3,stage4),col=c(2:5),lty=1,cex=0.6)cox.zph(fit9,transform=rank,global=F)library(MASS)Scope=list(upper=(trt+age+sex+ascites+hepato+spiders+edema+bili+chol+albumin+copper+alk.phos+ast+trig+platelet+protime+stage),lower=1)fit10=coxph(Surv(time,status)1,data=w)fit11=

26、stepAIC(fit6,Scope,direction=both)fit12=coxph(Surv(time,status)bili + stage + copper + albumin + protime + edema + alk.phos + ast + sex,data=w)cox.zph(fit12,transform=rank,global=F)c.bili=w$bili3.256c.bili=c.bili+0summary(w)fit9=coxph(Surv(time,status)strata(c.bili) + stage + copper + albumin + alk.

27、phos + ast + sex+protime,data=w)cox.zph(fit9,transform=rank,global=F)fit10=coxph(Surv(time,status)strata(c.bili) + stage + copper + albumin + alk.phos + ast + sex,data=w)cox.zph(fit10,transform=rank,global=F)wei=survreg(Surv(w$time/365.25, w$status)1,dist=w)kappa=wei$scalelambda=exp(-wei$coef1)kappazeit=seq(from=0,to=13,length.out=1000)s=exp(-lambda*zeitkappa)h=lambda*kappa*zeit(kappa-1)par(mfrow=c(1,2)plot(zeit,s,type=l,xlab=t,ylab=expression(hat(s)(t)plot(zeit,h,type=l,xlab=t,ylab=expression(hat(h)