[R] [newbie] Cox Baseline Hazard

[Daniel A. Powers]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._