[R] Getting estimates from survfit.coxph

[Mark Wardle]

Dear all,

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)

I'm sure I am missing something obvious. Apologies - but any help
greatfully received!

Best wishes,

Mark

P.S. I can get to diferrent estimates for median survival for
different groups using simpler mechanisms, but they yield different
estimates: From my data, so no reproducible (and ataxSurv() is a
wrapper function that calls plain Surv() after manipulating the data
simply):


# For an "average" patient: (doesn't make any sense biologically)
> survfit(surv.results$cox)
Call: survfit.coxph(object = surv.results$cox)

      n  events  median 0.95LCL 0.95UCL
    136      96       6       6       8

#
# predict a curve for a patient: (these are the answers I really want
to extract)
#
> survfit(surv.results$cox, newdata=data.frame(onset=50, ic.duration=10, simple.msa=c('MSA','Not MSA'), autoimmune=F, carcinoma=F))
Call: survfit.coxph(object = surv.results$cox, newdata = data.frame(onset = 50,
    ic.duration = 10, simple.msa = c("MSA", "Not MSA"), autoimmune = F,
    carcinoma = F))

       n events median 0.95LCL 0.95UCL
[1,] 136     96      8       7      11
[2,] 136     96      3       2       6

#
# without using Cox regression:
#
> survfit(ataxSurv(surv.support, surv.follow.up, surv.results$data) ~ simple.msa, data=surv.results$data)
Call: survfit(formula = ataxSurv(surv.support, surv.follow.up,
surv.results$data) ~
    simple.msa, data = surv.results$data)

   1 observation deleted due to missingness
                     n events median 0.95LCL 0.95UCL
simple.msa=Not MSA 120     80      8       6      11
simple.msa=MSA      19     17      2       1       4


-- 
Dr. Mark Wardle
Specialist registrar, Neurology
Cardiff, UK