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)
#
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
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