This function computes or extracts Bayes factors from fitted models.
The bf_*
function is an alias of the main function.
bayesfactor_models(..., denominator = 1, verbose = TRUE) bf_models(..., denominator = 1, verbose = TRUE) # S3 method for bayesfactor_models update(object, subset = NULL, reference = NULL, ...) # S3 method for bayesfactor_models as.matrix(x, ...)
...  Fitted models (see details), all fit on the same data, or a single


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

verbose  Toggle off warnings. 
object, x  A 
subset  Vector of model indices to keep or remove. 
reference  Index of model to reference to, or 
A data frame containing the models' formulas (reconstructed fixed and
random effects) and their log(BF)
s, that prints nicely.
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 followup
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 mostly a wraparound BayesFactor::extractBF()
.
BIC approximations are used to compute Bayes factors for all other model types (with a BIC method).
Note that BICs are extracted from models asis. So if for example you want to compare mixedmodels bases on ML instead of REML, you must supply models fit with ML.
In order to correctly and precisely estimate Bayes factors, a rule of thumb
are the 4 P's: Proper Priors and Plentiful
Posteriors. How many? 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, bayesfactor_models()
gives a warning.
See also the Bayes factors vignette.
There is also a plot()
method implemented in the seepackage.
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 nullmodel) (Wetzels et al. 2011).
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), 773795.
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), 779804.
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. doi: 10.1177/1745691611406923
Mattan S. BenShachar
# With lm objects: #  lm1 < lm(Sepal.Length ~ 1, data = iris) lm2 < lm(Sepal.Length ~ Species, data = iris) lm3 < lm(Sepal.Length ~ Species + Petal.Length, data = iris) lm4 < lm(Sepal.Length ~ Species * Petal.Length, data = iris) bayesfactor_models(lm1, lm2, lm3, lm4, denominator = 1)#> Bayes Factors for Model Comparison #> #> Model BF #> [lm2] Species 1.70e+29 #> [lm3] Species + Petal.Length 5.84e+55 #> [lm4] Species * Petal.Length 2.20e+54 #> #> * Against Denominator: [lm1] (Intercept only) #> * Bayes Factor Type: BIC approximationbayesfactor_models(lm2, lm3, lm4, denominator = lm1) # same result#> Bayes Factors for Model Comparison #> #> Model BF #> [lm2] Species 1.70e+29 #> [lm3] Species + Petal.Length 5.84e+55 #> [lm4] Species * Petal.Length 2.20e+54 #> #> * Against Denominator: [lm1] (Intercept only) #> * Bayes Factor Type: BIC approximationBFM < bayesfactor_models(lm1, lm2, lm3, lm4, denominator = lm1) # same result update(BFM, reference = "bottom")#> Bayes Factors for Model Comparison #> #> Model BF #> [lm2] Species 1.70e+29 #> [lm3] Species + Petal.Length 5.84e+55 #> [lm4] Species * Petal.Length 2.20e+54 #> #> * Against Denominator: [lm1] (Intercept only) #> * Bayes Factor Type: BIC approximation#> # Bayes Factors for Model Comparison #> #> Numerator #> Denominator #> #>  [1]  [2]  [3]  [4] #>  #> [1] (Intercept only)  1  1.70e+29  5.84e+55  2.20e+54 #> [2] Species  5.90e30  1  3.45e+26  1.30e+25 #> [3] Species + Petal.Length  1.71e56  2.90e27  1  0.038 #> [4] Species * Petal.Length  4.54e55  7.70e26  26.52  1if (FALSE) { # With lmerMod objects: #  if (require("lme4")) { lmer1 < lmer(Sepal.Length ~ Petal.Length + (1  Species), data = iris) lmer2 < lmer(Sepal.Length ~ Petal.Length + (Petal.Length  Species), data = iris) lmer3 < lmer( Sepal.Length ~ Petal.Length + (Petal.Length  Species) + (1  Petal.Width), data = iris ) bayesfactor_models(lmer1, lmer2, lmer3, denominator = 1) bayesfactor_models(lmer1, lmer2, lmer3, denominator = lmer1) } # rstanarm models #  # (note that a unique diagnostic_file MUST be specified in order to work) if (require("rstanarm")) { stan_m0 < stan_glm(Sepal.Length ~ 1, data = iris, family = gaussian(), diagnostic_file = file.path(tempdir(), "df0.csv") ) stan_m1 < stan_glm(Sepal.Length ~ Species, data = iris, family = gaussian(), diagnostic_file = file.path(tempdir(), "df1.csv") ) stan_m2 < 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) } # brms models #  # (note the save_pars MUST be set to save_pars(all = TRUE) in order to work) if (require("brms")) { brm1 < brm(Sepal.Length ~ 1, data = iris, save_all_pars = TRUE) brm2 < brm(Sepal.Length ~ Species, data = iris, save_all_pars = TRUE) brm3 < brm( Sepal.Length ~ Species + Petal.Length, data = iris, save_pars = save_pars(all = TRUE) ) bayesfactor_models(brm1, brm2, brm3, denominator = 1) } # BayesFactor #  if (require("BayesFactor")) { data(puzzles) BF < anovaBF(RT ~ shape * color + ID, data = puzzles, whichRandom = "ID", progress = FALSE ) BF bayesfactor_models(BF) # basically the same } }