Sampling Bayes factors under the null hypothesis to estimate a threshold of "significance" for hypothesis testing

Context: I have a psychology experiment with a 2 x 2 design (with Condition (label, no label) and ContrastType (head, tail) as my two factors) where I want to estimate whether the mean of each subgroup is different from zero, using Bayes factors. In the sample theory based framework I would do something like this in R:

library(lmer); library(emmeans)

model <- lmer(ChanceArcsin ~ ContrastType*Condition + (1 | Participant))
posthoc <- emmeans(model, ~ ContrastType | Condition,
options = list(infer = c(T, T),
null = 0,
level = .95))

Now, if I want to have Bayes factors instead of $$p$$-values, I would use brms::hypothesis like that:

library(brms)

p <- c(set_prior("uniform(-.8,.8)",
class = "Intercept"),
set_prior("normal(0,.5)", class = "b"))
model <- brm(ChanceArcsin ~ ContrastType*Condition + (1 | Participant),
prior = p, family = gaussian(),
chains = 4, cores = 4, iter = 2000,
save_all_pars = T)
hyp.test <- hypothesis(model,
c("Intercept > 0",
"Intercept + ContrastTypeTail > 0",
"Intercept + ConditionNo Label > 0",
paste("Intercept +",
"ConditionNo Label +",
"ContrastTypeTail +",
"ContrastTypeTail:ConditionNo Label",
"> 0")))

The problem here is that, with my data, I have Highest Posterior Density Intervals not far from zero and non-significant results from lmer + emmeans, but fairly high Bayes factors for a psychology experiment (between 10 and 46). I was a bit surprised and I ran some simulations on my design to get a better idea of the likelihood of observing such $$b$$-values under the null hypothesis (full code at the bottom). What I found was that the 95th percentile for Bayes factors was 14.5, and 25.5% of the Bayes factors from my simulations where above 3 (the somewhat common threshold value for results worth mentioning).

Questions: Is it okay/a good idea to use the 95th percentile (or whatever other value) as a threshold to interpret my null hypothesis test Bayes factors as being worth mentioning or not?

If I do so, is it better to run simulations with a larger sample size on my design, or use the actual sample size I have (that is, 69 observations, unequally divided between the different cells)?

I also saw this question but the link to the course is broken and I'm not very familiar with the idea of minimal Bayes factor. The reason why I want to report Bayes factors is that Bayesian analysis is not yet too widespread in psychology, and I'm not sure how publishable would be an article without any kind of explicit hypothesis testing, sadly...

Simulation code (a bit long):

# LIBRARY IMPORTS =======================================================================
library(lme4)
library(emmeans)
library(brms)
library(tidyverse)

# BAYESIAN POINT NULL HYPOTHESIS TESTING ================================================
H_naught.test <- function(N){
# Generate random data
# mean to 0, sd=.5 from real data
new_old.sims <- replicate(100,
list(tibble(ChanceArcsin = rnorm(N*2, 0, .5),
Participant = rep(1:N, each = 2),
Condition = factor(rep(c("Label", "No Label"),
times = N/2, each = 2)),
times = N)))))

# Define STB analysis function, returning emmeans analysis
stb.analysis <- function(df.list){
e <- lapply(df.list,
function(df){
m <- lmer(ChanceArcsin ~ ContrastType*Condition +
(1 | Participant),
data = df)
t <- emmeans(m, ~ ContrastType | Condition,
options = list(infer = c(T, T), null = 0,
level = .95)) %>%
as_tibble()
return(t)
})
return(bind_rows(e))
}

# Define Bayesian analysis function, returning hyp. test bf
bayesian.analysis <- function(df.list){
p <- c(set_prior("uniform(-.8,.8)",
class = "Intercept"),
set_prior("normal(0,.5)", class = "b"))
m.list <- brm_multiple(ChanceArcsin ~ ContrastType + Condition +
ContrastType:Condition +
(1 | Participant),
data = df.list, prior = p, family = gaussian(),
chains = 4, cores = 4, iter = 2000,
save_all_pars = T, combine = F)
bf.list <- lapply(m.list,
function(m){
bf <- hypothesis(m,
c("Intercept > 0",
"Intercept + ContrastTypeTail > 0",
"Intercept + ConditionNo Label > 0",
paste("Intercept +",
"ConditionNo Label +",
"ContrastTypeTail +",
"ContrastTypeTail:ConditionNo Label",
"> 0")))
return(as_tibble(bf\$hypothesis))
})
bf <- bind_rows(bf.list)
return(bf)
}

