PPC for Cox model

[fweber144]

A posterior predictive check (PPC) for the Cox model would be nice, but it would probably have to be conducted differently than in all other models since the spline estimate of the baseline hazard makes it hard to sample from the response distribution. Section 4.4 ("Standardised survival probabilities") of the preprint by Brilleman et al. (2020) describes a PPC which is feasible in this situation. In my understanding (and only under right-censoring, I guess), this PPC basically compares the Kaplan-Meier estimate of the survival curve (i.e. the estimated CCDF) for the observed data to the posterior CCDFs which are averages of the individual posterior CCDFs. Each individual posterior CCDF is obtained by plugging in the predictor values for one individual from the observed dataset. With the default settings, there would be 4000 posterior CCDFs, so as usual, only a random subset of them will be used for the overlay plot. The resulting plot should be similar to bayesplot::ppc_ecdf_overlay(), but taking the right censoring of the observed data into account and applying the "CCDF = 1 - CDF" transformation (although this transformation is not strictly necessary and would just be a convention from traditional survival analysis).

In principle, the posterior CCDFs could be computed in the generated quantities block, but for better generalizability (e.g. later also CCDF prediction for new data), it would probably be desirable to compute them in R.

The same style of PPC could also be used for the non-Cox time-to-event models already existing in brms (e.g. right-censored log-normal or right-censored Weibull), but with an easier implementation because sampling from the (uncensored) response distribution is easy for these distributions. (Sampling leads to an easier implementation because it doesn't need extra code to compute the posterior CCDFs exactly, but also because it doesn't need individual posterior CCDFs and their averaging.) However, this is a different topic, so I'll open a new issue for that.

[paul-buerkner]

I agree. As you point out, the key point here is to provide sampling from the time-to-event distribution and I don't know yet how to do this for the Cox model, but I haven't invested time in this topic yet. Once figured out, we simply write a posterior_predict_cox method and all the currently existing PPC stuff will work.

Edit: I will not do any special case coding for PPCs in the Cox model. This needs to go in more specialized packages, not in brms. For brms, we need to find a way to sample from the response distribution.

[fweber144]

Edit: I will not do any special case coding for PPCs in the Cox model. This needs to go in more specialized packages, not in brms. For brms, we need to find a way to sample from the response distribution.

OK, I totally understand that. But then I think this will be hard to realize. For example, inverse transform sampling would require the inverse CDF which in turn requires the inversion of the I-spline. This would probably have to be done numerically, e.g. using Stan's algebra_solver. I don't know if the computational burden of this would still be practical. Besides, inverse transform sampling also requires the evaluation of the I-spline at an arbitrary timepoint which will probably be hard to realize, too (at least here, where the results from splines2 have to be passed to the Stan model). I found this related thread in The Stan Forums, but it seems like they also didn't find an easy solution.

With these thoughts, I don't want to talk you into doing special case coding for Cox model PPCs. I just wanted to share some thoughts about the limitations of inverse transform sampling. But perhaps there are other sampling approaches out there which might be applicable here?

[fweber144]

Just to make this clear: A PPC for the Cox model would only be a "nice-to-have" feature, not a must-have. So if you want to close this issue, I'm ok with it. I also think that brms doesn't need to handle all kinds of special models in a comprehensive manner. The yet unpublished survival feature of rstanarm seems to be the better place for special time-to-event models. Btw, the Cox model is not my favorite time-to-event model anyway (I would prefer a parametric one, such as the log-normal or the Weibull model).

[paul-buerkner]

I realize my edit came out harsher then intended :-D Sorry for that. On my list of ideas for research projects, I actually have the point of figuring out how to sample from the PPD of cox models with spline baseline harzards (for posterior_predict) and for the mean of the PPD (for posterior_epred), and this issue is a good reminder to take some time to investigate this at some point (or find someone to investigate this).

[fweber144]

Yes, that's definitely an interesting research topic.