Xianqun (Wilson) Wang
2007-Nov-21 22:55 UTC
[R] survest and survfit.coxph returned different confidence intervals on estimation of survival probability at 5 year
I wonder if anyone know why survest (a function in Design package) and standard survfit.coxph (survival) returned different confidence intervals on survival probability estimation (say 5 year). I am trying to estimate the 5-year survival probability on a continuous predictor (e.g. Age in this case). Here is what I did based on an example in "help cph". The 95% confidence intervals returned by survfit.coxph and survest seem different. ### Create example dataset n <- 1000 set.seed(731) age <- 50 + 12*rnorm(n) label(age) <- "Age" sex <- factor(sample(c('Male','Female'), n, TRUE)) cens <- 15*runif(n) h <- .02*exp(.04*(age-50)+.8*(sex=='Female')) dt <- -log(runif(n))/h label(dt) <- 'Follow-up Time' e <- ifelse(dt <= cens,1,0) dt <- pmin(dt, cens) units(dt) <- "Year" dd <- datadist(age, sex) options(datadist='dd') Srv <- Surv(dt,e) ### Use Design package cph and survest to estimate 5-year survival and its 95% CI. library(Design) par(mfrow=c(1,2), pty="s") fd <- cph(Srv ~ age, x=TRUE, y=TRUE, surv=T) tmp <- survest(fd, newdata= age, times=5) dd.plot <- data.frame(age=age, lower=tmp$lower[,1], surv=tmp$surv[,1], upper=tmp$upper[,1]) oo <- order(age) dd.plot <- dd.plot[oo,] matplot(dd.plot$age, dd.plot[,-1], col=c(2,1,2), type="l", lty=c(4,1,4), xlab="Age", ylab="Survival Prob", main="Use cph, survest") ### Use coxph and survfit to estimate 5-year survival and its 95% CI. detach(package:Design) fs <- coxph(Srv ~ age, x=T, y=T) tmp <- survfit(fs, newdata= data.frame(age=age)) tm <- max((1:length(tmp$time))[tmp$time <= 5 + 1e-6]) ds.plot <- data.frame(age=age, lower=tmp$lower[tm, ], surv=tmp$surv[tm, ], upper=tmp$upper[tm, ]) oo <- order(age) ds.plot <- ds.plot[oo,] matplot(ds.plot$age, ds.plot[,-1], col=c(2,1,2), type="l", lty=c(4,1,4), xlab="Age", ylab="Survival Prob", main="Use coxph, survfit") #### The estimated 95% CI are different between survfit.coxph and survest> ii=(1:nrow(dd.plot))[(dd.plot$age>58 & dd.plot$age<59)] > dd.plot[ii,] ### Estimation by cph + survestage lower surv upper 95 58.18543 0.8098515 0.8144546 0.8189590 719 58.23063 0.8094270 0.8141218 0.8187143 573 58.34560 0.8083436 0.8132730 0.8180904 763 58.37432 0.8080720 0.8130605 0.8179343 33 58.51589 0.8067287 0.8120095 0.8171628 560 58.54394 0.8064616 0.8118006 0.8170096 11 58.56214 0.8062881 0.8116650 0.8169102 269 58.56323 0.8062777 0.8116569 0.8169042 993 58.60200 0.8059077 0.8113677 0.8166922 506 58.67991 0.8051621 0.8107854 0.8162654 433 58.72130 0.8047651 0.8104754 0.8160384 583 58.72179 0.8047604 0.8104717 0.8160357 958 58.72254 0.8047532 0.8104661 0.8160316 125 58.73593 0.8046245 0.8103657 0.8159581 42 58.76653 0.8043304 0.8101361 0.8157900 992 58.78134 0.8041878 0.8100249 0.8157086 988 58.85797 0.8034489 0.8094486 0.8152869 618 58.93011 0.8027511 0.8089047 0.8148891 199 58.94150 0.8026407 0.8088186 0.8148262 475 58.97280 0.8023370 0.8085821 0.8146533 ### Estimation by coxph + survfit> ds.plot[ii,]age lower surv upper 95 58.18543 0.7844845 0.8144546 0.8455696 719 58.23063 0.7840935 0.8141218 0.8453001 573 58.34560 0.7830952 0.8132730 0.8446138 763 58.37432 0.7828450 0.8130605 0.8444422 33 58.51589 0.7816067 0.8120095 0.8435949 560 58.54394 0.7813603 0.8118006 0.8434268 11 58.56214 0.7812004 0.8116650 0.8433177 269 58.56323 0.7811908 0.8116569 0.8433112 993 58.60200 0.7808496 0.8113677 0.8430786 506 58.67991 0.7801619 0.8107854 0.8426109 433 58.72130 0.7797957 0.8104754 0.8423622 583 58.72179 0.7797913 0.8104717 0.8423592 958 58.72254 0.7797847 0.8104661 0.8423547 125 58.73593 0.7796660 0.8103657 0.8422742 42 58.76653 0.7793946 0.8101361 0.8420902 992 58.78134 0.7792630 0.8100249 0.8420011 988 58.85797 0.7785811 0.8094486 0.8415397 618 58.93011 0.7779371 0.8089047 0.8411050 199 58.94150 0.7778352 0.8088186 0.8410363 475 58.97280 0.7775549 0.8085821 0.8408474 Wilson Wang AviaraDx Inc. xwang@aviaradx.com [[alternative HTML version deleted]]