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.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this question
up vote 6 down vote accepted

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...

share|improve this answer
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? – Dihan Dec 26 '12 at 11:11
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. – ocram Dec 26 '12 at 11:32

Your Answer


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.