Hello everybody. First of all, I would like to present myself. I'm a french student in public health and I like statistics though I'm not that good in mathematics (but I try to catch up). I've discovered R recently while trying to find a statistical program in order to avoid rebooting my computer under windows when I need to do some statistical work. And here is my first question. I'm actually trying to model a survival curve using a cox model. But I need the baseline cumulative hazard. I've been looking around on the web and the archives, and I found that it as to do with the nelson-aalen estimation. But I must not understand all the doc very well, because I can't get the baseline hazard. Could you help me, please. Best regards. PS : I hope you'll forgive my poor english writing. -- Blaise TRAMIER aka Meles. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Thu, 22 Feb 2001, Meles MELES wrote:> Hello everybody. > First of all, I would like to present myself. > I'm a french student in public health and I like statistics though I'm > not that good in mathematics (but I try to catch up). I've discovered R > recently while trying to find a statistical program in order to avoid > rebooting my computer under windows when I need to do some statistical > work. > > And here is my first question. > I'm actually trying to model a survival curve using a cox model. > But I need the baseline cumulative hazard. I've been looking around on > the web and the archives, and I found that it as to do with the > nelson-aalen estimation. > But I must not understand all the doc very well, because I can't get > the baseline hazard. > Could you help me, please.In the survival5 package the function survfit() will extract baseline survival from a Cox model. In order to get the baseline hazard you need to transform the baseline survival. There are two methods of calculating the baseline survival, the default one gives the baseline hazard estimator you want. It is attributed to Aalen, Breslow, or Peto, depending on who one asks. Using the example from help(survfit) #fit a cox proportional hazards model and plot the #predicted survival curve data(ovarian) fit <- coxph( Surv(futime,fustat)~resid.ds+rx+ecog.ps,data=ovarian) Compute the survival function ss<-survfit(fit) Now you can graph the baseline cumulative hazard with plot(ss,fun="cumhaz") or extract the information into a variable with cumulative.hazard<- -log(ss$surv) observation.times<-ss$time The survival object "ss" also contains upper and lower confidence bounds and other information about the survival curve -thomas Thomas Lumley Asst. Professor, Biostatistics tlumley at u.washington.edu University of Washington, Seattle -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Meles MELES <meles at free.fr> writes:> And here is my first question. > I'm actually trying to model a survival curve using a cox model. > But I need the baseline cumulative hazard. I've been looking around on > the web and the archives, and I found that it as to do with the > nelson-aalen estimation. > But I must not understand all the doc very well, because I can't get > the baseline hazard. > Could you help me, please.survfit(coxph(...)) gives you the baseline survival function, which by default (type="aalen") is calculated as S(t) = exp(-Lambda(t)) where Lambda(t) is the cumulative hazard. So all you have to do is to reverse the formula as Lambda(t) = -log(S(t)) The estimated S(t) is in the $surv component of the object returned by survfit(). e.g. example(survfit) s <- survfit(fit) plot(c(0,s$time),c(0,-log(s$surv)),type='s') -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Meles --
Here are a couple of routines. One computes Aalen's estimate of the
cumulative hazard and the other a log-likelihood based on it
Cheers,
Dan
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-Dan Powers
Associate Professor, Sociology
University of Texas at Austin
--snip
aalen <- function(coxfit)
{
#calculate Aalen's estimate of cum hazard
if(!inherits(coxfit,"coxph")) stop(
"argument must be a proportional hazard's fit")
event <- coxfit$y[,"status"] == 1
time <- coxfit$y[,"time"]
risk <- exp(coxfit$linear.predictor)
dt <- unique(time[event])
k <- length(dt)
lambda <- rep(0,k)
for(i in 1:k) {
lambda[i] <- sum(event[time==dt[i]])/sum(risk[time >= dt[i]])
}
data.frame(time=dt, lambda=lambda)
}
aalen.loglik <- function(coxfit)
{
# Calculate survival likelihood using Cox-Aalen Estimates
#
a <- aalen(coxfit)
event <- coxfit$y[,"status"] == 1
time <- coxfit$y[,"time"]
risk <- exp(coxfit$linear.predictors)
loglik <-0
for(i in 1:length(event)) {
if (event[i]) {
lambda <- risk[i] * a$lambda[a$time == time[i]]
loglik <- loglik + log(lambda)
}
loglik <- loglik - risk[i] * sum(a$lambda[a$time <= time[i]])
}
loglik
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._