bayesfactor_models: Bayes Factors (BF) for model comparison

In bayestestR: Understand and Describe Bayesian Models and Posterior Distributions

Bayes Factors (BF) for model comparison

Description

This function computes or extracts Bayes factors from fitted models.

The ⁠bf_*⁠ function is an alias of the main function.

For more info, see the Bayes factors vignette.

Usage

bayesfactor_models(..., denominator = 1, verbose = TRUE)

bf_models(..., denominator = 1, verbose = TRUE)

## Default S3 method:
bayesfactor_models(..., denominator = 1, verbose = TRUE)

Arguments

...	Fitted models (see details), all fit on the same data, or a single BFBayesFactor object (see 'Details'). Ignored in as.matrix(), update(). If the following named arguments are present, they are passed to insight::get_loglikelihood() (see details): estimator (defaults to "ML") check_response (defaults to FALSE)
denominator	Either an integer indicating which of the models to use as the denominator, or a model to be used as a denominator. Ignored for BFBayesFactor.
verbose	Toggle off warnings.

Details

If the passed models are supported by insight the DV of all models will be tested for equality (else this is assumed to be true), and the models' terms will be extracted (allowing for follow-up analysis with bayesfactor_inclusion).

• For brmsfit or stanreg models, Bayes factors are computed using the bridgesampling package.

  • brmsfit models must have been fitted with save_pars = save_pars(all = TRUE).

  • stanreg models must have been fitted with a defined diagnostic_file.

• For BFBayesFactor, bayesfactor_models() is a wraparound BayesFactor::extractBF().

• For all other model types, Bayes factors are computed using the BIC approximation. Note that BICs are extracted from using insight::get_loglikelihood, see documentation there for options for dealing with transformed responses and REML estimation.

Additional methods

The resulting output is supported by the following methods:

• as.matrix(): Extract a full matrix of (log-)Bayes factors between all models (using the transitivity of Bayes factors).

• update(): subset and/or re-reference the Bayes factors to a different model.

• as.numeric(): Extract the (possibly log-)Bayes factor values.

See examples and bayesfactor_methods.

Value

A data frame containing the models' formulas (reconstructed fixed and random effects) and their log(BF)s (Use as.numeric() to extract the non-log Bayes factors; see examples), that prints nicely.

For as.matrix() a square matrix of (log) Bayes factors, with rows as denominators and columns as numerators.

Prior and posterior considerations

In order to correctly and precisely estimate Bayes factors, a rule of thumb are the 4 P's: Proper Priors and Plentiful Posteriors.

For the computation of Bayes factors, the model priors must be proper priors (at the very least they should be not flat, and it is preferable that they be informative) (Note that by default, brms::brm() uses flat priors for fixed-effects); Wide priors result in smaller marginal likelihoods, and thus models with wider priors are trivially less likely than models with narrower priors - where, at the extreme, that a model with completely flat priors is infinitely less favorable than a point null model (this is called the Jeffreys-Lindley-Bartlett paradox). Thus, you should only ever try (or want) to compute a Bayes factor when you have an informed prior.

Additionally, for models using MCMC estimation the number of posterior samples needed for testing is substantially larger than for estimation (the default of 4000 samples may not be enough in many cases). A conservative rule of thumb is to obtain 10 times more samples than would be required for estimation (Gronau, Singmann, & Wagenmakers, 2017). If less than 40,000 samples are detected, a warning is issued.

Transitivity of Bayes factors

For multiple inputs (models or hypotheses), the function will return multiple Bayes factors between each model and the same reference model (the denominator or un-restricted model). However, we can take advantage of the transitivity of Bayes factors - where if we have two Bayes factors for Model A and model B against the same reference model C, we can obtain a Bayes factor for comparing model A to model B by dividing them:

BF_{AB} = \frac{BF_{AC}}{BF_{BC}} = \frac{\frac{ML_{A}}{ML_{C}}}{\frac{ML_{B}}{ML_{C}}} = \frac{ML_{A}}{ML_{B}}

(Where ML is the marginal likelihood.)

A full matrix comparing all models can be obtained with as.matrix().

Interpreting Bayes Factors

A Bayes factor greater than 1 can be interpreted as evidence against the null, at which one convention is that a Bayes factor greater than 3 can be considered as "substantial" evidence against the null (and vice versa, a Bayes factor smaller than 1/3 indicates substantial evidence in favor of the null-model). See also effectsize::interpret_bf().

Note

