[R] Getting estimates from survfit.coxph

[Frank E Harrell Jr]

Dieter Menne wrote:
> 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)

Aside from getting confidence limits, the Quantile function in Design 
(Quantile.cph) will create an R function with the analytic 
representation of the quantiles of the fitted survival distribution.

Frank

> 
> 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
> 
> ______________________________________________
> R-help at 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.
> 


-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University