I am trying to caculate survival function in a time dependent covariates Cox model from its baseline hazard function. However, my program gives a different result compared with survfit
. Based on the formula
$S(t)=\exp(-\int_{0}^{t}\lambda_{0}(\mu)\exp[\hat{\beta}Z(\mu)]d\mu$
result should match. I just need this program for illustration but speed. The data I used can be download here and program have been attached. Figure 1 shown the graph generated by the program below.
library(survival)
fit.cox <- coxph( Surv(t1,t2,event) ~ x, data = sim$data)
lambda0 <- basehaz(fit.cox, centered = F) # Estimated Baseline Hazard Function
# Caculate Survival Function
t.cut <- sim$t.cut
lambda0 <- rbind(lambda0, c(0,0) )
x <- sim$x[1, ]
beta <- fit.cox$coefficients
pred <- exp( x * beta )
#Baseline Hazard Function
lambda0.fun <- function(t) {
approx(x = lambda0$time, y = lambda0$hazard, xout = t)$y
}
#Hazard Function
lambda.fun <- function(t){
which.t <- sum(t >= t.cut)
lambda0.fun(t) * pred[which.t]
}
#Survival Function
survival.fun <- function(t){
cum.hazard <- integrate(Vectorize(lambda.fun), 0, t, subdivisions = 1e3L)
exp(- cum.hazard$value)
}
#Test Program
test <-sim$data[1,]
s.est <- Vectorize(survival.fun)(seq(0,2,0.1))
s.est2 <- survfit(fit.cox, newdata = test, id = id,
se.fit = F, type = "efron")
plot(s.est2, xlim = c(0,2))
points(seq(0,2,0.1), s.est)
survival:::survfit.coxph
function, or contact the package maintainer. A few minor quibbles: you should use piecewise constant interpolation of the hazard function, and summation is probably better than integration, but these changes won't fix your main problem. $\endgroup$ – Aniko Aug 22 '13 at 15:06