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