#' Bayes Factors (BF) for Order Restricted Models #' #' This method computes Bayes factors for comparing a model with an order restrictions on its parameters #' with the fully unrestricted model. \emph{Note that this method should only be used for confirmatory analyses}. #' \cr \cr #' \strong{For info on specifying correct priors for factors with more than 2 levels, see \href{https://easystats.github.io/bayestestR/articles/bayes_factors.html}{the Bayes factors vignette}.} #' \cr \cr #' For more info, see \href{https://easystats.github.io/bayestestR/articles/bayes_factors.html}{the Bayes factors vignette}. #' #' @param posterior A \code{stanreg} / \code{brmsfit} object, \code{emmGrid} or a data frame - representing a posterior distribution(s) from (see Details). #' @param hypothesis A character vector specifying the restrictions as logical conditions (see examples below). #' @param prior An object representing a prior distribution (see Details). #' @inheritParams hdi #' #' @details This method is used to compute Bayes factors for order-restricted models vs un-restricted #' models by setting an order restriction on the prior and posterior distributions #' (\cite{Morey & Wagenmakers, 2013}). #' #' (Though it is possible to use \code{bayesfactor_restricted} to test interval restrictions, #' it is more suitable for testing order restrictions (see examples)). #' #' When \code{posterior} is a model (\code{stanreg}, \code{brmsfit}), posterior and prior samples are #' extracted for each parameter, and Savage-Dickey Bayes factors are computed for each parameter. #' #' \strong{NOTE:} For \code{brmsfit} models, the model must have been fitted with \emph{custom (non-default)} priors. See example below. #' #' \subsection{Setting the correct \code{prior}}{ #' It is important to provide the correct \code{prior} for meaningful results. #' \itemize{ #' \item When \code{posterior} is a \code{data.frame}, \code{prior} should also be a \code{data.frame}, with matching column order. #' \item When \code{posterior} is a \code{stanreg} or \code{brmsfit} model: \itemize{ #' \item \code{prior} can be set to \code{NULL}, in which case prior samples are drawn internally. #' \item \code{prior} can also be a model equvilant to \code{posterior} but with samples from the priors \emph{only}. #' } #' \item When \code{posterior} is an \code{emmGrid} object: \itemize{ #' \item \code{prior} should be the \code{stanreg} or \code{brmsfit} model used to create the \code{emmGrid} objects. #' \item \code{prior} can also be an \code{emmGrid} object equvilant to \code{posterior} but created with a model of priors samples \emph{only}. #' } #' }} #' \subsection{Interpreting Bayes Factors}{ #' A Bayes factor greater than 1 can be interpereted 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-hypothesis) #' (\cite{Wetzels et al. 2011}). #' } #' #' @return A data frame containing the Bayes factor representing evidence \emph{against} the un-restricted model. #' #' @examples #' library(bayestestR) #' prior <- data.frame( #' X = rnorm(100), #' X1 = rnorm(100), #' X3 = rnorm(100) #' ) #' #' posterior <- data.frame( #' X = rnorm(100, .4), #' X1 = rnorm(100, -.2), #' X3 = rnorm(100) #' ) #' #' hyps <- c( #' "X > X1 & X1 > X3", #' "X > X1" #' ) #' #' bayesfactor_restricted(posterior, hypothesis = hyps, prior = prior) #' \dontrun{ #' # rstanarm models #' # --------------- #' library(rstanarm) #' fit_stan <- stan_glm(mpg ~ wt + cyl + am, #' data = mtcars #' ) #' hyps <- c( #' "am > 0 & cyl < 0", #' "cyl < 0", #' "wt - cyl > 0" #' ) #' bayesfactor_restricted(fit_stan, hypothesis = hyps) #' #' # emmGrid objects #' # --------------- #' library(emmeans) #' #' # replicating http://bayesfactor.blogspot.com/2015/01/multiple-comparisons-with-bayesfactor-2.html #' disgust_data <- read.table(url("http://www.learnbayes.org/disgust_example.txt"), header = TRUE) #' #' contrasts(disgust_data\$condition) <- contr.bayes # see vignette #' fit_model <- stan_glm(score ~ condition, data = disgust_data, family = gaussian()) #' #' em_condition <- emmeans(fit_model, ~condition) #' hyps <- c("lemon < control & control < sulfur") #' #' bayesfactor_restricted(em_condition, prior = fit_model, hypothesis = hyps) #' # > # Bayes Factor (Order-Restriction) #' # > #' # > Hypothesis P(Prior) P(Posterior) Bayes Factor #' # > lemon < control & control < sulfur 0.17 0.75 4.49 #' # > --- #' # > Bayes factors for the restricted movel vs. the un-restricted model. #' } #' #' @references #' \itemize{ #' \item Morey, R. D., & Wagenmakers, E. J. (2014). Simple relation between Bayesian order-restricted and point-null hypothesis tests. Statistics & Probability Letters, 92, 121-124. #' \item Morey, R. D., & Rouder, J. N. (2011). Bayes factor approaches for testing interval null hypotheses. Psychological methods, 16(4), 406. #' \item Morey, R. D. (Jan, 2015). Multiple Comparisons with BayesFactor, Part 2 – order restrictions. Retrived from https://richarddmorey.org/category/order-restrictions/. #' } #' #' @export bayesfactor_restricted <- function(posterior, hypothesis, prior = NULL, verbose = TRUE, ...) { UseMethod("bayesfactor_restricted") } #' @importFrom insight get_parameters #' @rdname bayesfactor_restricted #' @export bayesfactor_restricted.stanreg <- function(posterior, hypothesis, prior = NULL, verbose = TRUE, effects = c("fixed", "random", "all"), ...) { effects <- match.arg(effects) # Get Priors if (is.null(prior)) { prior <- posterior } prior <- .update_to_priors(prior, verbose = verbose) prior <- insight::get_parameters(prior, effects = effects) posterior <- insight::get_parameters(posterior, effects = effects) # Get savage-dickey BFs bayesfactor_restricted.data.frame( posterior = posterior, prior = prior, hypothesis = hypothesis ) } #' @rdname bayesfactor_restricted #' @export bayesfactor_restricted.brmsfit <- bayesfactor_restricted.stanreg #' @importFrom stats update #' @importFrom insight get_parameters #' @rdname bayesfactor_restricted #' @export bayesfactor_restricted.emmGrid <- function(posterior, hypothesis, prior = NULL, verbose = TRUE, ...) { if (!requireNamespace("emmeans")) { stop("Package \"emmeans\" needed for this function to work. Please install it.") } if (is.null(prior)) { prior <- posterior warning( "Prior not specified! ", "Please provide the original model to get meaningful results." ) } else if (!inherits(prior, "emmGrid")) { # then is it a model prior <- .update_to_priors(prior, verbose = verbose) prior <- insight::get_parameters(prior, effects = "fixed") prior <- stats::update(posterior, post.beta = as.matrix(prior)) } prior <- as.data.frame(as.matrix(emmeans::as.mcmc.emmGrid(prior, names = FALSE))) posterior <- as.data.frame(as.matrix(emmeans::as.mcmc.emmGrid(posterior, names = FALSE))) bayesfactor_restricted.data.frame( posterior = posterior, prior = prior, hypothesis = hypothesis ) } #' @export bayesfactor_restricted.data.frame <- function(posterior, hypothesis, prior = NULL, ...) { p_hypothesis <- parse(text = hypothesis) if (is.null(prior)) { prior <- posterior warning( "Prior not specified! ", "Please specify priors (with column names matching 'posterior')", " to get meaningful results." ) } .get_prob <- function(x, data) { x_logical <- try(eval(x, envir = data), silent = TRUE) if (inherits(x_logical, "try-error")) { cnames <- colnames(data) is_name <- make.names(cnames) == cnames cnames[!is_name] <- paste0("", cnames[!is_name], "") cnames <- paste0(cnames, collapse = ", ") stop(x_logical, "Available parameters are: ", cnames) } else if (!all(is.logical(x_logical))) { stop("Hypotheses must be logical") } mean(x_logical) } posterior_p <- sapply(p_hypothesis, .get_prob, data = posterior) prior_p <- sapply(p_hypothesis, .get_prob, data = prior) BF <- posterior_p / prior_p res <- data.frame( Hypothesis = hypothesis, Prior_prob = prior_p, Posterior_prob = posterior_p, BF = BF ) class(res) <- unique(c( "bayesfactor_restricted", class(res) )) res }