# Survival analysis - plot empirical hazard function first?

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

• You can just use a reverse trapezoid rule. When you run basehaz you return an expression of time versus the cumulative hazard. Apply consecutive differences between cumulative hazards and divide by their respective time intervals, apply to the later time point Feb 2 at 21:03

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

• Thank you, but I was more interested in the hazard function over time. Nov 25, 2022 at 7:01
• what about plot(fit, cumhaz = TRUE)? Nov 28, 2022 at 21:29

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.

• ok thank you - I can take all of that on board. But surely there is still a way to do this in R? It's a simple command in Stata... Nov 27, 2022 at 4:33

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