There is also a plot()-method implemented in the see-package.

Author(s)

Mattan S. Ben-Shachar

References

• Gronau, Q. F., Singmann, H., & Wagenmakers, E. J. (2017). Bridgesampling: An R package for estimating normalizing constants. arXiv preprint arXiv:1710.08162.

• Kass, R. E., and Raftery, A. E. (1995). Bayes Factors. Journal of the American Statistical Association, 90(430), 773-795.

• Robert, C. P. (2016). The expected demise of the Bayes factor. Journal of Mathematical Psychology, 72, 33–37.

• Wagenmakers, E. J. (2007). A practical solution to the pervasive problems of p values. Psychonomic bulletin & review, 14(5), 779-804.

• Wetzels, R., Matzke, D., Lee, M. D., Rouder, J. N., Iverson, G. J., and Wagenmakers, E.-J. (2011). Statistical Evidence in Experimental Psychology: An Empirical Comparison Using 855 t Tests. Perspectives on Psychological Science, 6(3), 291–298. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/1745691611406923")}

See Also

bayesfactor_inclusion() for testing predictors across Bayesian models.

Other Bayes factors: bayesfactor_inclusion(), bayesfactor_parameters(), bayesfactor_restricted()

Examples

# With lm objects:
# ----------------
lm1 <- lm(mpg ~ 1, data = mtcars)
lm2 <- lm(mpg ~ hp, data = mtcars)
lm3 <- lm(mpg ~ hp + drat, data = mtcars)
lm4 <- lm(mpg ~ hp * drat, data = mtcars)
(BFM <- bayesfactor_models(lm1, lm2, lm3, lm4, denominator = 1))
# bayesfactor_models(lm2, lm3, lm4, denominator = lm1) # same result
# bayesfactor_models(lm1, lm2, lm3, lm4, denominator = lm1) # same result

update(BFM, reference = "bottom")
as.matrix(BFM)
as.numeric(BFM)

lm2b <- lm(sqrt(mpg) ~ hp, data = mtcars)
# Set check_response = TRUE for transformed responses
bayesfactor_models(lm2b, denominator = lm2, check_response = TRUE)


# With lmerMod objects:
# ---------------------
lmer1 <- lme4::lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
lmer2 <- lme4::lmer(Sepal.Length ~ Petal.Length + (Petal.Length | Species), data = iris)
lmer3 <- lme4::lmer(
  Sepal.Length ~ Petal.Length + (Petal.Length | Species) + (1 | Petal.Width),
  data = iris
)
bayesfactor_models(lmer1, lmer2, lmer3,
  denominator = 1,
  estimator = "REML"
)

# rstanarm models
# ---------------------
# (note that a unique diagnostic_file MUST be specified in order to work)
stan_m0 <- suppressWarnings(rstanarm::stan_glm(Sepal.Length ~ 1,
  data = iris,
  family = gaussian(),
  diagnostic_file = file.path(tempdir(), "df0.csv")
))
stan_m1 <- suppressWarnings(rstanarm::stan_glm(Sepal.Length ~ Species,
  data = iris,
  family = gaussian(),
  diagnostic_file = file.path(tempdir(), "df1.csv")
))
stan_m2 <- suppressWarnings(rstanarm::stan_glm(Sepal.Length ~ Species + Petal.Length,
  data = iris,
  family = gaussian(),
  diagnostic_file = file.path(tempdir(), "df2.csv")
))
bayesfactor_models(stan_m1, stan_m2, denominator = stan_m0, verbose = FALSE)


# brms models
# --------------------
# (note the save_pars MUST be set to save_pars(all = TRUE) in order to work)
brm1 <- brms::brm(Sepal.Length ~ 1, data = iris, save_pars = save_pars(all = TRUE))
brm2 <- brms::brm(Sepal.Length ~ Species, data = iris, save_pars = save_pars(all = TRUE))
brm3 <- brms::brm(
  Sepal.Length ~ Species + Petal.Length,
  data = iris,
  save_pars = save_pars(all = TRUE)
)

bayesfactor_models(brm1, brm2, brm3, denominator = 1, verbose = FALSE)


# BayesFactor
# ---------------------------
data(puzzles)
BF <- BayesFactor::anovaBF(RT ~ shape * color + ID,
  data = puzzles,
  whichRandom = "ID", progress = FALSE
)
BF
bayesfactor_models(BF) # basically the same

bayestestR documentation built on May 22, 2026, 1:06 a.m.