https://github.com/cran/bayestestR
Revision d8462ad2168ad7ee61c0d7e679174e775f01a9be authored by Dominique Makowski on 18 January 2020, 07:10:02 UTC, committed by cran-robot on 18 January 2020, 07:10:02 UTC
1 parent 4936034
Tip revision: d8462ad2168ad7ee61c0d7e679174e775f01a9be authored by Dominique Makowski on 18 January 2020, 07:10:02 UTC
version 0.5.0
version 0.5.0
Tip revision: d8462ad
bayesfactor_models.R
#' Bayes Factors (BF) for model comparison
#'
#' @description This function computes or extracts Bayes factors from fitted models.
#' \cr \cr
#' The \code{bf_*} function is an alias of the main function.
#'
#' @author Mattan S. Ben-Shachar
#'
#' @param ... Fitted models (see details), all fit on the same data, or a single \code{BFBayesFactor} object (see 'Details').
#' @param 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 \code{BFBayesFactor}.
#' @inheritParams hdi
#'
#' @details
#' If the passed models are supported by \pkg{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 \code{bayesfactor_inclusion}).
#'
#' \itemize{
#' \item For \code{brmsfit} or \code{stanreg} models, Bayes factors are computed using the \CRANpkg{bridgesampling} package.
#' \itemize{
#' \item \code{brmsfit} models must have been fitted with \code{save_all_pars = TRUE}.
#' \item \code{stanreg} models must have been fitted with a defined \code{diagnostic_file}.
#' }
#' \item For \code{BFBayesFactor}, \code{bayesfactor_models()} is mostly a wraparoud \code{BayesFactor::extractBF()}.
#' \item For all other model types (supported by \CRANpkg{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: \strong{P}roper \strong{P}riors and \strong{P}lentiful \strong{P}osterior
#' (i.e. probably at leat 40,000 samples instead of the default of 4,000).
#' \cr \cr
#' A Bayes factor greater than 1 can be interpereted as evidence against the compared-to
#' model (the denominator). One convention is that a Bayes factor greater than 3 can be considered
#' as "substantial" evidence against the denominator model (and vice versa, a Bayes factor
#' smaller than 1/3 indicates substantial evidence in favor of the denominator model)
#' (\cite{Wetzels et al. 2011}).
#' \cr \cr
#' See also \href{https://easystats.github.io/bayestestR/articles/bayes_factors.html}{the Bayes factors vignette}.
#'
#' @return A data frame containing the models' formulas (reconstructed fixed and random effects) and their BFs, that prints nicely.
#'
#' @examples
#' # 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
#' bayesfactor_models(lm1, lm2, lm3, lm4, denominator = lm1) # same result
#'
#' # With lmerMod objects:
#' # ---------------------
#' library(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)
#' \dontrun{
#' # rstanarm models
#' # ---------------------
#' # (note that a unique diagnostic_file MUST be specified in order to work)
#' library(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_all_pars MUST be set to TRUE in order to work)
#' library(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_all_pars = TRUE
#' )
#'
#' bayesfactor_models(brm1, brm2, brm3, denominator = 1)
#'
#'
#' # BayesFactor
#' # ---------------------------
#' library(BayesFactor)
#' data(puzzles)
#' BF <- anovaBF(RT ~ shape * color + ID,
#' data = puzzles,
#' whichRandom = "ID", progress = FALSE
#' )
#' BF
#' bayesfactor_models(BF) # basically the same
#' }
#'
#' @references
#' \itemize{
#' \item Gronau, Q. F., Wagenmakers, E. J., Heck, D. W., and Matzke, D. (2019). A simple method for comparing complex models: Bayesian model comparison for hierarchical multinomial processing tree models using Warp-III bridge sampling. Psychometrika, 84(1), 261-284.
#' \item Kass, R. E., and Raftery, A. E. (1995). Bayes Factors. Journal of the American Statistical Association, 90(430), 773-795.
#' \item Robert, C. P. (2016). The expected demise of the Bayes factor. Journal of Mathematical Psychology, 72, 33–37.
#' \item Wagenmakers, E. J. (2007). A practical solution to the pervasive problems of p values. Psychonomic bulletin & review, 14(5), 779-804.
#' \item 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}
#' }
#'
#' @importFrom insight get_response is_model_supported
#' @export
bayesfactor_models <- function(..., denominator = 1, verbose = TRUE) {
UseMethod("bayesfactor_models")
}
#' @rdname bayesfactor_models
#' @export
bf_models <- bayesfactor_models
#' @importFrom stats BIC
#' @export
bayesfactor_models.default <- function(..., denominator = 1, verbose = TRUE) {
# Organize the models and their names
mods <- list(...)
mnames <- sapply(match.call(expand.dots = F)$`...`, .safe_deparse)
if (!is.numeric(denominator)) {
model_name <- .safe_deparse(match.call()[["denominator"]])
denominator_model <- which(mnames == model_name)
if (length(denominator_model) == 0) {
mods <- c(mods, list(denominator))
mnames <- c(mnames, model_name)
denominator <- length(mods)
} else {
denominator <- denominator_model
}
}
# supported models
supported_models <- all(sapply(mods, insight::is_model_supported))
if (supported_models) {
# Test that all is good:
resps <- lapply(mods, insight::get_response)
if (!all(sapply(resps[-denominator], function(x) identical(x, resps[[denominator]])))) {
stop("Models were not computed from the same data.")
}
# Get formula
mforms <- sapply(mods, .find_full_formula)
} else {
if (verbose) {
object_names <- match.call(expand.dots = FALSE)$`...`
warning(sprintf(
"Unable to extract terms or validate that all models use the same data - following objects are not (supported) model objects (yet!): \n%s",
paste0(mnames[!supported_models], collapse = ", ")
), call. = FALSE)
}
mforms <- mnames
}
# Get BF
mBIC <- sapply(mods, BIC)
mBFs <- (mBIC - mBIC[denominator]) / (-2)
res <- data.frame(
Model = mforms,
BF = exp(mBFs),
stringsAsFactors = FALSE
)
attr(res, "denominator") <- denominator
attr(res, "BF_method") <- "BIC approximation"
attr(res, "unsupported_models") <- !supported_models
class(res) <- c("bayesfactor_models", "see_bayesfactor_models", class(res))
res
}
#' @importFrom insight get_response find_algorithm
.bayesfactor_models_stan <- function(..., denominator = 1, verbose = TRUE) {
if (!requireNamespace("bridgesampling")) {
stop("Package 'bridgesampling' required for this function to work. Please install it by running `install.packages('bridgesampling')`.")
}
# Orgenize the models
mods <- list(...)
# Warn
n_samps <- sapply(mods, function(x) {
alg <- insight::find_algorithm(x)
(alg$iterations - alg$warmup) * alg$chains
})
if (any(n_samps < 4e4)) {
warning(
"Bayes factors might not be precise.\n",
"For precise Bayes factors, it is recommended sampling at least 40,000 posterior samples."
)
}
if (!is.numeric(denominator)) {
model_name <- .safe_deparse(match.call()[["denominator"]])
arg_names <- sapply(match.call(expand.dots = F)$`...`, .safe_deparse)
denominator_model <- which(arg_names == model_name)
if (length(denominator_model) == 0) {
mods <- c(mods, list(denominator))
denominator <- length(mods)
} else {
denominator <- denominator_model
}
}
# Test that all is good:
resps <- lapply(mods, insight::get_response)
if (!all(sapply(resps[-denominator], function(x) identical(x, resps[[denominator]])))) {
stop("Models were not computed from the same data.")
}
# Get BF
if (verbose) {
message("Computation of Bayes factors: estimating marginal likelihood, please wait...")
}
mML <- lapply(mods, function(x) {
bridgesampling::bridge_sampler(x, silent = TRUE)
})
mBFs <- sapply(mML, function(x) {
bridgesampling::bf(x, mML[[denominator]], log = TRUE)[["bf"]]
})
# Get formula
mforms <- sapply(mods, .find_full_formula)
res <- data.frame(
Model = mforms,
BF = exp(mBFs),
stringsAsFactors = FALSE
)
attr(res, "denominator") <- denominator
attr(res, "BF_method") <- "marginal likelihoods (bridgesampling)"
attr(res, "unsupported_models") <- FALSE
class(res) <- c("bayesfactor_models", "see_bayesfactor_models", class(res))
res
}
#' @export
bayesfactor_models.stanreg <- function(..., denominator = 1, verbose = TRUE) {
if (!requireNamespace("rstanarm")) {
stop("Package 'rstanarm' required for this function to work. Please install it by running `install.packages('rstanarm')`.")
}
.bayesfactor_models_stan(..., denominator = denominator)
}
#' @export
bayesfactor_models.brmsfit <- function(..., denominator = 1, verbose = TRUE) {
if (!requireNamespace("brms")) {
stop("Package 'brms' required for this function to work. Please install it by running `install.packages('brms')`.")
}
if (!("brms" %in% .packages())) {
stop("This function requires package 'brms' to be loaded. Please run `library(brms)`.")
}
.bayesfactor_models_stan(..., denominator = denominator)
}
#' @export
bayesfactor_models.BFBayesFactor <- function(..., verbose = TRUE) {
models <- c(...)
if (!requireNamespace("BayesFactor")) {
stop("Package 'BayesFactor' required for this function to work. Please install it by running `install.packages('BayesFactor')`.")
}
mBFs <- c(0, BayesFactor::extractBF(models, TRUE, TRUE))
mforms <- sapply(c(models@denominator, models@numerator), function(x) x@shortName)
mforms[mforms == "Intercept only"] <- "1"
res <- data.frame(
Model = unname(mforms),
BF = exp(mBFs),
stringsAsFactors = FALSE
)
attr(res, "denominator") <- 1
attr(res, "BF_method") <- "JZS (BayesFactor)"
attr(res, "unsupported_models") <- !"BFlinearModel" %in% class(models@denominator)
class(res) <- c("bayesfactor_models", "see_bayesfactor_models", class(res))
res
}
#' @keywords internal
#' @importFrom insight find_formula
.find_full_formula <- function(mod) {
formulas <- insight::find_formula(mod)
conditional <- random <- NULL
if (!is.null(formulas$conditional)) {
conditional <- as.character(formulas$conditional)[3]
}
if (!is.null(formulas$random)) {
if (!is.list(formulas$random)) {
formulas$random <- list(formulas$random)
}
random <- sapply(formulas$random, function(x) {
paste0("(", as.character(x)[2], ")")
})
}
paste(c(conditional, random), collapse = " + ")
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...