# How to estimate baseline hazard function in Cox Model with R

I need to estimate baseline hazard function $\lambda_0(t)$ in a time dependent Cox model

$\lambda(t) = \lambda_0(t) \exp(Z(t)'\beta)$

While I took Survival course, I remember that the direct derivative of cumulative hazard function ($\lambda_0(t) dt = d\Lambda_0(t)$) would not be a good estimator because Breslow estimator gives a step function.

So, is there any function in R that I could use directly ? Or any reference on this topic ?

I am not sure if it is worth to open another question, so I just add some background why baseline hazard function is important for me. The formula below estimates the probability that the survival time for one subject is larger than another,. Under a Cox model setting, baseline hazard function $\lambda_0(t)$ is required.

$P(T_1 > T_2 ) = - \int_0^\infty S_1(t) dS_2(t) = - \int_0^\infty S_1(t)S_2(t)\lambda_2(t)dt$

A Cox model was explicitly designed to be able to estimate the hazard ratios without having to estimate the baseline hazard function. This is a strength and a weakness. The strength is that you cannot make errors in functions you don't estimate. This is a real strength and is the reason why people refer to it as "semi-parametric" and is to a large extent responsible for its popularity. However, it is also a real weakness, in that once you want to know something other than the hazard ratio, you will often require the baseline hazard function and that defeats the very purpose of a Cox model.

So I tend to use Cox models only when I am interested in hazard ratios and nothing else. If I want to know other things, I typically move on to other models like the ones discussed here: http://www.stata.com/bookstore/flexible-parametric-survival-analysis-stata/

The baseline hazard function can be estimated in R using the "basehaz" function. The "help" file states that it is the "predicted survival" function which it's clearly not. If one inspects the code, it's clearly the cumulative hazard function from a survfit object. For further silliness, the default setting is centered=TRUE which a) is not a baseline hazard function (as the name would suggest), and b) employs prediction-at-the-means which is wildly discredited as valid in any practical sense.

And to your earlier point: yes this function makes use of the step function. You can transform that output to a hazard function using smoothing. The worst part of it all, what's the uncertainty interval for that prediction? You may get a Fields medal if you can derive it. I don't think we even know whether bootstrapping works or not.

As an example:

set.seed(1234)
x <- rweibull(1000, 2, 3)
coxfit <- coxph(Surv(x) ~ 1)
bhest <- basehaz(coxfit)
haz <- exp(diff(bhest[, 1])*diff(bhest[, 2]))
time <- (bhest[-1,2] + bhest[-1000, 2])/2
b <- 2^-3

curve(3*b*x, from=0, to=max(x), xlab='Survival time', ylab='Weibull hazard')
points(t <- bhest[-1,2], h <- diff(bhest[, 1])/diff(bhest[, 2]), col='grey')
smooth <- loess.smooth(t, h)
lines(smooth$x, smooth$y, col='red')
legend('topright', lty=c(1,1,0), col=c('black', 'red', 'grey'), pch=c(NA,NA,1), c('Actual hazard fun', 'Smoothed hazard fun', 'Stepped discrete-time hazards'), bg='white')