https://github.com/cran/bayestestR
Tip revision: d73e4a2afcfbd6402c11716877e8f7466f309ef4 authored by Dominique Makowski on 22 October 2020, 13:40:02 UTC
version 0.7.5
version 0.7.5
Tip revision: d73e4a2
bayesfactor_parameters.R
#' Bayes Factors (BF) for a Single Parameter
#'
#' This method computes Bayes factors against the null (either a point or an interval),
#' based on prior and posterior samples of a single parameter. This Bayes factor indicates
#' the degree by which the mass of the posterior distribution has shifted further away
#' from or closer to the null value(s) (relative to the prior distribution), thus indicating
#' if the null value has become less or more likely given the observed data.
#' \cr \cr
#' When the null is an interval, the Bayes factor is computed by comparing the prior
#' and posterior odds of the parameter falling within or outside the null interval
#' (Morey & Rouder, 2011; Liao et al., 2020); When the null is a point, a Savage-Dickey
#' density ratio is computed, which is also an approximation of a Bayes factor comparing
#' the marginal likelihoods of the model against a model in which the tested parameter
#' has been restricted to the point null (Wagenmakers et al., 2010; Heck, 2019).
#' \cr \cr
#' Note that the \code{logspline} package is used for estimating densities and probabilities,
#' and must be installed for the function to work.
#' \cr \cr
#' \code{bayesfactor_pointnull()} and \code{bayesfactor_rope()} are wrappers around
#' \code{bayesfactor_parameters} with different defaults for the null to be tested against
#' (a point and a range, respectively). Aliases of the main functions are prefixed
#' with \code{bf_*}, like \code{bf_parameters()} or \code{bf_pointnull()}
#' \cr \cr
#' \strong{For more info, in particular 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}.}
#'
#' @param posterior A numerical vector, \code{stanreg} / \code{brmsfit} object, \code{emmGrid}
#' or a data frame - representing a posterior distribution(s) from (see 'Details').
#' @param prior An object representing a prior distribution (see 'Details').
#' @param direction Test type (see 'Details'). One of \code{0}, \code{"two-sided"} (default, two tailed),
#' \code{-1}, \code{"left"} (left tailed) or \code{1}, \code{"right"} (right tailed).
#' @param null Value of the null, either a scalar (for point-null) or a range
#' (for a interval-null).
#' @param ... Arguments passed to and from other methods.
#' (Can be used to pass arguments to internal \code{\link[logspline]{logspline}}.)
#' @inheritParams hdi
#'
#' @return A data frame containing the Bayes factor representing evidence \emph{against} the null.
#'
#' @note There is also a \href{https://easystats.github.io/see/articles/bayestestR.html}{\code{plot()}-method} implemented in the \href{https://easystats.github.io/see/}{\pkg{see}-package}.
#'
#' @details This method is used to compute Bayes factors based on prior and posterior distributions.
#' \cr\cr
#' For the computation of Bayes factors, the model priors must be proper priors (at the very least
#' they should be \emph{not flat}, and it is preferable that they be \emph{informative}); As the priors for
#' the alternative get wider, the likelihood of the null value(s) increases, to the extreme that for completely
#' flat priors the null is infinitely more favorable than the alternative (this is called \emph{the Jeffreys-Lindley-Bartlett
#' paradox}). Thus, you should only ever try (or want) to compute a Bayes factor when you have an informed prior.
#' \cr\cr
#' (Note that by default, \code{brms::brm()} uses flat priors for fixed-effects; 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 numerical vector, \code{prior} should also be a numerical vector.
#' \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 equivalent 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 equivalent to \code{posterior} but created with a model of priors samples \emph{only}.
#' \item \strong{Note:} When the \code{emmGrid} has undergone any transformations (\code{"log"}, \code{"response"}, etc.), or \code{regrid}ing, then \code{prior} must be an \code{emmGrid} object, as stated above.
#' }
#' }}
#' \subsection{One-sided Tests (setting an order restriction)}{
#' One sided tests (controlled by \code{direction}) are conducted by restricting the prior and
#' posterior of the non-null values (the "alternative") to one side of the null only
#' (\cite{Morey & Wagenmakers, 2014}). For example, if we have a prior hypothesis that the
#' parameter should be positive, the alternative will be restricted to the region to the right
#' of the null (point or interval).
#' }
#' \subsection{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)
#' (\cite{Wetzels et al. 2011}).
#' }
#'
#' @examples
#' library(bayestestR)
#'
#' prior <- distribution_normal(1000, mean = 0, sd = 1)
#' posterior <- distribution_normal(1000, mean = .5, sd = .3)
#'
#' bayesfactor_parameters(posterior, prior)
#' \dontrun{
#' # rstanarm models
#' # ---------------
#' if (require("rstanarm") && require("emmeans")) {
#' contrasts(sleep$group) <- contr.bayes # see vingette
#' stan_model <- stan_lmer(extra ~ group + (1 | ID), data = sleep)
#' bayesfactor_parameters(stan_model)
#' bayesfactor_parameters(stan_model, null = rope_range(stan_model))
#'
#' # emmGrid objects
#' # ---------------
#' group_diff <- pairs(emmeans(stan_model, ~group))
#' bayesfactor_parameters(group_diff, prior = stan_model)
#' }
#'
#' # brms models
#' # -----------
#' if (require("brms")) {
#' contrasts(sleep$group) <- contr.bayes # see vingette
#' my_custom_priors <-
#' set_prior("student_t(3, 0, 1)", class = "b") +
#' set_prior("student_t(3, 0, 1)", class = "sd", group = "ID")
#'
#' brms_model <- brm(extra ~ group + (1 | ID),
#' data = sleep,
#' prior = my_custom_priors
#' )
#' bayesfactor_parameters(brms_model)
#' }
#' }
#' @references
#' \itemize{
#' \item Wagenmakers, E. J., Lodewyckx, T., Kuriyal, H., and Grasman, R. (2010). Bayesian hypothesis testing for psychologists: A tutorial on the Savage-Dickey method. Cognitive psychology, 60(3), 158-189.
#' \item Heck, D. W. (2019). A caveat on the Savage–Dickey density ratio: The case of computing Bayes factors for regression parameters. British Journal of Mathematical and Statistical Psychology, 72(2), 316-333.
#' \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 Liao, J. G., Midya, V., & Berg, A. (2020). Connecting and contrasting the Bayes factor and a modified ROPE procedure for testing interval null hypotheses. The American Statistician, 1-19.
#' \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}
#' }
#'
#' @author Mattan S. Ben-Shachar
#'
#' @export
bayesfactor_parameters <- function(posterior, prior = NULL, direction = "two-sided", null = 0, verbose = TRUE, ...) {
UseMethod("bayesfactor_parameters")
}
#' @rdname bayesfactor_parameters
#' @export
bayesfactor_pointull <- function(posterior, prior = NULL, direction = "two-sided", null = 0, verbose = TRUE, ...) {
if (length(null) > 1) {
message("'null' is a range - computing a ROPE based Bayes factor.")
}
bayesfactor_parameters(
posterior = posterior,
prior = prior,
direction = direction,
null = null,
verbose = verbose,
...
)
}
#' @rdname bayesfactor_parameters
#' @export
bayesfactor_rope <- function(posterior, prior = NULL, direction = "two-sided", null = rope_range(posterior), verbose = TRUE, ...) {
if (length(null) < 2) {
message("'null' is a point - computing a Savage-Dickey (point null) Bayes factor.")
}
bayesfactor_parameters(
posterior = posterior,
prior = prior,
direction = direction,
null = null,
verbose = verbose,
...
)
}
#' @rdname bayesfactor_parameters
#' @export
bf_parameters <- bayesfactor_parameters
#' @rdname bayesfactor_parameters
#' @export
bf_pointull <- bayesfactor_pointull
#' @rdname bayesfactor_parameters
#' @export
bf_rope <- bayesfactor_rope
#' @rdname bayesfactor_parameters
#' @export
bayesfactor_parameters.numeric <- function(posterior, prior = NULL, direction = "two-sided", null = 0, verbose = TRUE, ...) {
# nm <- .safe_deparse(substitute(posterior)
if (is.null(prior)) {
prior <- posterior
if (verbose) {
warning(
"Prior not specified! ",
"Please specify a prior (in the form 'prior = distribution_normal(1000, 0, 1)')",
" to get meaningful results."
)
}
}
prior <- data.frame(X = prior)
posterior <- data.frame(X = posterior)
# colnames(posterior) <- colnames(prior) <- nm
# Get BFs
sdbf <- bayesfactor_parameters.data.frame(
posterior = posterior, prior = prior,
direction = direction, null = null, ...
)
sdbf$Parameter <- NULL
sdbf
}
#' @importFrom insight clean_parameters
#' @rdname bayesfactor_parameters
#' @export
bayesfactor_parameters.stanreg <- function(posterior,
prior = NULL,
direction = "two-sided",
null = 0,
verbose = TRUE,
effects = c("fixed", "random", "all"),
component = c("conditional", "zi", "zero_inflated", "all"),
parameters = NULL,
...) {
cleaned_parameters <- insight::clean_parameters(posterior)
effects <- match.arg(effects)
component <- match.arg(component)
samps <- .clean_priors_and_posteriors(posterior, prior,
verbose = verbose,
effects = effects, component = component,
parameters = parameters)
# Get BFs
temp <- bayesfactor_parameters.data.frame(
posterior = samps$posterior, prior = samps$prior,
direction = direction, null = null, ...
)
bf_val <- .prepare_output(temp, cleaned_parameters, inherits(posterior, "stanmvreg"))
class(bf_val) <- class(temp)
attr(bf_val, "hypothesis") <- attr(temp, "hypothesis") # don't change the name of this attribute - it is used only internally for "see" and printing
attr(bf_val, "direction") <- attr(temp, "direction")
attr(bf_val, "plot_data") <- attr(temp, "plot_data")
bf_val
}
#' @rdname bayesfactor_parameters
#' @export
bayesfactor_parameters.brmsfit <- bayesfactor_parameters.stanreg
#' @rdname bayesfactor_parameters
#' @export
bayesfactor_parameters.emmGrid <- function(posterior,
prior = NULL,
direction = "two-sided",
null = 0,
verbose = TRUE,
...) {
samps <- .clean_priors_and_posteriors(posterior, prior,
verbose = verbose)
# Get BFs
bayesfactor_parameters.data.frame(
posterior = samps$posterior, prior = samps$prior,
direction = direction, null = null, ...
)
}
#' @export
bayesfactor_parameters.emm_list <- bayesfactor_parameters.emmGrid
#' @rdname bayesfactor_parameters
#' @export
bayesfactor_parameters.data.frame <- function(posterior,
prior = NULL,
direction = "two-sided",
null = 0,
verbose = TRUE,
...) {
# find direction
direction <- .get_direction(direction)
if (is.null(prior)) {
prior <- posterior
warning(
"Prior not specified! ",
"Please specify priors (with column order matching 'posterior')",
" to get meaningful results."
)
}
sdbf <- numeric(ncol(posterior))
for (par in seq_along(posterior)) {
sdbf[par] <- .bayesfactor_parameters(
posterior[[par]],
prior[[par]],
direction = direction,
null = null,
...
)
}
bf_val <- data.frame(
Parameter = colnames(posterior),
BF = sdbf,
stringsAsFactors = FALSE
)
class(bf_val) <- unique(c(
"bayesfactor_parameters",
"see_bayesfactor_parameters",
class(bf_val)
))
attr(bf_val, "hypothesis") <- null # don't change the name of this attribute - it is used only internally for "see" and printing
attr(bf_val, "direction") <- direction
attr(bf_val, "plot_data") <- .make_BF_plot_data(posterior, prior, direction, null, ...)
bf_val
}
#' @keywords internal
#' @importFrom insight print_color
.bayesfactor_parameters <- function(posterior, prior, direction = 0, null = 0, ...) {
if (isTRUE(all.equal(posterior, prior))) {
return(1)
}
if (!requireNamespace("logspline")) {
stop("Package \"logspline\" needed for this function to work. Please install it.")
}
if (length(null) == 1) {
relative_density <- function(samples) {
f_samples <- .logspline(samples, ...)
d_samples <- logspline::dlogspline(null, f_samples)
if (direction < 0) {
norm_samples <- logspline::plogspline(null, f_samples)
} else if (direction > 0) {
norm_samples <- 1 - logspline::plogspline(null, f_samples)
} else {
norm_samples <- 1
}
d_samples / norm_samples
}
return(relative_density(prior) /
relative_density(posterior))
} else if (length(null) == 2) {
null <- sort(null)
null[is.infinite(null)] <- 1.797693e+308 * sign(null[is.infinite(null)])
f_prior <- .logspline(prior, ...)
f_posterior <- .logspline(posterior, ...)
h0_prior <- diff(logspline::plogspline(null, f_prior))
h0_post <- diff(logspline::plogspline(null, f_posterior))
BF_null_full <- h0_post / h0_prior
if (direction < 0) {
h1_prior <- logspline::plogspline(min(null), f_prior)
h1_post <- logspline::plogspline(min(null), f_posterior)
} else if (direction > 0) {
h1_prior <- 1 - logspline::plogspline(max(null), f_prior)
h1_post <- 1 - logspline::plogspline(max(null), f_posterior)
} else {
h1_prior <- 1 - h0_prior
h1_post <- 1 - h0_post
}
BF_alt_full <- h1_post / h1_prior
return(BF_alt_full / BF_null_full)
} else {
stop("'null' must be of length 1 or 2")
}
}
# Bad Methods -------------------------------------------------------------
#' @export
bayesfactor_parameters.bayesfactor_models <- function(...) {
stop(
"Oh no, 'bayesfactor_parameters()' does not know how to deal with multiple models :(\n",
"You might want to use 'bayesfactor_inclusion()' here to test specific terms across models."
)
}
#' @export
bayesfactor_parameters.sim <- function(...) {
stop(
"Bayes factors are based on the shift from a prior to a posterior. ",
"Since simulated draws are not based on any priors, computing Bayes factors does not make sense :(\n",
"You might want to try `rope`, `ci`, `pd` or `pmap` for posterior-based inference."
)
}
#' @export
bayesfactor_parameters.sim.merMod <- bayesfactor_parameters.sim