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

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

