Let's say I have a "kidney catheter" data set. I'm trying to model a survival curve using a Cox model. If I consider a Cox model: $$h(t,Z) = h_0 \exp(b'Z),$$ I need the estimate of the baseline hazard. By using the built-in survival package R function basehaz(), I can easily do it like this:


fit <- coxph(Surv(time, status) ~ age , kidney)

But if I want to write a step by step function of the baseline hazard for a given estimate of parameter b how can I proceed? I tried:

bhaz <- function(beta, time, status, x) {

    data <- data.frame(time,status,x)
    data <- data[order(data$time), ]
    dt   <- data$time
    k    <- length(dt)
    risk <- exp(data.matrix(data[,-c(1:2)]) %*% beta)
    h    <- rep(0,k)

    for(i in 1:k) {
        h[i] <- data$status[data$time==dt[i]] / sum(risk[data$time>=dt[i]])          

    return(data.frame(h, dt))

h0 <- bhaz(fit$coef, kidney$time, kidney$status, kidney$age)

But this does not give the same result as basehaz(fit). What is the problem?

  • $\begingroup$ @gung could you help with this question? I struggled for couple of days... $\endgroup$
    – Haitao Du
    Commented Mar 22, 2017 at 20:46

1 Answer 1


Apparently, basehaz() actually computes a cumulative hazard rate, rather than the hazard rate itself. The formula is as follows: $$ \hat{H}_0(t) = \sum_{y_{(l)} \leq t} \hat{h}_0(y_{(l)}), $$ with $$ \hat{h}_0(y_{(l)}) = \frac{d_{(l)}}{\sum_{j \in R(y_{(l)})} \exp(\mathbf{x}^{\prime}_j \mathbf{\beta})} $$ where $y_{(1)} < y_{(2)} < \cdots$ denote the distinct event times, $d_{(l)}$ is the number of events at $y_{(l)}$, and $R(y_{(l)})$ is the risk set at $y_{(l)}$ containing all individuals still susceptible to the event at $y_{(l)}$.

Let's try this. (The following code is there for illustration only and is not intended to be very well written.)


#------some data------

tab <- data.frame(table(kidney[kidney$status == 1, "time"])) 
y <- as.numeric(levels(tab[, 1]))[tab[, 1]] #ordered distinct event times
d <- tab[, 2]                               #number of events

#------Cox model------
fit<-coxph(Surv(time, status)~age, data=kidney)

#------cumulative hazard obtained from basehaz()------
H0 <- basehaz(fit, centered=FALSE)
H0 <- H0[H0[, 2] %in% y, ] #only keep rows where events occurred

#------my quick implementation------
betaHat <- fit$coef

h0 <- rep(NA, length(y))
for(l in 1:length(y))
  h0[l] <- d[l] / sum(exp(kidney[kidney$time >= y[l], "age"] * betaHat))

cbind(H0, cumsum(h0))

partial output:

       hazard time cumsum(h0)
1  0.01074980    2 0.01074980
5  0.03399089    7 0.03382306
6  0.05790570    8 0.05757756
7  0.07048941    9 0.07016127
8  0.09625105   12 0.09573508
9  0.10941921   13 0.10890324
10 0.13691424   15 0.13616338

I suspect that the slight difference might be due to the approximation of the partial likelihood in coxph() due to ties in the data...

  • $\begingroup$ Thanks a lot. Yes, there are slight difference for approximation method. But there are 76 time points with ties, if I want to find the baseline hazard for every time point. What can i do? What type of modification in R code is needed? $\endgroup$
    – Dihan
    Commented Dec 26, 2012 at 11:11
  • 1
    $\begingroup$ The discretised hazard is zero, except at event times. This indeed gives the largest contribution to the likelihood if a discrete hazard function is supposed. You might want to interpolate between any two estimates assuming, for example, that the hazard stays constant. $\endgroup$
    – ocram
    Commented Dec 26, 2012 at 11:32
  • $\begingroup$ Method of Breslow (1974) $\endgroup$
    – tomka
    Commented Feb 27, 2017 at 17:49
  • $\begingroup$ I need to note some problems with this implementation. Using kidney$time >= y[l] can run into numerical problems when time is numeric due to the tabulation in creating $y$. Furthermore, the way you define your risk set is inaccurate, because if there is a tie of two observations, one with status=0 and one with status=1, then $d=2$ but your code gives $d=1$ as you exclude all status=0 observations. The latter problem applies higher numbers of ties likewise. $\endgroup$
    – tomka
    Commented Feb 28, 2017 at 18:10
  • 1
    $\begingroup$ As \@ocram defined the baseline hazard, you can find this formula in the book Applied Survival Analysis Using R-2016 by Dirk F. Moore (page 64), either. I applied the formula for my dataset and it works, in fact, it estimated the baseline hazrad correctly as I checked it couple of times for some different covariates. And also, the slight difference in the results because of method is used. As @mr.bjerre stated, we will get the same results when we choose the method='breslow' in fitting cox proportional hazard. $\endgroup$
    – mhy
    Commented Mar 21, 2020 at 13:07

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.