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

1.0femaleman0.20.6s10246zeit81012 4、参数模型的选择

选入的变量有stage、copper、albumin、sex，这些变量是在Cox比例风险模型中极为显著的，参数模型的形式选择了4种，最终得到每种模型的AIC值如下：

其中最小的是log-logistic模型，输出的模型如下所示：

- 16 -

五、代码

library(survival) data=pbc

w=data[1:312,2:20]

w$status[w$status!=0]=1 sum(complete.cases(w)) library(missForest) w=missForest(w)$ximp

write.table(w,\

w=read.table(\

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, main=\estimate with 95% confidence bounds\ylab=expression(hat(S)(t)))

fit2=survfit(Surv(time,status)~trt,data=w)

plot(fit2,conf.int=FALSE,mark.time=TRUE,col=c(\ab=\

legend(\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,conf.int=FALSE,mark.time=TRUE,col=c(\b=\

legend(\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=coxph(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)

- 17 -

s0=exp(-cumsum(h0)) beta=c(fit5$coef)

xf=c(1)-mean(as.numeric(w$sex)) sf=s0^exp(t(beta)%*%xf)

xm=c(0)-mean(as.numeric(w$sex)) sm=s0^exp(t(beta)%*%xm)

plot(times,sf,type='s',xlab=\lines(times,sm,col=2,type='s')

legend(\cox.zph(fit5,transform=\

fit6=coxph(Surv(time/365.25,status)~sex,data=w) g<-levels(w$sex)

plot(survfit(fit6,newdata=data.frame(sex=g[1])),conf.int=F,col=2,xlab=\)))

lines(survfit(fit6, newdata=data.frame(sex=g[2])),conf.int=F,col=3) km=survfit(Surv(w$time/365.25, w$status)~w$sex) lines(km,col=c(\

legend(\col=c(\

fit7=coxph(Surv(time,status)~edema,data=w) summary(fit7)

g<-levels(w$edema)

plot(survfit(fit7,newdata=data.frame(edema=g[1])),conf.int=F,col=2,xlab=\S)(t)))

lines(survfit(fit7,newdata=data.frame(edema=g[2])),conf.int=F,col=3) lines(survfit(fit7,newdata=data.frame(edema=g[3])),conf.int=F,col=4)

legend(\cox.zph(fit7,transform=\

w2=survSplit(w,cut=c(2000),end=\w2$gt=(w2$start==2000)+0

fit8<- 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$stage) plot(survfit(fit9,

newdata=data.frame(stage=g[1])),conf.int=F,col=2,xlab=\lines(survfit(fit9, newdata=data.frame(stage=g[2])),conf.int=F,col=3) lines(survfit(fit9,newdata=data.frame(stage=g[3])),conf.int=F,col=4) lines(survfit(fit9,newdata=data.frame(stage=g[4])),conf.int=F,col=5)

- 18 -

legend(\cox.zph(fit9,transform=\

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=stepAIC(fit6,Scope,direction=\

fit12=coxph(Surv(time,status)~bili + stage + copper + albumin + protime + edema + alk.phos + ast + sex,data=w)

cox.zph(fit12,transform=\

c.bili=w$bili<3.256 c.bili=c.bili+0 summary(w)

fit9=coxph(Surv(time,status)~strata(c.bili) + stage + copper + albumin + alk.phos + ast + sex+protime,data=w)

cox.zph(fit9,transform=\

fit10=coxph(Surv(time,status)~strata(c.bili) + stage + copper + albumin + alk.phos + ast + sex,data=w)

cox.zph(fit10,transform=\

wei=survreg(Surv(w$time/365.25, w$status)~1,dist=\kappa=wei$scale

lambda=exp(-wei$coef[1])^kappa

zeit=seq(from=0,to=13,length.out=1000) s=exp(-lambda*zeit^kappa)

h=lambda*kappa*zeit^(kappa-1)

par(mfrow=c(1,2))

plot(zeit,s,type=\plot(zeit,h,type=\

plot(survfit(Surv(w$time/365.25,

w$status)~1),xlab='t',ylab=expression(hat(s)(t)),lty=2,conf.int=F,col=3) lines(zeit,s,col=2)

legend(\

wei=survreg(Surv(w$time/365.25, w$status)~w$sex,dist=\kappa=wei$scale

lambda1=exp(-wei$coeff[1])^kappa

lambda2=exp(-wei$coeff[1]-wei$coeff[2])^kappa

- 19 -

zeit=seq(from=0,to=13,length.out=1000) s1=exp(-lambda1*zeit^kappa) s2=exp(-lambda2*zeit^kappa)

plot(zeit,s1,conf.int=F,col=1,type=\lines(zeit,s2,col=2,lwd=2)

legend(\

wei=survreg(Surv(time,status)~stage + copper + albumin + sex,data=w,dist=\

exp=survreg(Surv(time,status)~stage + copper + albumin + sex,data=w,dist=\lgl=survreg(Surv(time,status)~stage + copper + albumin + sex,data=w,dist=\lgn=survreg(Surv(time,status)~stage + copper + albumin + sex,data=w,dist=\extractAIC(wei)[2] extractAIC(exp)[2] extractAIC(lgl)[2] extractAIC(lgn)[2] lgl

- 20 -