Using stan_jm for parametric proportional hazards regression only


Might be, but I wasn’t aware of it. Thanks anyways.


I think I oversaw something here: This works if the Poisson regression is used with identity link and doesn’t work with logarithmic link, unless you find a way to write

\log{\left[e^{x_1' \beta} + \cdots + e^{x_n'\beta}\right]} = f(x_1,\cdots,x_n) \beta

for any choice of \beta, n and x_i, and f is a (row-)vector valued function.


@ermeel Yeah, specifying assoc = NULL in stan_jm means that the likelihood for the longitudinal and survival submodels are effectively independent, so in theory you should get estimates for the survival submodel that are equivalent to fitting the survival model alone. But, you must provide at least one longitudinal measurement per individual, so you would have to supply dummy data, and the joint posterior of all parameters is being used so there is a lot of unnecessary work happening (maybe even influencing the quality of the estimates?). In short, it would be an undesirable way to fit a standard survival model!

Re Poisson time splitting, yeah, you can get an equivalent model to the Cox model by splitting at all event times as you suggested. Here is another reference on this approach. It’s been around a long time, but the splitting means you have lots of rows of data, and so the computations weren’t really feasible until more recently I guess. But I find it is a pretty clunky approach, having to split/expand your data like that.

Royston-Parmar flexible parametric models are much better, I think. They model the log cumulative baseline hazard as a smooth function using cubic splines, and they don’t require the splitting of your data. I’ve wanted to implement them (and standard parametric survival models) either in rstanarm or a new package, but haven’t got started yet. I think I am tending towards a separate package from rstanarm, because there is scope to include competing risks, multi-state models, etc, and incorporating all of it in rstanarm starts to seem infeasible.


Either way is fine by me. If it doesn’t jive nicely with your code already in place for jm and mvmer then a separate package could certainly be more feasible. But if you end up wanting to add them to rstanarm that’s fine too.


Hmm… you are right. I recall that this approach is somehow “between the lines” in the text books, but never explicitly discussed (at least in my text books on the matter).

Interesting approach. However, while it is common practice to model hazards, these are conditional on post-baseline risk sets. Thus, their causal interpretation can be challenged.

A package like this is really missing in the Stan universe! If you have the means for it, then go for it… it’s probably not a huge user base, but those would be very thankful for this.


I’d be interested in this and would be happy to collaborate/contribute.

What precisely do you mean with post-baseline risk sets? Sorry in advance if this is again an effectively well-known nomenclature/concept which I’m not aware of.


The hazard is defined based on conditioning on the risk set at some time point t past baseline. We all know that messing with post-baseline stuff is not a good idea usually. It can be OK, but it does rely on assumptions - and I would argue that most of us keep forgetting this fact in TTE analyses. It is a bit like using post-baseline time-varying covariates in your analysis which makes your analyses rest on assumptions (assumptions are fine as long as you know them and know their limitations, etc.).


OK. So you are actually referring to something like the partial likelihood in (9.38) chapter 9.4 page 145 in Computer Age Statistical Inference, which is based on the Lemma above Eq. 9.37? See also The Efficiency of Cox’s Likelihood Function for Censored Data.


You are right. Sorry for any confusion (the parametric approaches above should not be affected by this…but they are heavily criticized for being parametric, usually).


Yes! I’ve been lamenting the lack of user-friendly survival models in Stan. There are a lot of potential models to add (parametric AFT and PH, flexible baseline hazard PH models, competing risks, semicompeting risks/multistate models, etc) and enough machinery shared among them that shoving it all into rstanarm seems infeasible. I wish there were a better way to borrow awesome recentering and other tricks from rstanarm without adding a new model there.

@ermeel, if you’re willing to consider non-Stan options and want an answer today, the BayesSurv_HReg function of the SemiCompRisks package implements a Bayesian Cox-like proportional hazards model with a piecewise exponential baseline hazard. (Gory details on the flexibility: it’s a time-homogenous poisson process, meaning it assumes the number of intervals for chunking the baseline hazard time scale is poisson-distributed. Then the values of the log-baseline hazard on those intervals are MVN-ICAR, which smoothes everything a bit by making intervals closer together in time have more similar baseline hazard values than intervals farther apart in time. I do not think Stan could fit this model because of the discrete parameter for the number of intervals, but splines for the cumulative hazard could give similar answers.)

@wds15, are you referring to the “hazard of hazard ratios”? I see no inherent problem with using hazard models, even for causal contrasts. The problem is interpreting hazard ratios causally because they have selection bias baked in through that troublesome conditioning-on-survival in the definition. Even so, you can still use a hazard model to estimate other (causal) quantities.


Funny way of saying this! So yes, this selection bias thing is what I am concerned about. That does not defeat a causal interpretation if you get to estimate them right.

A Stan survival package would be awesome, but a major effort!


This has been on my to-do list for a while but, realistically, I probably won’t have time to do it anytime soon. brms can actually already fit some survival models, but I was hoping to have precompiled versions at some point. The functions in the survival package can be used to parse formulas and get the data prepared, and these models can be precompiled, so that all fits well with the rstanarm paradigm. Someone just needs to add it to rstanarm or copy the rstanarm strategy in a separate package. If anyone is interested in working on that we can definitely provide guidance along the way.

Until then there are a few options that don’t do exactly this but are related:


Actually I would be interested to work on it (as a side project). It would be a nice way for me to give something back to the Stan community and learn more about survival models. And I’d really appreciate any guidance!


Just to add that survHE is also available as a development version from The version on CRAN is OK for the purposes of solely estimating the survival parameters. There are other functions (that are specific to health economic evaluation) that have been updated slightly in the development version, though.

In fact, we’re working to expanding survHE so comments and/or suggestions more than welcome!


I would definitely be interested in working on this. I have been talking about working on it for about 18 months now but have never gotten around to it. Happy to help out in any way.


@kaybenleroll and @ermeel or anyone else interested, can you start a separate topic for this? We can discuss the best way to collaborate on this and whether it should be part of an existing package or a new one.

@Gianluca_Baio, Thanks for the update. I’ll definitely take a closer look at survHE! I looked enough that I felt comfortable recommending it here but I’m not super familiar with it yet.

Survival models in rstanarm
Survival models in rstanarm

@jonah, could you point me to an example how to do this? I would like to make my recent i-spline-based Royston & Parmar Stan code somewhat easier to use through this.


I haven’t really tried it with the survival package, but from what I recall it’s similar to lme4 and other packages we’re emulating in rstanarm in that you can retrieve model matrices and all sorts of metadata from the objects it returns. The survreg function (and others from the survival package) probably does a lot of the same preprocessing that is needed.


Thanks @jonah . Indeed looks very helpful!


Great. Our new package (or rstanarm) can certainly import the survival package in order to use it internally. I haven’t looked at its license yet, but we can also probably copy and tweak its code and use it internally, acknowledging the appropriate copyrights (see e.g. the authors section in the rstanarm DESCRIPTION file). So you’re not restricted to just calling its functions, you can probably copy chunks of code too.