How to do bridge sampling and calculate Bayes Factor with brms with cmdstanr backend?

I’m using `brm()` with `cmdstanr` as its backend (i.e. I set `backend = "cmdstanr"` in formulae) and trying to calculate the Bays Factor using `bayes_factor()` comparing multiple models.

However, `bayes_factor()` returns `Error: Backend 'rstan' is required for this method` and no computation shows progress. Actually, I really need to use `cmdstanr` to analyse more complicated data repetitively with speed, and I dis-prefer to using `rstan` this time. Is there any way to compute BF with brms with cmdstanr backend?

MWE

``````## For sleepstudy data
library(lme4)

## Null model
fm0 <- brm(
Reaction ~ 0 + Intercept + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
prior(normal(0, 1), class = b, coef = Intercept),
prior(normal(0, 1), class = sd),
prior(lkj(2),       class = cor)
),
save_pars = save_pars(all = TRUE),
backend = "cmdstanr"
)

## Alternative model
fm1 <- brm(
Reaction ~ 0 + Intercept + Days + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
prior(normal(0, 1), class = b, coef = Intercept),
prior(normal(0, 1), class = b, coef = Days),
prior(normal(0, 1), class = sd),
prior(lkj(2),       class = cor)
),
save_pars = save_pars(all = TRUE),
backend = "cmdstanr"
)

## The following commands return an error suggesting us to use rstan
bayes_factor(bridge_sampler(fm1), bridge_sampler(fm0))
bayes_factor(fm1, fm0)
``````

I’m running on:

• Operating System: Windows 10 x64 (build 18362)
• CmdStan Version: 0.1.3
• brms Version: 2.14.0

Any suggestion is appreciated.

That currently only works for the subset of models that can also be compiled via `rstan`, see the discussion here: https://github.com/quentingronau/bridgesampling/issues/27

Thank you for your quick reply and letting me know the thread to discuss the issues. I will try the transformation of the model calculated with cmdstanr, recompilation with rstan, and the bridge simpling.

I noticed that I don’t know how to save (or export) the sampling results from `brm` object. In the MWE in my previous post, I can save the fitting results of `brm` to `.rds` by adding the `file` option, i.e.:

``````fm0 <- brm(
Reaction ~ 0 + Intercept + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
prior(normal(0, 1), class = b, coef = Intercept),
prior(normal(0, 1), class = sd),
prior(lkj(2),       class = cor)
),
save_pars = save_pars(all = TRUE),
backend = "cmdstanr",
warmup = 2000,
iter = 10000,
chains = 8,
save_model = "fm0",
file = "fm0" #rds file
)
``````

Obviously, this doesn’t save the sampling result as csv and, therefore, `rstan::read_stan_csv(fm0\$output_files())` is not applicable.

So, is there any way (1) to save the sampling results of `brm` to a `.csv` file, or (2) to make `rstan` read a `.rds` file? I would like to ask my question to @paul.buerkner , too, and I tagged him here.

I am not sure for what you want to use this approach? brms stores the rstan file in the \$fit slot, which you can try to reuse if needed.

Thank you for your reply. I will try \$fit later. Actually, I would like to convert cmdstanr object made via brm() to rstan object, so that I can compile models faster with the cmdstanr and calculate Bayes Factor with rstan, as suggested by @Henrik_Singmann ( https://github.com/quentingronau/bridgesampling/issues/27 )

Would `brm()` with cmdstanr backend save the sampling results in `\$fit`? Are such sampling results readable/reusable to `rstan`? I would need the sampling results primarily to calculate Bayes Factor, rather than the stan code/stan file.

Yes it is stored in \$fit. They are not usable for bayes factors at the moment which will be a feature implemented in cmdstan® in the future. There is nothing i can do from the brms side for now.

OK, I appreciate your letting me know that! This time, I will use brm() with rstan background, looking forward to the implemention of the Bayes Factor calculation in cmdstanr…

I think I finally managed to calculate the Bayes Factor with `brm()` and `cmdstanr` backend bypassing `rstan`. In the following MWE, (1) I calculated the parameters of the model (`fm0` and `fm1`) by `brm()` with `cmdstanr` backend, (2) compiled the same model with `rstan` backend but without any sampling (`rrfm0` and `rrfm1`), and finally (3) compute the Bayes Factor with `bridgesampling::bridge_sampler()`, extracting the samples by `fmX\$fit` and the parameters by `rrfmX\$fit`.

The result of the Bayes Factor test seems to be reasonable for me, since the model with the effect of days (`fm1`) is preferred to the null model, as expected, by BF = 67.8.

Do you think this calculation is valid? Would you find any erroneous computation in the process?

``````library(lme4) # to get sleep study data
fm0 <- brm(
Reaction ~ 0 + Intercept + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
prior(normal(0, 1), class = b, coef = Intercept),
prior(normal(0, 1), class = sd),
prior(lkj(2),       class = cor)
),
save_pars = save_pars(all = TRUE),
backend = "cmdstanr",
warmup = 2000,
iter = 10000,
chains = 8
#save_model = "fm0",
#file = "fm0" #rds file
)

rrfm0 <- brm(
Reaction ~ 0 + Intercept + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
# global parameters
prior(normal(0, 1), class = b, coef = Intercept),
prior(normal(0, 1), class = sd),
prior(lkj(2),       class = cor)
),
save_pars = save_pars(all = TRUE),
backend = "rstan",
iter = 0,
chains = 0
)
bridgesampling::bridge_sampler(
samples = fm0\$fit,
stanfit_model = rrfm0\$fit
)

fm1 <- brm(
Reaction ~ 0 + Intercept + Days + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
# global parameters
prior(normal(0, 1), class = b, coef = Intercept),
prior(normal(0, 1), class = b, coef = Days),
prior(normal(0, 1), class = sd),
prior(lkj(2),       class = cor)
),
save_pars = save_pars(all = TRUE),
backend = "cmdstanr"
)

rrfm1 <- brm(
Reaction ~ 0 + Intercept + Days + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
# global parameters
prior(normal(0, 1), class = b, coef = Intercept),
prior(normal(0, 1), class = b, coef = Days),
prior(normal(0, 1), class = sd),
prior(lkj(2),       class = cor)
),
save_pars = save_pars(all = TRUE),
backend = "rstan",
iter = 0,
chains = 0
)

bridgesampling::bridge_sampler(
samples = fm1\$fit,
stanfit_model = rrfm1\$fit
)

bayes_factor(
bridgesampling::bridge_sampler(
samples = fm1\$fit,
stanfit_model = rrfm1\$fit
),
bridgesampling::bridge_sampler(
samples = fm0\$fit,
stanfit_model = rrfm0\$fit
)
)
``````