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 + survest
age 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.
[EMAIL PROTECTED]
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
