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
I suspect the slight difference to be due to the approximation of the partial likelihood in coxph() due to ties in the data...