R help - Nov 2012 - Fitting and plotting a coxph with survfit, package(surv)

Hi Dear R-users

I have a database with 18000 observations and 20 variables. I am running
cox regression on five variables and trying to use survfit to plot the
survival based on a specific variable without success.

Lets say I have the following coxph:
>library(survival)
>fit <- coxph(Surv(futime, fustat) ~ age + rx, data = ovarian)
>fit
what I am trying to do is plot a survival comparing objects based on rx.
Using this
>plot(survfit(fit, newdata=data.frame(rx =c(1:2), age=c(60)),
             xscale=365.25, xlab = "Years",
ylab="Survival"))
I get the survival for patients at 60, but is there an option to get a
survfit for the patients regardless of the value in variable "age"?

Thanks in advance
Maziar Mohaddes
M.D.
Gothenburg, Sweden

	[[alternative HTML version deleted]]

Your model is additive so the effect of rx is the same at every age.  There is
not one survival curve for all ages (unless the beta for age is 0).  The curves
will shift up and down as you vary age, but they will retain the same relation. 
A common approach is to use the sample mean of age.  Alternatively you can
compute the survival curves (or the cumulative hazard curves) for many ages and
compute a weighted average of them (weights depending on some reference
population perhaps). Finally, if you leave age out of the model you don't
even need to specify age to draw the survival curves.  It kind of depends on
what you are trying to do.

Chris

Ps last time I posted code it all ran together.  We'll see what happens this
time.

library(survival)
data(kidney)

mod1 <- coxph(Surv(time, status) ~ sex + age, data=kidney)
mod2 <- coxph(Surv(time, status) ~ sex, data=kidney)
meanage <- mean(kidney$age)

plot(survfit(mod1, newdata=data.frame(sex =c(1:2), age=meanage)), xscale=365.25,
xlab = "Years", ylab="Survival", col=1:2)
plot(survfit(mod2, newdata=data.frame(sex =c(1:2))), xscale=365.25, xlab =
"Years", ylab="Survival", col=1:2)


ttt <- hist(kidney$age)
weights <- ttt$counts/ sum(ttt$counts)

survcurves <- survfit(mod, newdata=expand.grid(sex =c(1:2), age=ttt$mids),
se.fit=FALSE)
str(survcurves)
sc2 <- survcurves
sc2$surv <- t(apply(survcurves$surv, 1, function(x) 
c(mean(x[seq(1,length(x),2)]), mean(x[seq(2,length(x),2)]))))
plot(sc2, xscale=365.25, xlab = "Years", ylab="Survival",
col=1:2)


-----Original Message-----
From: Maziar Mohaddes [mailto:maziar.mohaddes at gmail.com] 
Sent: Tuesday, November 27, 2012 12:45 PM
To: r-help at r-project.org
Subject: [R] Fitting and plotting a coxph with survfit, package(surv)

Hi Dear R-users

I have a database with 18000 observations and 20 variables. I am running cox
regression on five variables and trying to use survfit to plot the survival
based on a specific variable without success.

Lets say I have the following coxph:
>library(survival)
>fit <- coxph(Surv(futime, fustat) ~ age + rx, data = ovarian) fit
what I am trying to do is plot a survival comparing objects based on rx.
Using this
>plot(survfit(fit, newdata=data.frame(rx =c(1:2), age=c(60)),
             xscale=365.25, xlab = "Years",
ylab="Survival")) I get the survival for patients at 60, but is there
an option to get a survfit for the patients regardless of the value in variable
"age"?

Thanks in advance
Maziar Mohaddes
M.D.
Gothenburg, Sweden

	[[alternative HTML version deleted]]


**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be used
for urgent or sensitive issues

I answered a similar question yesterday:
"The survfit routine will produce predicted survival curves for any
requested combination
of the covariates in the original model.  This is not the same thing as an
"adjusted"
survival curve.  Confusion on this is prevalent, however.  True adjustment
requires a
population average over the confounding factors and is closely related to the
standardized
incidence ratio concept found in epidemiology."

To do this you need to define a poplation of ages.  See chapter 10 of the book
by Therneau
and Grambsch for an explantion of the issues and examples of how to get the
population
value.  It's hard to distill 20 pages down into an email message.

Terry Therneau

---------- begin included message ---------

I have a database with 18000 observations and 20 variables. I am running
cox regression on five variables and trying to use survfit to plot the
survival based on a specific variable without success.

Lets say I have the following coxph:
 >library(survival)
 >fit <- coxph(Surv(futime, fustat) ~ age + rx, data = ovarian)
 >fit
what I am trying to do is plot a survival comparing objects based on rx.
Using this
 >plot(survfit(fit, newdata=data.frame(rx =c(1:2), age=c(60)),
              xscale=365.25, xlab = "Years",
ylab="Survival"))
I get the survival for patients at 60, but is there an option to get a
survfit for the patients regardless of the value in variable "age"?

The nice thing with R (in contrary  to point and click statistical
software) and this community is that you learn alot. Well u are forced to
in order ro be able to ask the question :-)
I am not only refering to codes provided helping me in the analysis but
also pure statistical learning.
What I am trying to plot is a true adjusted survival curve and I realise
that survfit from coxph prob is not the right way to go.

Thanks Chris for the code and thank you Terry. Will google to find the book
and read the chapter you suggested.

	[[alternative HTML version deleted]]

On Nov 28, 2012, at 2:49 PM, Maziar Mohaddes wrote:
> The nice thing with R (in contrary  to point and click statistical
> software) and this community is that you learn alot. Well u are forced to
> in order ro be able to ask the question :-)
> I am not only refering to codes provided helping me in the analysis but
> also pure statistical learning.
Yes. That is very true. The statisticians using R do not generally think that
all the methods promulgated by the SAS and SPSS manuals are correct. I find it a
challenge to keep the boundary between what I do know and do not know clear to
myself. From time to time it becomes clear that I have stepped over than line
and I generally hear about it quickly. This I think is good thing.
> What I am trying to plot is a true adjusted survival curve and I realise
> that survfit from coxph prob is not the right way to go.
There was an article a couple of years ago in JAMA from a group in Canada
advocating for one version of a "true adjusted survival curve" rather
than what might be called the naive "survival curve of the means." I
don't think it matters much, myself.
> 
> Thanks Chris for the code and thank you Terry. Will google to find the book
> and read the chapter you suggested.
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> 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.
David Winsemius, MD
Alameda, CA, USA