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


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.

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)
  • $\begingroup$ You are right, something seems strange. You might have to dig into the 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, 2013 at 15:06

3 Answers 3


After looking for the source code of basehaz.S, I've got the reason why I am wrong here. First basehaz simply compute cumulate hazard function $\Lambda_0(t)$ by using survfit instead of instantaneous hazard function $\lambda_0(t)$

Second the main code for basehaz.S is

    H<- -log(sfit$surv)

Then it's wrong to use this function to estimate a time dependent covariates Cox model (which require id option in survfit).

I think we should avoid to use the basehaz function becuase it exists only because Prof. Therneau try to comfort SAS programmers as he described in the document of basehaz function.


They give the same results. Try this:

    H0 <- basehaz(fit.cox, centered=T)
    h0 <- unique(-log(survfit(fit.cox)$surv))

You will see that the two results are the same. If you want to get the survival instead of the cumulative hazard:

    S0 <- exp(-H0)
    s0 <- survfit(fit.cox)$surv

I found that doing exp(-H0) where H0=basehaz(fit.cox, centered=T) also changed the time values.

a solution to this was:

H0=basehaz(fit.cox, centered=T)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.