library(survival) attach(aml) aml.surv <- Surv(aml$time,aml$status) aml1.surv <- Surv(aml$time[aml$x=="Maintained"], aml$status[aml$x=="Maintained"]) aml0.surv <- Surv(aml$time[aml$x=="Nonmaintained"],aml$status[aml$x=="Nonmaintained"]) group <- ifelse(aml$x=="Nonmaintained",rep(0,23),rep(1,23)) aml1.surv aml0.surv summary(coxph(aml.surv~group)) summary(survfit(coxph(aml.surv~ group),data.frame(group=0))) summary(survfit(coxph(aml.surv~ group))) #---- Motorette example --- library(MASS) motorette <- motors[11:40,] attach(motorette) Surv(time[temp==170],cens[temp==170]) Surv(time[temp==190],cens[temp==190]) Surv(time[temp==220],cens[temp==220]) x <- 1000 / (273.2+motorette$temp) motorette.surv <- Surv(motorette$time,motorette$cens) summary(coxph(motorette.surv~x)) par(mfcol=c(1,2)) plot(survfit(motorette.surv ~ x,type="kaplan-meier"),conf.int=F,main="KM fit") plot(survfit(coxph(motorette.surv ~ x),data.frame(x=unique(x)[1])),conf.int=F,main="Cox PH model") lines(survfit(coxph(motorette.surv ~ x),data.frame(x=unique(x)[2])),conf.int=F) lines(survfit(coxph(motorette.surv ~ x),data.frame(x=unique(x)[3])),conf.int=F) #-------------- martingale and deviance residuals ---- par(mfrow=c(1,3)) scatter.smooth(1:30,fit.x$res,xlab="index",ylab="residuals",main="Martingale residuals") abline(0,0) scatter.smooth(1:30,resid(fit.x,type="deviance"),xlab="index",ylab="residuals",main="Deviance residuals") abline(0,0) scatter.smooth(x,resid(fit.x,type="deviance"),xlab="X",ylab="residuals",main="Deviance residuals") abline(0,0) index <- 1:30 summary(coxph(motorette.surv~ x + index)) #---------------- Schoenfeld res and dfbetas to examine fit and outlying covariate values ---- fit.temp <- coxph(motorette.surv ~ motorette$temp) par(mfcol=c(1,2)) detail <- coxph.detail(fit.x) tt <- detail$y[,2] ss <- detail$y[,3] sch <- resid(fit.temp,type="schoenfeld") sch <- resid(fit.x,type="schoenfeld") plot(tt[ss==1],sch,xlab="ordered survival times",ylab="schoenfeld res for temp") ind <- as.numeric(dimnames(detail$y)[[1]]) plot(ind[ss==1],sch,xlab="Index",ylab="schoenfeld res for temp") par(mfcol=c(1,2)) plot(motorette$temp,motorette$time) plot(motorette$temp,motorette$time,log="y") motorette par(mfcol=c(1,1)) bresid <- resid(fit.x,type="dfbeta") plot(11:40,bresid,type="h",ylab="scales change in coeff.",xlab="observation") summary(coxph(motorette.surv~x)) summary(coxph(motorette.surv~x,subset=((1:30) != 11)))