Apparently, `basehaz` actually computes a cumulative hazard rate, rather than the hazard rate itself. The formula is a 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. #------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