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

Interfaces CmdStan
1

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.

2

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

3

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.

4

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.

5

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.

6

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 )

7

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.

8

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.

9

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…

10

@Henrik_Singmann and @paul.buerkner

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
  )
)