Survival analysis - plot empirical hazard function first?

Question

I am wondering how analysts might generally approach a survival analysis problem? I used to naively think just to throw a Cox model at everything but am becoming attracted to the idea of using parametric survival models on a more regular basis because they return the baseline hazard. Of course this means that you should really have some idea of what the actual shape of the hazard function is as that will dictate what parametric distribution you choose to model survival time with. But even if you plan to use a Cox model there must still be utility in visualising the empirical hazard functions grouped by whatever main exposure of interest you have - as this will assist in determining if hazards are proportional even if the model doesn't directly estimate them.

So, if people generally agree with that logic (please explain if you don't) I want to know how to plot empirical hazard functions, as a first step in approaching survival analyses going forward.

In Stata, I think I have worked out that you can simply do this with sts graph, e.g. sts graph, hazard by(x1) if you want two curves according to the two levels of x1.

But how does one do this simply in R? (basehaz doesn't provide the hazard rate but instead the cumulative hazard rate.)

Answer 1

Use survival::survfit to fit a Kaplan-Meier curve. As a simple example of comparing an empirical and fitted curve:

library(survival)
set.seed(4949)
x <- rexp(n = 20, rate = 0.2)

sur <- Surv(time = x)
fit <- survfit(sur ~ 1) # fit Kaplan Meier
reg <- survreg(sur ~ 1, dist = "exp")
plot(fit) # empirical
lines( # fitted
  x = seq(0.01, 14, by = 0.01)
  , y = pexp(q = seq(0.01, 14, by = 0.01)
             , rate = 1/exp(coef(reg))
             , lower.tail = FALSE)
  , col = "red"
)

Answer 2

The cumulative hazard function is just the integral of the instantaneous hazard function over time. Anything that you want to do with the hazard function over time can be done--probably more effectively--with the cumulative hazard function.

Hazard functions over time per se aren't very useful in Cox models, as they can jump around a lot. The estimated cumulative baseline hazard $\hat H_0(t)$ is typically much more useful, as it represents overall trends a bit more smoothly and can still be used to compare against the cumulative baseline hazard of any parametric form that you'd like to consider. You can also use cumulative hazard curves grouped by exposures of interest to evaluate proportional hazards (PH) graphically in simple cases. In more complicated models, evaluating PH is probably better done via plotting scaled Schoenfeld residuals over time after fitting a proportional-hazards model.

Answer 3

I think the closest I've come to a solution (and this seems to work), is the bshazard package.

https://journal.r-project.org/archive/2014/RJ-2014-028/RJ-2014-028.pdf