From: Daniel A. Powers (dpowers@mail.la.utexas.edu)
Date: Fri 23 Feb 2001 - 13:53:15 EST
Message-ID: <Pine.GSO.4.05.10102222150230.7015-100000@dahlia.la.utexas.edu>
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@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
This archive was generated by hypermail 2b30 : Fri 22 Jun 2001 - 18:58:33 EST