Getting estimates from survfit.coxph
Mark Wardle <mark <at> wardle.org> writes:
I'm having difficulty getting access to data generated by survfit and
print.survfit when they are using with a Cox model (survfit.coxph).
I would like to programmatically access the median survival time for
each strata together with the 95% confidence interval. I can get it on
screen, but can't get to it algorithmically. I found myself examining
the source of print.survfit to try to work out how it is done
internally, but is there a better way?
An example (and I realise that estimating survival curses from an
average status and time is incorrect in this instance, but it keeps
this example simple):
test1 <- list(time= c(4, 3,1,1,2,2,3),
status=c(1,NA,1,0,1,1,0),
x= c(0, 2,1,1,1,0,0),
sex= c(0, 0,0,0,1,1,1))
c1 <- coxph( Surv(time, status) ~ x + strata(sex), test1) #stratified model
f1 <- survfit(c1)
sf1 <- summary(f1)
str(f1)
print(f1)
print(sf1)
str(sf1)
(Disclaimer: there may be a better way got get it with library Design by
Frank Harrell, but let's assume we have to do it the hard way)
Looks like it is a bit hidden. f1 is of class(print.survfit), as str(f1)
tells us. So let's try getAnyhwere(print.survfit). In the lower part you
find line like the following:
x1 <- pfun(nsubjects, stime, surv, x$n.risk, x$n.event,
x$lower, x$upper)
if (is.matrix(x1)) {
if (is.null(x$lower))
dimnames(x1) <- list(NULL, plab)
else dimnames(x1) <- list(NULL, c(plab, plab2))
}
else {
if (is.null(x$lower))
names(x1) <- plab
else names(x1) <- c(plab, plab2)
}
if (show.rmean)
print(x1)
Make a copy of that function under a new name, and return x1.
Dieter