Difference between brms `unupdate()` and setting priors manually for calculating Bayes Factors

Interfaces brms
, , ,
1

Hello,

I’m calculating Bayes Factors for a brms model using bayesfactor_parameters() and I get wildly different results when manually setting a normal distribution, versus when I allow the function to calculate it’s own priors from the model. By what I can tell from the documentation, it is doing this by using the unupdate() function on the model posteriors.

I manually set priors to my model before running (see below). My question is, should I allow the bayesfactor_parametrs() function to do it’s thing on it’s own, or should I set priors again manually, either with a normal distribution or with my previous priors from the model?

This is for a basic, one predictor mixed effects model.
Formula: DV ~ IV + (1 + IV || participant)

Here are my model and priors

# set priors
fcP <- c(prior(normal(0,1), class='Intercept'), prior(normal(0,1), class='b'),
         prior(gamma(1,1), class='sd'), prior(gamma(1,1), class='phi'))


# make model
fcM <- brm(fcF, family=Beta('logit'), data=beta_p_play_means, fcP,
           chains=8, cores=8, init=0, control=list(adapt_delta=.99),
           threads=threading(2), backend='cmdstanr', normalize=F)

When I use bayesfactor_parameters(fcM) I get a BF of 0.15

But when I extract the psoteriors myself and use a normal distribution prior (note: “slope” is the name of the IV):

fcposterior <- posterior_samples(fcM)
as.numeric(bayesfactor_parameters(fcposterior$b_slope, prior=distribution_normal(8000,0,1)))

I get a BF of 3.47

I know BFs are sensitive to priors but this difference is enormous.
I tried passing my fcP prior df to bayesfactor_parameters() but it won’t accept it as it tries to use unupdate() whenever it detects a brmsfit object.

I am leaning towards feeling that there is no “significant” effect here and so am more inclined to trust bayesfactor_parameters(fcM) without manually specifying priors. But am still a bit confused by the massive difference. I have a similar model with a different predictor (but same priors) and the Bayes Factors from the two methods above are much more (reassuringly) similar.

Hope the info I provided is detailed enough and appreciate any thoughts! 🙂

2

Calling @mattansb

3

Hi @Alex-A14 ,

A few things:

  1. I would use more samples when computing any Bayes factors.

  2. But even with 40,000 samples, results are still highly variable, regardless of method:

library(brms)
library(bayestestR)

mod <- brm(
  mpg ~ hp,
  data = mtcars,
  prior = prior(normal(-0.05, 0.01), class = 'b'),
  
  chains = 10,
  iter = 5000,
  warmup = 1000
)

bayesfactor_parameters(mod)
#> Bayes Factor (Savage-Dickey density ratio) 
#> 
#> Parameter   |       BF
#> ----------------------
#> (Intercept) | 4.37e+20
#> hp          |    64.76
#> 
#> * Evidence Against The Null: 0

samps <- as.data.frame(mod)
bayesfactor_parameters(
  posterior = samps$b_hp, 
  prior = distribution_normal(40000, -0.05, 0.01)
)
#> Bayes Factor (Savage-Dickey density ratio)
#> 
#> BF   
#> -----
#> 54.13
#> 
#> * Evidence Against The Null: 0

bayesfactor_parameters(
  posterior = samps$b_hp, 
  prior = distribution_normal(40000, -0.05, 0.01, random = TRUE)
)
#> Bayes Factor (Savage-Dickey density ratio)
#> 
#> BF    
#> ------
#> 164.25
#> 
#> * Evidence Against The Null: 0

These differences seem big, but look what we get at 3 more runs of the first method (each of these re-calls unupdate() to get new MCMC samples of the prior):

bayesfactor_parameters(mod)
#> Bayes Factor (Savage-Dickey density ratio) 
#> 
#> Parameter   |       BF
#> ----------------------
#> (Intercept) | 4.48e+20
#> hp          |   321.14
#> 
#> * Evidence Against The Null: 0

bayesfactor_parameters(mod)
#> Bayes Factor (Savage-Dickey density ratio) 
#> 
#> Parameter   |       BF
#> ----------------------
#> (Intercept) | 4.44e+20
#> hp          |   304.70
#> 
#> * Evidence Against The Null: 0

bayesfactor_parameters(mod)
#> Bayes Factor (Savage-Dickey density ratio) 
#> 
#> Parameter   |       BF
#> ----------------------
#> (Intercept) | 4.40e+20
#> hp          |    91.23
#> 
#> * Evidence Against The Null: 0

Note that when using unupdate() we are using MCMC and then logspline() to estimate the prior density at the null.

We can also get these exact densities for the prior in this case (but not always, which is by

library(logspline)

d_post <- logspline(samps$b_hp) |> dlogspline(q = 0)
d_prior <- dnorm(0, -0.05, 0.01)

d_prior/d_post
#> [1] 30.05626

Note that this result still suffers from instability in estimating the posterior density of the null value.

When running the model again to get new posterior samples, we again get a different result:

mod2 <- update(mod) # re run MCMC
samps2 <- as.data.frame(mod2)

d_post2 <- logspline(samps2$b_hp) |> dlogspline(q = 0)

d_post2
#> [1] 179.2607

To summarize, density ratios based on MCMC samples are unstable (high Monte Carlo error - MCE), and this instability combines when it is for two densities. In this example we get ratios within one order of magnitude (again, with 40,000 samples, not the default 4,000!).

You can probably decrease the MCE by further increasing the number of samples, but at this point I would probably do a proper BF by comparing the marginal likelihoods of two nested models (with {bridgesampling} directly or bayesfactor_models() convenient wrapper).

Alternatively, Bayes factors for order restrictions are probably more stable (no density estimation involved):

bayesfactor_restricted(mod, hypothesis = c("b_hp < -0.05", "b_hp > -0.05"))
#> Bayes Factor (Order-Restriction)
#> 
#> Hypothesis   P(Prior) P(Posterior)    BF
#> b_hp < -0.05     0.50         0.88  1.76
#> b_hp > -0.05     0.50         0.12 0.235
#> 
#> * Bayes factors for the restricted model vs. the un-restricted model.

bayesfactor_restricted(mod, hypothesis = c("b_hp < -0.05", "b_hp > -0.05"))
#> Bayes Factor (Order-Restriction)
#> 
#> Hypothesis   P(Prior) P(Posterior)    BF
#> b_hp < -0.05     0.50         0.88  1.78
#> b_hp > -0.05     0.50         0.12 0.232
#> 
#> * Bayes factors for the restricted model vs. the un-restricted model.

We can divide both BFs to a get a BF between the two restricted models. Se here we would have a BF of 1.75/0.23 = 7 (data are 7 times more compatible with the slope being smaller than -0.05 compared to it being larger than -0.05).

Or ditch BFs are use posterior bayes inference (pd, ROPE, etc…)