# Get evidence summary from simulations
t <- proc.time()
new_old.sims.stb.evid <- stb.analysis(new_old.sims)
stb.time <- proc.time() - t
t <- proc.time()
new_old.sims.bayesian.evid <- bayesian.analysis(new_old.sims) %>%
mutate(Condition = factor(ifelse(grepl("NoLabel", Hypothesis), "No Label", "Label")),
ContrastType = factor(ifelse(grepl("Tail", Hypothesis), "Tail", "Head")))
bayesian.time <- proc.time() - t

return(list(BayesianEvidence = new_old.sims.bayesian.evid,
SampleTheoryEvidence = new_old.sims.stb.evid,
}

new_old.sims.results.200 <- H_naught.test(N = 200)
bayesian.95 <- quantile(new_old.sims.results.200$$BayesianEvidence$$Evid.Ratio, .95)
stb.95 <- quantile(new_old.sims.results.200$$SampleTheoryEvidence$$p.value, .05)
bayesian.above3 <- sum(new_old.sims.results.200$$BayesianEvidence$$Evid.Ratio > 3) / 400
stb.bellow_dot05 <- sum(new_old.sims.results.200$$SampleTheoryEvidence$$p.value<.05) / 400
• There is an error in the code for the Bayesian hypothesis testing. It's necessary to make two changes: (1) Set sample_prior = "yes" in the brm_multiple call; and (2) Specify the intercept prior as set_prior("uniform(-0.8, 0.8)", class = "b", coef = "Intercept"). Take a look at the brm_multiple documentation for details. Commented Oct 13, 2023 at 15:04
• You'd still get a large proportion of Bayes Factors > 3 as this is not particularly strong evidence. See Held, Leonhard; Ott, Manuela (2018). On p-Values and Bayes Factors. Annual Review of Statistics and Its Application, 5(1):593-419. Perhaps "Bayes Factor worth mentioning" is the Bayesian equivalent of "p-value trending towards significance". Neither is saying very much. Commented Oct 13, 2023 at 15:07
• I would suggest to use posterior probabilities instead or even better -- posterior credible intervals. That is, switch from "an evidence framework" to "an estimation framework". I think the results would be easier to interpret and explain that way as well. Commented Oct 13, 2023 at 15:09
• I'm fairly sure that the prior syntax in this was correct when I posted (5 years ago), though it's very likely it has changed since then. Not sure about the sample_prior option. And agreed on using credible intervals instead but as mentioned the context in my field did not allow for that 5 years ago (unsure about now, probably still true). I'm no longer working in this field but thanks for the article, looks interesting! You should post this as an answer in case other people stumble upon this question with the same issue :) Commented Oct 18, 2023 at 14:41

There is no universally accepted "strength of evidence" scale for Bayes factors.

Here is a side-by-side comparison of three proposed scales:

Image from Held, Leonhard; Ott, Manuela (2018). On p-Values and Bayes Factors. Annual Review of Statistics and Its Application, 5(1):593-419. https://doi.org/10.1146/annurev-statistics-031017-100307

The two more recent scales, by Goodman and Held & Ott, agree the evidence is (at least) strong if the BF is ≥ 30. This would suggest to interpret the BFs reported by the OP (between 10 and 46) more cautiously.

That being said, I think that Bayes factors are something of a cop-out if one is interested in doing Bayesian analysis. BFs try to be a compromise between frequentist and Bayesian methods; in a way, to make the Bayesian approach more palatable to those who are used to classical statistics. (The OP suggests as much himself.)

It may be more straightforward to plot the posterior distribution of the quantities of interest — here the four linear effects — and calculate relevant posterior probabilities; for example: $$\operatorname{Pr}\left\{\beta_{kl} > 0\right\}$$ and/or $$\operatorname{Pr}\left\{- b < \beta_{kl} < b \right\}$$ where the region $$[-b,b]$$ designates a region of "practical equivalence" around the null hypothesis of 0. We can think of this as switching from evidence-based framework to estimation-based framework.

PS. The specification of the priors and the model needs to be updated to test both one-sided (a > b) and two sided (a = b) hypotheses with newer brms versions. (I have 2.20.4.)

prior <- c(
set_prior("normal(0, 0.5)", class = "b"),
set_prior("uniform(-0.8, 0.8)", class = "b", coef = "Intercept")
)
m.list <- brm_multiple(
Y ~ 0 + Intercept + X1 * X2 + (1 | Participant),
family = gaussian(), data = df.list,
prior = prior, sample_prior = "yes",
chains = 4, cores = 4, iter = 2000,
combine = F # save_all_pars = T,
)

There a couple of changes:

• 0 + Intercept turns off centering of the design matrix — the default behavior — so that a prior can be defined on the real intercept directly.
• Since the intercept is now another population-level parameter, its prior is of class = "b", not class = "Intercept".
• sample_prior is set to "yes" to draw samples from the priors for the hypothesis testing.
• save_all_pars is deprecated, so it's commented out.