Cox Baseline Hazard

Question

Let's say I have a 'kidney catheter' data set. I'm actually trying to model a survival curve using a Cox model. If I consider a Cox model

                                h(t,Z)=h0(t)exp(b'Z)

I need the estimate baseline hazard. By using built in survival package R code basehaz(), I can easily do it like this

R code

library(survival)
data(kidney)
fit<-coxph(Surv(time, status) ~ age , kidney)
           basehaz(fit)

But if I want to write a step by step function of baseline hazard for given estimate of parameter "b", how can I proceed? I am trying like,

R code

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, it does not give the same result like basehaz(fit). What is the problem?

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

#------package------
library(survival)
#-------------------

#------some data------
data(kidney)
#---------------------

#------preparation------
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))
}
#-----------------------------------

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