[R] Getting estimates from survfit.coxph

[Dieter Menne]

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