Dear all,

I’m interested in replicating a frequentist analysis using Bayesian methods (first using essentially non-informative priors, second using informative priors and comparing). I am planning to use the brms package, but open to other options if necessary.

In brief, the original analysis uses a Cox proportional hazards model (without violations of the proportional hazards assumption). The model is adjusted, but uses no time-transformations or other advanced features.

While I am aware of the possibility of using parametric survival models, this is not what I am interested in here. I am ware that the Cox PH model models relative differences (using partial MLE and ignoring the baseline hazards); however, it seems that the frequentist survival R package includes the capability of predicting absolute risks for new patients at specific time-points; I’ve looked at the code, but have not been able to dis-entangle the relevant parts to manually do this using, e.g., brms, which supports the Cox model, but does not include functionality to generate posterior predictions.

I’ve searched the forums and found no previous answer that matches this question, so I hope that somebody on here can help.

Here’s a simple reproducible example showing how it’s done using the survival package (and fitting of the model using brms):

```
library(survival)
# Simple randomly generated dataset - in the real use-case, adjustments would be needed, but no time-varying effects, etc.
set.seed(12345)
dt <- data.frame(trt = 1:200 > 100,
time = sample(1:90, size = 200, replace = TRUE),
event = c(1:100 < 35, 1:100 < 25))
fit <- coxph(Surv(time, event) ~ trt, data = dt)
summary(fit)
# Predictions for new patients at the maximum follow-up time - works
predict(fit, type = "survival", newdata = data.frame(trt = c(F, T), time = c(90, 90), event = c(T, T)))
# Using brms
library(brms)
bfit <- brm(time | cens(1-event) ~ trt, data = dt, family = cox(), cores = 4, backend = "cmdstanr")
summary(bfit)
# Model works, but predictions are not supported, i.e., the following leads to an error:
#posterior_epred(bfit, newdata = data.frame(trt = c(F, T), time = c(90, 90), event = c(T, T)))
```

If anybody has figured out a way to similarly generate predictions from the Bayesian model, it’d be greatly appreciated, as would any suggestions for how this could possibly be approached.

So far, my best solution has been to fit the model using both brms and survival::coxph with identical data/formulas, and then for each set of posterior draws overwrite the coefficients of the frequentist model and use that for predictions; however, in examples similar to that above, all patients in the reference group get similar values as the intercepts from the brms model cannot be used to overwrite any parameters in the frequentist model, so it doesn’t solve the issue of obtaining full posteriors for each group (and if I use an informative prior, e.g. a neutral but relatively narrow one for the treatment effect, it will only pull the non-reference group towards the reference, instead of pulling each two closer to each other, as normally happens).

- Operating System: Windows 10 x64
- brms Version: 2.17.0
- R version: 4.2.1
- cmdstan version: 2.29.2

Thanks in advance and best regards,

Anders