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

Question

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)),
                                        ContrastType = factor(rep(c("Head", "Tail"),
                                                                  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

Answer 1

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.

See also Bayes factor (B) vs p-values: sensitive (H0/H1) vs insensitive data.

---

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.