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

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