Hi,
I fit a Cox PH model to estimate the cause-specific hazards (in a competing
risks setting). Then , I compute the survival estimates for all the
individuals in my data set using the `survfit' function. I am currently
playing with a data set that has about 6000 observations and 12 covariates. I
am finding that the survfit function is very slow.
Here is a simple simulation example (modified from Frank Harrell's example for
`cph') that illustrates the problem:
#n <- 500
set.seed(4321)
age <- 50 + 12*rnorm(n)
sex <- factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4)))
cens <- 5 * runif(n)
h <- 0.02 * exp(0.04 * (age-50) + 0.8 * (sex=='Female'))
dt <- -log(runif(n))/h
e <- ifelse(dt <= cens, 1, 0)
dt <- pmin(dt, cens)
Srv <- Surv(dt, e)
f <- coxph(Srv ~ age + sex, x=TRUE, y=TRUE)
system.time(ans <- survfit(f, type="aalen", se.fit=FALSE, newdata=f$x))
When I run the above code with sample sizes, n, taking on values of 500, 1000,
2000, and 4000, the time it takes for survfit to run are as follows:
# n <- 500
> system.time(ans <- survfit(f, type="aalen", se.fit=FALSE, newdata=f$x))
user system elapsed
0.16 0.00 0.15
# n <- 1000
> system.time(ans <- survfit(f, type="aalen", se.fit=FALSE, newdata=f$x))
user system elapsed
1.45 0.00 1.48
# n <- 2000
> system.time(ans <- survfit(f, type="aalen", se.fit=FALSE, newdata=f$x))
user system elapsed
10.19 0.00 10.25
# n <- 4000
> system.time(ans <- survfit(f, type="aalen", se.fit=FALSE, newdata=f$x))
user system elapsed
72.87 0.05 74.87
I eventually want to use `survfit' on a data set with roughly 50K observations,
which I am afraid is going to be painfully slow. I would much appreciate
hints/suggestions on how to make `survfit' faster or any other faster
alternatives.
Thanks.
Best,
Ravi.
____________________________________________________________________
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
______________________________________________
R-help@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.
