bayesfactor_models: 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.

Usage

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, ...)

Arguments

...

Fitted models (see details), all fit on the same data, or a single BFBayesFactor object (see 'Details'). Ignored in as.matrix(), update().

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.

object, x

A bayesfactor_models object.

subset

Vector of model indices to keep or remove.

reference

Index of model to rereference to, or "top" to reference to the best model, or "bottom" to reference to the worst model.

Value

A data frame containing the models' formulas (reconstructed fixed and random effects) and their BFs, that prints nicely.

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) (Wetzels et al. 2011).

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 mostly a wraparound BayesFactor::extractBF().
For all other model types (supported by insight), BIC approximations are used to compute Bayes factors.

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.

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<U+2013>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<U+2013>298. 10.1177/1745691611406923

Examples

Run this code

# NOT RUN {
# 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)
bayesfactor_models(lm2, lm3, lm4, denominator = lm1) # same result
BFM <- bayesfactor_models(lm1, lm2, lm3, lm4, denominator = lm1) # same result

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

# }
# NOT RUN {
# 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
}
# }

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning