Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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)
share|improve this question
    
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. –  Aniko Aug 22 '13 at 15:06

1 Answer 1

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

    sfit<-survfit(fit)
    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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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