Personally I think the smooth approximations to the baseline hazard are preferable to the piecewise constant approach. In the smooth spline-based approximations there is likely to be less sensitivity to the choice of the number of df or knots. Whereas the number of discrete segments for the piecewise constant baseline hazard likely needs to be chosen by the user, is more context-specific, and the choice is likely to have a bigger influence on the results. At least that is my feeling.
@ermeel @avehtari @jonah Fyi, I abandoned the repo I was working on in my own GitHub page, and instead started a survival branch on the rstanarm repo. Perhaps worth adding this to the wiki @avehtari? It’s still work in progress, but has options for M-splines on the baseline hazard (closed form for cumulative hazard, i.e. fast) or B-splines on the log baseline hazard (uses quadrature, i.e. slow, with 15 quadrature nodes its about 20 times slower than using the M-splines). I’m going to switch the B-splines out for R&P restricted cubic splines instead though, since they are likely to be more stable I think (as shown in the simulation study in this paper).
I think once you introduce spline-based time-dependent effects, the need for quadrature becomes near unavoidable, hence why I included it as an option in the stan_surv function on that new rstanarm branch. But, one suggestion from a colleague of mine was to use the approach described in Section 2.2 of this paper where they use the R&P approach and just replace the hazard with a small value (e.g. 1E-8) when the hazard strays negative. The approach feels a little funky to me, since you are basically fudging the likelihood, but they seem to claim it works ok, so I might try it out. But I would be interested to know whether you guys think it would work in a Bayesian setting.