Bayesian implementation of the Royston-Parmar survival model

Does anyone know if there is any Bayesian implementation of the Royston-Parmar survival model?

Have a look here. Failing that, @ermeel and @sambrilleman (and possibly @jonah) are the right people to ask.

1 Like

Check out our arXiv paper here:
It describes the two approaches we took: 1) use M-splines on the hazard, or 2) use B-splines on the log-hazard. The M-splines approach is faster to estimate because it allows you to calculate the cumulative hazard in closed form and therefore no quadrature is needed.

We looked into implementing the Royston Parmar model but, from memory, the difficulty is that the restricted cubic splines do not constrain the (log) cumulative hazard to be monotonically increasing – instead this requirement is maintained because the data helps to ensures it, rather than because it is explicitly a constraint in the model formulation. I think Royston or Parmar mention that somewhere in their paper. This is fine in a maximum likelihood setting, but in a Bayesian world we need to explore the entire parameter space for the model – and unfortunately when the cumulative hazard is allowed to decrease, then the hazard is negative and the log-likelihood isn’t defined (i.e. log of a negative number).

So MCMC doesn’t really work, unless you can figure out a way to get around that issue. Hence why we went for M-splines on the hazard or B-splines on the log hazard – they aren’t the Royston Parmar model but they are a flexible parametric equivalent.

Hope that helps…


Here is some Stan code I wrote some time ago, that is very close to Royston& Parmar (it uses M-splines instead of restricted cubic splines and a QR transformation for any covariates). I can recommend the arXiv paper @sambrilleman mentioned, for a practical implementation and the thread @maxbiostat mentioned contains a lot of the reasoning that influenced the stan_surv implementation.


I’ve got a working Bayesian implementation of the Royston-Parmar survival model done in JAGS
based on code in Freeman & Carpenter (

There’s a little trick to deal with the “restricted cubic splines do not constrain the (log) cumulative hazard to be monotonically increasing” issue. I’m not sure if this little trick is strictly allowed, but it seems to work.

Please take a look and let me know your thoughts:

thoughts? @maxbiostat, @ermeel @sambrilleman, @jonah?


Hello, I am a clinician. It would be nice for me to be able to select the degrees of freedom of the baseline hazard and covariate effects to allow for ver flexible curves (i.e., crossing survival curves). It would also be nice to have posterior probabilities of therapeutic effect, as well as survival differences, with their respective HDIs, to test posterior probabilities as a function of time.

Bit late here but check out the survHE package:
survHE: Survival Analysis for Health Economic Evaluation and Cost-Effectiveness Modeling | Journal of Statistical Software (
It has Royston-Parmar (just the Weibull models though).

1 Like