Package 'bayestestR'

Description

Based on the DescTools AUC function. It can calculate the area under the curve with a naive algorithm or a more elaborated spline approach. The curve must be given by vectors of xy-coordinates. This function can handle unsorted x values (by sorting x) and ties for the x values (by ignoring duplicates).

Usage

area_under_curve(x, y, method = c("trapezoid", "step", "spline"), ...)

auc(x, y, method = c("trapezoid", "step", "spline"), ...)

Arguments

See Also

DescTools

Examples

library(bayestestR)
posterior <- distribution_normal(1000)

dens <- estimate_density(posterior)
dens <- dens[dens$x > 0, ]
x <- dens$x
y <- dens$y

area_under_curve(x, y, method = "trapezoid")
area_under_curve(x, y, method = "step")
area_under_curve(x, y, method = "spline")

Description

Coerce to a Data Frame

Usage

## S3 method for class 'density'
as.data.frame(x, ...)

Arguments

Description

Convert to Numeric

Usage

## S3 method for class 'map_estimate'
as.numeric(x, ...)

## S3 method for class 'p_direction'
as.numeric(x, ...)

## S3 method for class 'p_map'
as.numeric(x, ...)

## S3 method for class 'p_significance'
as.numeric(x, ...)

Arguments

Description

This function compte the Bayes factors (BFs) that are appropriate to the input. For vectors or single models, it will compute BFs for single parameters(), or is hypothesis is specified, BFs for restricted models(). For multiple models, it will return the BF corresponding to comparison between models() and if a model comparison is passed, it will compute the inclusion BF().

For a complete overview of these functions, read the Bayes factor vignette.

Usage

bayesfactor(
  ...,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  hypothesis = NULL,
  effects = c("fixed", "random", "all"),
  verbose = TRUE,
  denominator = 1,
  match_models = FALSE,
  prior_odds = NULL
)

Arguments

Value

Some type of Bayes factor, depending on the input. See bayesfactor_parameters(), bayesfactor_models() or bayesfactor_inclusion()

Note

There is also a plot()-method implemented in the see-package.

Examples

library(bayestestR)

prior <- distribution_normal(1000, mean = 0, sd = 1)
posterior <- distribution_normal(1000, mean = .5, sd = .3)

bayesfactor(posterior, prior = prior, verbose = FALSE)


# rstanarm models
# ---------------
model <- suppressWarnings(rstanarm::stan_lmer(extra ~ group + (1 | ID), data = sleep))
bayesfactor(model, verbose = FALSE)

# Frequentist models
# ---------------
m0 <- lm(extra ~ 1, data = sleep)
m1 <- lm(extra ~ group, data = sleep)
m2 <- lm(extra ~ group + ID, data = sleep)

comparison <- bayesfactor(m0, m1, m2)
comparison

bayesfactor(comparison)

Description

The ⁠bf_*⁠ function is an alias of the main function.

For more info, see the Bayes factors vignette.

Usage

bayesfactor_inclusion(models, match_models = FALSE, prior_odds = NULL, ...)

bf_inclusion(models, match_models = FALSE, prior_odds = NULL, ...)

Arguments

Details

Inclusion Bayes factors answer the question: Are the observed data more probable under models with a particular effect, than they are under models without that particular effect? In other words, on average - are models with effect $X$ more likely to have produced the observed data than models without effect $X$?

Match Models

If match_models=FALSE (default), Inclusion BFs are computed by comparing all models with a term against all models without that term. If TRUE, comparison is restricted to models that (1) do not include any interactions with the term of interest; (2) for interaction terms, averaging is done only across models that containe the main effect terms from which the interaction term is comprised.

Value

a data frame containing the prior and posterior probabilities, and log(BF) for each effect (Use as.numeric() to extract the non-log Bayes factors; see examples).

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

Note

Random effects in the lmer style are converted to interaction terms: i.e., (X|G) will become the terms 1:G and X:G.

Author(s)

Mattan S. Ben-Shachar

References

• Hinne, M., Gronau, Q. F., van den Bergh, D., and Wagenmakers, E. (2019, March 25). A conceptual introduction to Bayesian Model Averaging. doi:10.31234/osf.io/wgb64

• Clyde, M. A., Ghosh, J., & Littman, M. L. (2011). Bayesian adaptive sampling for variable selection and model averaging. Journal of Computational and Graphical Statistics, 20(1), 80-101.

• Mathot, S. (2017). Bayes like a Baws: Interpreting Bayesian Repeated Measures in JASP. Blog post.

See Also

weighted_posteriors() for Bayesian parameter averaging.

Examples

library(bayestestR)

# Using bayesfactor_models:
# ------------------------------
mo0 <- lm(Sepal.Length ~ 1, data = iris)
mo1 <- lm(Sepal.Length ~ Species, data = iris)
mo2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)
mo3 <- lm(Sepal.Length ~ Species * Petal.Length, data = iris)

BFmodels <- bayesfactor_models(mo1, mo2, mo3, denominator = mo0)
(bf_inc <- bayesfactor_inclusion(BFmodels))

as.numeric(bf_inc)


# BayesFactor
# -------------------------------
BF <- BayesFactor::generalTestBF(len ~ supp * dose, ToothGrowth, progress = FALSE)
bayesfactor_inclusion(BF)

# compare only matched models:
bayesfactor_inclusion(BF, match_models = TRUE)

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)

## Default S3 method:
bayesfactor_models(..., denominator = 1, verbose = TRUE)

## S3 method for class 'bayesfactor_models'
update(object, subset = NULL, reference = NULL, ...)

## S3 method for class 'bayesfactor_models'
as.matrix(x, ...)

Arguments

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, Bayes factors are computed using the BIC approximation. Note that BICs are extracted from using insight::get_loglikelihood, see documentation there for options for dealing with transformed responses and REML estimation.

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.

Value

A data frame containing the models' formulas (reconstructed fixed and random effects) and their log(BF)s (Use as.numeric() to extract the non-log Bayes factors; see examples), 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).

Note

There is also a plot()-method implemented in the see-package.

Author(s)

Mattan S. Ben-Shachar

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–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–298. doi:10.1177/1745691611406923

Examples

# With lm objects:
# ----------------
lm1 <- lm(mpg ~ 1, data = mtcars)
lm2 <- lm(mpg ~ hp, data = mtcars)
lm3 <- lm(mpg ~ hp + drat, data = mtcars)
lm4 <- lm(mpg ~ hp * drat, data = mtcars)
(BFM <- 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

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

lm2b <- lm(sqrt(mpg) ~ hp, data = mtcars)
# Set check_response = TRUE for transformed responses
bayesfactor_models(lm2b, denominator = lm2, check_response = TRUE)


# With lmerMod objects:
# ---------------------
lmer1 <- lme4::lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
lmer2 <- lme4::lmer(Sepal.Length ~ Petal.Length + (Petal.Length | Species), data = iris)
lmer3 <- lme4::lmer(
  Sepal.Length ~ Petal.Length + (Petal.Length | Species) + (1 | Petal.Width),
  data = iris
)
bayesfactor_models(lmer1, lmer2, lmer3,
  denominator = 1,
  estimator = "REML"
)

# rstanarm models
# ---------------------
# (note that a unique diagnostic_file MUST be specified in order to work)
stan_m0 <- suppressWarnings(rstanarm::stan_glm(Sepal.Length ~ 1,
  data = iris,
  family = gaussian(),
  diagnostic_file = file.path(tempdir(), "df0.csv")
))
stan_m1 <- suppressWarnings(rstanarm::stan_glm(Sepal.Length ~ Species,
  data = iris,
  family = gaussian(),
  diagnostic_file = file.path(tempdir(), "df1.csv")
))
stan_m2 <- suppressWarnings(rstanarm::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, verbose = FALSE)


# brms models
# --------------------
# (note the save_pars MUST be set to save_pars(all = TRUE) in order to work)
brm1 <- brms::brm(Sepal.Length ~ 1, data = iris, save_pars = save_pars(all = TRUE))
brm2 <- brms::brm(Sepal.Length ~ Species, data = iris, save_pars = save_pars(all = TRUE))
brm3 <- brms::brm(
  Sepal.Length ~ Species + Petal.Length,
  data = iris,
  save_pars = save_pars(all = TRUE)
)

bayesfactor_models(brm1, brm2, brm3, denominator = 1, verbose = FALSE)


# BayesFactor
# ---------------------------
data(puzzles)
BF <- BayesFactor::anovaBF(RT ~ shape * color + ID,
  data = puzzles,
  whichRandom = "ID", progress = FALSE
)
BF
bayesfactor_models(BF) # basically the same

Description

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.

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

Note that the logspline package is used for estimating densities and probabilities, and must be installed for the function to work.

bayesfactor_pointnull() and bayesfactor_rope() are wrappers around 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 ⁠bf_*⁠, like bf_parameters() or bf_pointnull().

For more info, in particular on specifying correct priors for factors with more than 2 levels, see the Bayes factors vignette.

Usage

bayesfactor_parameters(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  ...,
  verbose = TRUE
)

bayesfactor_pointnull(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  ...,
  verbose = TRUE
)

bayesfactor_rope(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = rope_range(posterior, verbose = FALSE),
  ...,
  verbose = TRUE
)

bf_parameters(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  ...,
  verbose = TRUE
)

bf_pointnull(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  ...,
  verbose = TRUE
)

bf_rope(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = rope_range(posterior, verbose = FALSE),
  ...,
  verbose = TRUE
)

## S3 method for class 'numeric'
bayesfactor_parameters(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  ...,
  verbose = TRUE
)

## S3 method for class 'stanreg'
bayesfactor_parameters(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "location", "smooth_terms", "sigma", "zi",
    "zero_inflated", "all"),
  parameters = NULL,
  ...,
  verbose = TRUE
)

## S3 method for class 'brmsfit'
bayesfactor_parameters(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "location", "smooth_terms", "sigma", "zi",
    "zero_inflated", "all"),
  parameters = NULL,
  ...,
  verbose = TRUE
)

## S3 method for class 'blavaan'
bayesfactor_parameters(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  ...,
  verbose = TRUE
)

## S3 method for class 'data.frame'
bayesfactor_parameters(
  posterior,
  prior = NULL,
  direction = "two-sided",
  null = 0,
  rvar_col = NULL,
  ...,
  verbose = TRUE
)

Arguments

Details

This method is used to compute Bayes factors based on prior and posterior distributions.

One-sided & Dividing Tests (setting an order restriction)

One sided tests (controlled by direction) are conducted by restricting the prior and posterior of the non-null values (the "alternative") to one side of the null only (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). For example, for a Bayes factor comparing the "null" of 0-0.1 to the alternative ⁠>0.1⁠, we would set bayesfactor_parameters(null = c(0, 0.1), direction = ">").

It is also possible to compute a Bayes factor for dividing hypotheses - that is, for a null and alternative that are complementary, opposing one-sided hypotheses (Morey & Wagenmakers, 2014). For example, for a Bayes factor comparing the "null" of ⁠<0⁠ to the alternative ⁠>0⁠, we would set bayesfactor_parameters(null = c(-Inf, 0)).

Value

A data frame containing the (log) Bayes factor representing evidence against the null (Use as.numeric() to extract the non-log Bayes factors; see examples).

Setting the correct prior

For the computation of Bayes factors, the model priors must be proper priors (at the very least they should be not flat, and it is preferable that they be 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 the Jeffreys-Lindley-Bartlett paradox). Thus, you should only ever try (or want) to compute a Bayes factor when you have an informed prior.

(Note that by default, brms::brm() uses flat priors for fixed-effects; See example below.)

It is important to provide the correct prior for meaningful results, to match the posterior-type input:

• A numeric vector - prior should also be a numeric vector, representing the prior-estimate.

• A data frame - prior should also be a data frame, representing the prior-estimates, in matching column order.

  • If rvar_col is specified, prior should be the name of an rvar column that represents the prior-estimates.

• Supported Bayesian model (stanreg, brmsfit, etc.)

  • prior should be a model an equivalent model with MCMC samples from the priors only. See unupdate().

  • If prior is set to NULL, unupdate() is called internally (not supported for brmsfit_multiple model).

• Output from a {marginaleffects} function - prior should also be an equivalent output from a {marginaleffects} function based on a prior-model (See unupdate()).

• Output from an {emmeans} function

  • prior should also be an equivalent output from an {emmeans} function based on a prior-model (See unupdate()).

  • prior can also be the original (posterior) model, in which case the function will try to "unupdate" the estimates (not supported if the estimates have undergone any transformations – "log", "response", etc. – or any regriding).

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

Note

There is also a plot()-method implemented in the see-package.

Author(s)

Mattan S. Ben-Shachar

References

• 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.

• 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.

• Morey, R. D., & Wagenmakers, E. J. (2014). Simple relation between Bayesian order-restricted and point-null hypothesis tests. Statistics & Probability Letters, 92, 121-124.

• Morey, R. D., & Rouder, J. N. (2011). Bayes factor approaches for testing interval null hypotheses. Psychological methods, 16(4), 406.

• 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.

• 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

Examples

library(bayestestR)
prior <- distribution_normal(1000, mean = 0, sd = 1)
posterior <- distribution_normal(1000, mean = .5, sd = .3)
(BF_pars <- bayesfactor_parameters(posterior, prior, verbose = FALSE))

as.numeric(BF_pars)



# rstanarm models
# ---------------
contrasts(sleep$group) <- contr.equalprior_pairs # see vingette
stan_model <- suppressWarnings(stan_lmer(
  extra ~ group + (1 | ID),
  data = sleep,
  refresh = 0
))
bayesfactor_parameters(stan_model, verbose = FALSE)
bayesfactor_parameters(stan_model, null = rope_range(stan_model))

# emmGrid objects
# ---------------
group_diff <- pairs(emmeans(stan_model, ~group, data = sleep))
bayesfactor_parameters(group_diff, prior = stan_model, verbose = FALSE)

# Or
# group_diff_prior <- pairs(emmeans(unupdate(stan_model), ~group))
# bayesfactor_parameters(group_diff, prior = group_diff_prior, verbose = FALSE)



# brms models
# -----------
## Not run: 
contrasts(sleep$group) <- contr.equalprior_pairs # 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 <- suppressWarnings(brm(extra ~ group + (1 | ID),
  data = sleep,
  prior = my_custom_priors,
  refresh = 0
))
bayesfactor_parameters(brms_model, verbose = FALSE)

## End(Not run)

Description

This method computes Bayes factors for comparing a model with an order restrictions on its parameters with the fully unrestricted model. Note that this method should only be used for confirmatory analyses.

The ⁠bf_*⁠ function is an alias of the main function.

For more info, in particular on specifying correct priors for factors with more than 2 levels, see the Bayes factors vignette.

Usage

bayesfactor_restricted(posterior, ...)

bf_restricted(posterior, ...)

## S3 method for class 'stanreg'
bayesfactor_restricted(
  posterior,
  hypothesis,
  prior = NULL,
  verbose = TRUE,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  ...
)

## S3 method for class 'brmsfit'
bayesfactor_restricted(
  posterior,
  hypothesis,
  prior = NULL,
  verbose = TRUE,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  ...
)

## S3 method for class 'blavaan'
bayesfactor_restricted(
  posterior,
  hypothesis,
  prior = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'emmGrid'
bayesfactor_restricted(
  posterior,
  hypothesis,
  prior = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'data.frame'
bayesfactor_restricted(
  posterior,
  hypothesis,
  prior = NULL,
  rvar_col = NULL,
  ...
)

## S3 method for class 'bayesfactor_restricted'
as.logical(x, which = c("posterior", "prior"), ...)

Arguments

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 (Morey & Wagenmakers, 2013).

(Though it is possible to use bayesfactor_restricted() to test interval restrictions, it is more suitable for testing order restrictions; see examples).

Value

A data frame containing the (log) Bayes factor representing evidence against the un-restricted model (Use as.numeric() to extract the non-log Bayes factors; see examples). (A bool_results attribute contains the results for each sample, indicating if they are included or not in the hypothesized restriction.)

Setting the correct prior

For the computation of Bayes factors, the model priors must be proper priors (at the very least they should be not flat, and it is preferable that they be 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 the Jeffreys-Lindley-Bartlett paradox). Thus, you should only ever try (or want) to compute a Bayes factor when you have an informed prior.

(Note that by default, brms::brm() uses flat priors for fixed-effects; See example below.)

It is important to provide the correct prior for meaningful results, to match the posterior-type input:

• A numeric vector - prior should also be a numeric vector, representing the prior-estimate.

• A data frame - prior should also be a data frame, representing the prior-estimates, in matching column order.

  • If rvar_col is specified, prior should be the name of an rvar column that represents the prior-estimates.

• Supported Bayesian model (stanreg, brmsfit, etc.)

  • prior should be a model an equivalent model with MCMC samples from the priors only. See unupdate().

  • If prior is set to NULL, unupdate() is called internally (not supported for brmsfit_multiple model).

• Output from a {marginaleffects} function - prior should also be an equivalent output from a {marginaleffects} function based on a prior-model (See unupdate()).

• Output from an {emmeans} function

  • prior should also be an equivalent output from an {emmeans} function based on a prior-model (See unupdate()).

  • prior can also be the original (posterior) model, in which case the function will try to "unupdate" the estimates (not supported if the estimates have undergone any transformations – "log", "response", etc. – or any regriding).

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

References

• Morey, R. D., & Wagenmakers, E. J. (2014). Simple relation between Bayesian order-restricted and point-null hypothesis tests. Statistics & Probability Letters, 92, 121-124.

• Morey, R. D., & Rouder, J. N. (2011). Bayes factor approaches for testing interval null hypotheses. Psychological methods, 16(4), 406.

• Morey, R. D. (Jan, 2015). Multiple Comparisons with BayesFactor, Part 2 – order restrictions. Retrieved from https://richarddmorey.org/category/order-restrictions/.

Examples

set.seed(444)
library(bayestestR)
prior <- data.frame(
  A = rnorm(500),
  B = rnorm(500),
  C = rnorm(500)
)

posterior <- data.frame(
  A = rnorm(500, .4, 0.7),
  B = rnorm(500, -.2, 0.4),
  C = rnorm(500, 0, 0.5)
)

hyps <- c(
  "A > B & B > C",
  "A > B & A > C",
  "C > A"
)


(b <- bayesfactor_restricted(posterior, hypothesis = hyps, prior = prior))

bool <- as.logical(b, which = "posterior")
head(bool)



see::plots(
  plot(estimate_density(posterior)),
  # distribution **conditional** on the restrictions
  plot(estimate_density(posterior[bool[, hyps[1]], ])) + ggplot2::ggtitle(hyps[1]),
  plot(estimate_density(posterior[bool[, hyps[2]], ])) + ggplot2::ggtitle(hyps[2]),
  plot(estimate_density(posterior[bool[, hyps[3]], ])) + ggplot2::ggtitle(hyps[3]),
  guides = "collect"
)



# rstanarm models
# ---------------
data("mtcars")

fit_stan <- rstanarm::stan_glm(mpg ~ wt + cyl + am,
  data = mtcars, refresh = 0
)
hyps <- c(
  "am > 0 & cyl < 0",
  "cyl < 0",
  "wt - cyl > 0"
)

bayesfactor_restricted(fit_stan, hypothesis = hyps)




# emmGrid objects
# ---------------
# replicating http://bayesfactor.blogspot.com/2015/01/multiple-comparisons-with-bayesfactor-2.html
data("disgust")
contrasts(disgust$condition) <- contr.equalprior_pairs # see vignette
fit_model <- rstanarm::stan_glm(score ~ condition, data = disgust, family = gaussian())

em_condition <- emmeans::emmeans(fit_model, ~condition, data = disgust)
hyps <- c("lemon < control & control < sulfur")

bayesfactor_restricted(em_condition, prior = fit_model, hypothesis = hyps)
# > # Bayes Factor (Order-Restriction)
# >
# >                          Hypothesis P(Prior) P(Posterior)   BF
# >  lemon < control & control < sulfur     0.17         0.75 4.49
# > ---
# > Bayes factors for the restricted model vs. the un-restricted model.

Description

Compute the Bias Corrected and Accelerated Interval (BCa) of posterior distributions.

Usage

bci(x, ...)

bcai(x, ...)

## S3 method for class 'numeric'
bci(x, ci = 0.95, verbose = TRUE, ...)

## S3 method for class 'data.frame'
bci(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...)

## S3 method for class 'MCMCglmm'
bci(x, ci = 0.95, verbose = TRUE, ...)

## S3 method for class 'sim.merMod'
bci(
  x,
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'sim'
bci(x, ci = 0.95, parameters = NULL, verbose = TRUE, ...)

## S3 method for class 'emmGrid'
bci(x, ci = 0.95, verbose = TRUE, ...)

## S3 method for class 'slopes'
bci(x, ci = 0.95, verbose = TRUE, ...)

## S3 method for class 'stanreg'
bci(
  x,
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'brmsfit'
bci(
  x,
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'BFBayesFactor'
bci(x, ci = 0.95, verbose = TRUE, ...)

## S3 method for class 'get_predicted'
bci(x, ci = 0.95, use_iterations = FALSE, verbose = TRUE, ...)

Arguments

Details

Unlike equal-tailed intervals (see eti()) that typically exclude ⁠2.5%⁠ from each tail of the distribution and always include the median, the HDI is not equal-tailed and therefore always includes the mode(s) of posterior distributions. While this can be useful to better represent the credibility mass of a distribution, the HDI also has some limitations. See spi() for details.

The ⁠95%⁠ or ⁠89%⁠ Credible Intervals (CI) are two reasonable ranges to characterize the uncertainty related to the estimation (see here for a discussion about the differences between these two values).

The ⁠89%⁠ intervals (ci = 0.89) are deemed to be more stable than, for instance, ⁠95%⁠ intervals (Kruschke, 2014). An effective sample size of at least 10.000 is recommended if one wants to estimate ⁠95%⁠ intervals with high precision (Kruschke, 2014, p. 183ff). Unfortunately, the default number of posterior samples for most Bayes packages (e.g., rstanarm or brms) is only 4.000 (thus, you might want to increase it when fitting your model). Moreover, 89 indicates the arbitrariness of interval limits - its only remarkable property is being the highest prime number that does not exceed the already unstable ⁠95%⁠ threshold (McElreath, 2015).

However, ⁠95%⁠ has some advantages too. For instance, it shares (in the case of a normal posterior distribution) an intuitive relationship with the standard deviation and it conveys a more accurate image of the (artificial) bounds of the distribution. Also, because it is wider, it makes analyses more conservative (i.e., the probability of covering 0 is larger for the ⁠95%⁠ CI than for lower ranges such as ⁠89%⁠), which is a good thing in the context of the reproducibility crisis.

A ⁠95%⁠ equal-tailed interval (ETI) has ⁠2.5%⁠ of the distribution on either side of its limits. It indicates the 2.5th percentile and the 97.5h percentile. In symmetric distributions, the two methods of computing credible intervals, the ETI and the HDI, return similar results.

This is not the case for skewed distributions. Indeed, it is possible that parameter values in the ETI have lower credibility (are less probable) than parameter values outside the ETI. This property seems undesirable as a summary of the credible values in a distribution.

On the other hand, the ETI range does change when transformations are applied to the distribution (for instance, for a log odds scale to probabilities): the lower and higher bounds of the transformed distribution will correspond to the transformed lower and higher bounds of the original distribution. On the contrary, applying transformations to the distribution will change the resulting HDI.

Value

A data frame with following columns:

• Parameter The model parameter(s), if x is a model-object. If x is a vector, this column is missing.

• CI The probability of the credible interval.

• CI_low, CI_high The lower and upper credible interval limits for the parameters.

References

DiCiccio, T. J. and B. Efron. (1996). Bootstrap Confidence Intervals. Statistical Science. 11(3): 189–212. 10.1214/ss/1032280214

See Also

Other ci: ci(), cwi(), eti(), hdi(), si(), spi()

Examples

posterior <- rnorm(1000)
bci(posterior)
bci(posterior, ci = c(0.80, 0.89, 0.95))

Description

The difference between two Bayesian information criterion (BIC) indices of two models can be used to approximate Bayes factors via:

$BF_{10} = e^{(BIC_0 - BIC_1)/2}$

Usage

bic_to_bf(bic, denominator, log = FALSE)

Arguments

Value

The Bayes Factors corresponding to the BIC values against the denominator.

References

Wagenmakers, E. J. (2007). A practical solution to the pervasive problems of p values. Psychonomic bulletin & review, 14(5), 779-804

Examples

bic1 <- BIC(lm(Sepal.Length ~ 1, data = iris))
bic2 <- BIC(lm(Sepal.Length ~ Species, data = iris))
bic3 <- BIC(lm(Sepal.Length ~ Species + Petal.Length, data = iris))
bic4 <- BIC(lm(Sepal.Length ~ Species * Petal.Length, data = iris))

bic_to_bf(c(bic1, bic2, bic3, bic4), denominator = bic1)

Description

Performs a simple test to check whether the prior is informative to the posterior. This idea, and the accompanying heuristics, were discussed in Gelman et al. 2017.

Usage

check_prior(model, method = "gelman", simulate_priors = TRUE, ...)

Arguments

Value

A data frame with two columns: The parameter names and the quality of the prior (which might be "informative", "uninformative") or "not determinable" if the prior distribution could not be determined).

References

Gelman, A., Simpson, D., and Betancourt, M. (2017). The Prior Can Often Only Be Understood in the Context of the Likelihood. Entropy, 19(10), 555. doi:10.3390/e19100555

Examples

library(bayestestR)
model <- rstanarm::stan_glm(mpg ~ wt + am, data = mtcars, chains = 1, refresh = 0)
check_prior(model, method = "gelman")
check_prior(model, method = "lakeland")

# An extreme example where both methods diverge:
model <- rstanarm::stan_glm(mpg ~ wt,
  data = mtcars[1:3, ],
  prior = normal(-3.3, 1, FALSE),
  prior_intercept = normal(0, 1000, FALSE),
  refresh = 0
)
check_prior(model, method = "gelman")
check_prior(model, method = "lakeland")
# can provide visual confirmation to the Lakeland method
plot(si(model, verbose = FALSE))

Description

Compute Confidence/Credible/Compatibility Intervals (CI) or Support Intervals (SI) for Bayesian and frequentist models. The Documentation is accessible for:

Usage

ci(x, ...)

## S3 method for class 'numeric'
ci(x, ci = 0.95, method = "ETI", verbose = TRUE, BF = 1, ...)

## S3 method for class 'data.frame'
ci(x, ci = 0.95, method = "ETI", BF = 1, rvar_col = NULL, verbose = TRUE, ...)

## S3 method for class 'sim.merMod'
ci(
  x,
  ci = 0.95,
  method = "ETI",
  effects = c("fixed", "random", "all"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'sim'
ci(x, ci = 0.95, method = "ETI", parameters = NULL, verbose = TRUE, ...)

## S3 method for class 'stanreg'
ci(
  x,
  ci = 0.95,
  method = "ETI",
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  verbose = TRUE,
  BF = 1,
  ...
)

## S3 method for class 'brmsfit'
ci(
  x,
  ci = 0.95,
  method = "ETI",
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  parameters = NULL,
  verbose = TRUE,
  BF = 1,
  ...
)

## S3 method for class 'BFBayesFactor'
ci(x, ci = 0.95, method = "ETI", verbose = TRUE, BF = 1, ...)

## S3 method for class 'MCMCglmm'
ci(x, ci = 0.95, method = "ETI", verbose = TRUE, ...)

Arguments

Details

• Bayesian models

• Frequentist models

Value

A data frame with following columns:

• Parameter The model parameter(s), if x is a model-object. If x is a vector, this column is missing.

• CI The probability of the credible interval.

• CI_low, CI_high The lower and upper credible interval limits for the parameters.

Note

When it comes to interpretation, we recommend thinking of the CI in terms of an "uncertainty" or "compatibility" interval, the latter being defined as "Given any value in the interval and the background assumptions, the data should not seem very surprising" (Gelman & Greenland 2019).

There is also a plot()-method implemented in the see-package.

References

Gelman A, Greenland S. Are confidence intervals better termed "uncertainty intervals"? BMJ 2019;l5381. 10.1136/bmj.l5381

See Also

Other ci: bci(), cwi(), eti(), hdi(), si(), spi()

Examples

library(bayestestR)

posterior <- rnorm(1000)
ci(posterior, method = "ETI")
ci(posterior, method = "HDI")

df <- data.frame(replicate(4, rnorm(100)))
ci(df, method = "ETI", ci = c(0.80, 0.89, 0.95))
ci(df, method = "HDI", ci = c(0.80, 0.89, 0.95))

model <- suppressWarnings(
  stan_glm(mpg ~ wt, data = mtcars, chains = 2, iter = 200, refresh = 0)
)
ci(model, method = "ETI", ci = c(0.80, 0.89))
ci(model, method = "HDI", ci = c(0.80, 0.89))


bf <- ttestBF(x = rnorm(100, 1, 1))
ci(bf, method = "ETI")
ci(bf, method = "HDI")


model <- emtrends(model, ~1, "wt", data = mtcars)
ci(model, method = "ETI")
ci(model, method = "HDI")

Description

Build contrasts for factors with equal marginal priors on all levels. The 3 functions give the same orthogonal contrasts, but are scaled differently to allow different prior specifications (see 'Details'). Implementation from Singmann & Gronau's bfrms, following the description in Rouder, Morey, Speckman, & Province (2012, p. 363).

Usage

contr.equalprior(n, contrasts = TRUE, sparse = FALSE)

contr.equalprior_pairs(n, contrasts = TRUE, sparse = FALSE)

contr.equalprior_deviations(n, contrasts = TRUE, sparse = FALSE)

Arguments

Details

When using stats::contr.treatment, each dummy variable is the difference between each level and the reference level. While this is useful if setting different priors for each coefficient, it should not be used if one is trying to set a general prior for differences between means, as it (as well as stats::contr.sum and others) results in unequal marginal priors on the means the the difference between them.

library(brms)

data <- data.frame(
  group = factor(rep(LETTERS[1:4], each = 3)),
  y = rnorm(12)
)

contrasts(data$group) # R's default contr.treatment
#>   B C D
#> A 0 0 0
#> B 1 0 0
#> C 0 1 0
#> D 0 0 1

model_prior <- brm(
  y ~ group, data = data,
  sample_prior = "only",
  # Set the same priors on the 3 dummy variable
  # (Using an arbitrary scale)
  prior = set_prior("normal(0, 10)", coef = c("groupB", "groupC", "groupD"))
)

est <- emmeans::emmeans(model_prior, pairwise ~ group)

point_estimate(est, centr = "mean", disp = TRUE)
#> Point Estimate
#>
#> Parameter |  Mean |    SD
#> -------------------------
#> A         | -0.01 |  6.35
#> B         | -0.10 |  9.59
#> C         |  0.11 |  9.55
#> D         | -0.16 |  9.52
#> A - B     |  0.10 |  9.94
#> A - C     | -0.12 |  9.96
#> A - D     |  0.15 |  9.87
#> B - C     | -0.22 | 14.38
#> B - D     |  0.05 | 14.14
#> C - D     |  0.27 | 14.00

We can see that the priors for means aren't all the same (A having a more narrow prior), and likewise for the pairwise differences (priors for differences from A are more narrow).

The solution is to use one of the methods provided here, which do result in marginally equal priors on means differences between them. Though this will obscure the interpretation of parameters, setting equal priors on means and differences is important for they are useful for specifying equal priors on all means in a factor and their differences correct estimation of Bayes factors for contrasts and order restrictions of multi-level factors (where k>2). See info on specifying correct priors for factors with more than 2 levels in the Bayes factors vignette.

NOTE: When setting priors on these dummy variables, always:

1. Use priors that are centered on 0! Other location/centered priors are meaningless!

2. Use identically-scaled priors on all the dummy variables of a single factor!

contr.equalprior returns the original orthogonal-normal contrasts as described in Rouder, Morey, Speckman, & Province (2012, p. 363). Setting contrasts = FALSE returns the $I_{n} - \frac{1}{n}$ matrix.

contr.equalprior_pairs

Useful for setting priors in terms of pairwise differences between means - the scales of the priors defines the prior distribution of the pair-wise differences between all pairwise differences (e.g., A - B, B - C, etc.).

contrasts(data$group) <- contr.equalprior_pairs
contrasts(data$group)
#>         [,1]       [,2]       [,3]
#> A  0.0000000  0.6123724  0.0000000
#> B -0.1893048 -0.2041241  0.5454329
#> C -0.3777063 -0.2041241 -0.4366592
#> D  0.5670111 -0.2041241 -0.1087736

model_prior <- brm(
  y ~ group, data = data,
  sample_prior = "only",
  # Set the same priors on the 3 dummy variable
  # (Using an arbitrary scale)
  prior = set_prior("normal(0, 10)", coef = c("group1", "group2", "group3"))
)

est <- emmeans(model_prior, pairwise ~ group)

point_estimate(est, centr = "mean", disp = TRUE)
#> Point Estimate
#>
#> Parameter |  Mean |    SD
#> -------------------------
#> A         | -0.31 |  7.46
#> B         | -0.24 |  7.47
#> C         | -0.34 |  7.50
#> D         | -0.30 |  7.25
#> A - B     | -0.08 | 10.00
#> A - C     |  0.03 | 10.03
#> A - D     | -0.01 |  9.85
#> B - C     |  0.10 | 10.28
#> B - D     |  0.06 |  9.94
#> C - D     | -0.04 | 10.18

All means have the same prior distribution, and the distribution of the differences matches the prior we set of "normal(0, 10)". Success!

contr.equalprior_deviations

Useful for setting priors in terms of the deviations of each mean from the grand mean - the scales of the priors defines the prior distribution of the distance (above, below) the mean of one of the levels might have from the overall mean. (See examples.)

Value

A matrix with n rows and k columns, with k=n-1 if contrasts is TRUE and k=n if contrasts is FALSE.

References

Rouder, J. N., Morey, R. D., Speckman, P. L., & Province, J. M. (2012). Default Bayes factors for ANOVA designs. Journal of Mathematical Psychology, 56(5), 356-374. https://doi.org/10.1016/j.jmp.2012.08.001

Examples

contr.equalprior(2) # Q_2 in Rouder et al. (2012, p. 363)

contr.equalprior(5) # equivalent to Q_5 in Rouder et al. (2012, p. 363)

## check decomposition
Q3 <- contr.equalprior(3)
Q3 %*% t(Q3) ## 2/3 on diagonal and -1/3 on off-diagonal elements

Description

Refit Bayesian model as frequentist. Can be useful for comparisons.

Usage

convert_bayesian_as_frequentist(model, data = NULL, REML = TRUE)

bayesian_as_frequentist(model, data = NULL, REML = TRUE)

Arguments

Examples

# Rstanarm ----------------------
# Simple regressions
model <- rstanarm::stan_glm(Sepal.Length ~ Species,
  data = iris, chains = 2, refresh = 0
)
bayesian_as_frequentist(model)

model <- rstanarm::stan_glm(vs ~ mpg,
  family = "binomial",
  data = mtcars, chains = 2, refresh = 0
)
bayesian_as_frequentist(model)

# Mixed models
model <- rstanarm::stan_glmer(
  Sepal.Length ~ Petal.Length + (1 | Species),
  data = iris, chains = 2, refresh = 0
)
bayesian_as_frequentist(model)

model <- rstanarm::stan_glmer(vs ~ mpg + (1 | cyl),
  family = "binomial",
  data = mtcars, chains = 2, refresh = 0
)
bayesian_as_frequentist(model)

Description

Compute the Curvewise interval (CWI) (also called the "simultaneous interval" or "joint interval") of posterior distributions using ggdist::curve_interval(). Whereas the more typical "pointwise intervals" contain xx% of the posterior for a single parameter, joint/curvewise intervals contain xx% of the posterior distribution for all parameters.

Usage

cwi(x, ...)

## S3 method for class 'data.frame'
cwi(x, ci = 0.95, ...)

Arguments

Details

Applied model predictions, pointwise intervals contain xx% of the predicted response values conditional on specific predictor values. In contrast, curvewise intervals contain xx% of the predicted response values across all predictor values. Put another way, curvewise intervals contain xx% of the full prediction lines from the model.

For more details, see the ggdist documentation on curvewise intervals.

Value

A data frame with following columns:

• Parameter The model parameter(s), if x is a model-object. If x is a vector, this column is missing.

• CI The probability of the credible interval.

• CI_low, CI_high The lower and upper credible interval limits for the parameters.

See Also

Other ci: bci(), ci(), eti(), hdi(), si(), spi()

Examples

library(bayestestR)

if (require("ggplot2") && require("rstanarm") && require("ggdist")) {
  # Generate data =============================================
  k <- 11 # number of curves (iterations)
  n <- 201 # number of rows
  data <- data.frame(x = seq(-15, 15, length.out = n))

  # Simulate iterations as new columns
  for (i in 1:k) {
    data[paste0("iter_", i)] <- dnorm(data$x, seq(-5, 5, length.out = k)[i], 3)
  }

  # Note: first, we need to transpose the data to have iters as rows
  iters <- datawizard::data_transpose(data[paste0("iter_", 1:k)])

  # Compute Median
  data$Median <- point_estimate(iters)[["Median"]]

  # Compute Credible Intervals ================================

  # Compute ETI (default type of CI)
  data[c("ETI_low", "ETI_high")] <- eti(iters, ci = 0.5)[c("CI_low", "CI_high")]

  # Compute CWI
  # ggdist::curve_interval(reshape_iterations(data), iter_value .width = 0.5)

  # Visualization =============================================
  ggplot(data, aes(x = x, y = Median)) +
    geom_ribbon(aes(ymin = ETI_low, ymax = ETI_high), fill = "red", alpha = 0.3) +
    geom_line(linewidth = 1) +
    geom_line(
      data = reshape_iterations(data),
      aes(y = iter_value, group = iter_group),
      alpha = 0.3
    )
}

Description

Compute the density value at a given point of a distribution (i.e., the value of the y axis of a value x of a distribution).

Usage

density_at(posterior, x, precision = 2^10, method = "kernel", ...)

Arguments

Examples

library(bayestestR)
posterior <- distribution_normal(n = 10)
density_at(posterior, 0)
density_at(posterior, c(0, 1))

Description

Compute indices relevant to describe and characterize the posterior distributions.

Usage

describe_posterior(posterior, ...)

## S3 method for class 'numeric'
describe_posterior(
  posterior,
  centrality = "median",
  dispersion = FALSE,
  ci = 0.95,
  ci_method = "eti",
  test = c("p_direction", "rope"),
  rope_range = "default",
  rope_ci = 0.95,
  keep_iterations = FALSE,
  bf_prior = NULL,
  BF = 1,
  verbose = TRUE,
  ...
)

## S3 method for class 'data.frame'
describe_posterior(
  posterior,
  centrality = "median",
  dispersion = FALSE,
  ci = 0.95,
  ci_method = "eti",
  test = c("p_direction", "rope"),
  rope_range = "default",
  rope_ci = 0.95,
  keep_iterations = FALSE,
  bf_prior = NULL,
  BF = 1,
  rvar_col = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'stanreg'
describe_posterior(
  posterior,
  centrality = "median",
  dispersion = FALSE,
  ci = 0.95,
  ci_method = "eti",
  test = c("p_direction", "rope"),
  rope_range = "default",
  rope_ci = 0.95,
  keep_iterations = FALSE,
  bf_prior = NULL,
  diagnostic = c("ESS", "Rhat"),
  priors = FALSE,
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  BF = 1,
  verbose = TRUE,
  ...
)

## S3 method for class 'brmsfit'
describe_posterior(
  posterior,
  centrality = "median",
  dispersion = FALSE,
  ci = 0.95,
  ci_method = "eti",
  test = c("p_direction", "rope"),
  rope_range = "default",
  rope_ci = 0.95,
  keep_iterations = FALSE,
  bf_prior = NULL,
  diagnostic = c("ESS", "Rhat"),
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all", "location",
    "distributional", "auxiliary"),
  parameters = NULL,
  BF = 1,
  priors = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Details

One or more components of point estimates (like posterior mean or median), intervals and tests can be omitted from the summary output by setting the related argument to NULL. For example, test = NULL and centrality = NULL would only return the HDI (or CI).

References

• Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. (2019). Indices of Effect Existence and Significance in the Bayesian Framework. Frontiers in Psychology 2019;10:2767. doi:10.3389/fpsyg.2019.02767

• Region of Practical Equivalence (ROPE)

• Bayes factors

Examples

library(bayestestR)

if (require("logspline")) {
  x <- rnorm(1000)
  describe_posterior(x, verbose = FALSE)
  describe_posterior(x,
    centrality = "all",
    dispersion = TRUE,
    test = "all",
    verbose = FALSE
  )
  describe_posterior(x, ci = c(0.80, 0.90), verbose = FALSE)

  df <- data.frame(replicate(4, rnorm(100)))
  describe_posterior(df, verbose = FALSE)
  describe_posterior(
    df,
    centrality = "all",
    dispersion = TRUE,
    test = "all",
    verbose = FALSE
  )
  describe_posterior(df, ci = c(0.80, 0.90), verbose = FALSE)

  df <- data.frame(replicate(4, rnorm(20)))
  head(reshape_iterations(
    describe_posterior(df, keep_iterations = TRUE, verbose = FALSE)
  ))
}

# rstanarm models
# -----------------------------------------------
if (require("rstanarm") && require("emmeans")) {
  model <- suppressWarnings(
    stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
  )
  describe_posterior(model)
  describe_posterior(model, centrality = "all", dispersion = TRUE, test = "all")
  describe_posterior(model, ci = c(0.80, 0.90))
  describe_posterior(model, rope_range = list(c(-10, 5), c(-0.2, 0.2), "default"))

  # emmeans estimates
  # -----------------------------------------------
  describe_posterior(emtrends(model, ~1, "wt"))
}

# BayesFactor objects
# -----------------------------------------------
if (require("BayesFactor")) {
  bf <- ttestBF(x = rnorm(100, 1, 1))
  describe_posterior(bf)
  describe_posterior(bf, centrality = "all", dispersion = TRUE, test = "all")
  describe_posterior(bf, ci = c(0.80, 0.90))
}

Description

Returns a summary of the priors used in the model.

Usage

describe_prior(model, ...)

## S3 method for class 'brmsfit'
describe_prior(
  model,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all", "location",
    "distributional", "auxiliary"),
  parameters = NULL,
  ...
)

Arguments

Examples

library(bayestestR)

# rstanarm models
# -----------------------------------------------
if (require("rstanarm")) {
  model <- rstanarm::stan_glm(mpg ~ wt + cyl, data = mtcars)
  describe_prior(model)
}

# brms models
# -----------------------------------------------
if (require("brms")) {
  model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
  describe_prior(model)
}

# BayesFactor objects
# -----------------------------------------------
if (require("BayesFactor")) {
  bf <- ttestBF(x = rnorm(100, 1, 1))
  describe_prior(bf)
}

Description

Returns the accumulated log-posterior, the average Metropolis acceptance rate, divergent transitions, treedepth rather than terminated its evolution normally.

Usage

diagnostic_draws(posterior, ...)

Arguments

Examples

set.seed(333)

if (require("brms", quietly = TRUE)) {
  model <- suppressWarnings(brm(mpg ~ wt * cyl * vs,
    data = mtcars,
    iter = 100, control = list(adapt_delta = 0.80),
    refresh = 0
  ))
  diagnostic_draws(model)
}

Description

Extract diagnostic metrics (Effective Sample Size (ESS), Rhat and Monte Carlo Standard Error MCSE).

Usage

diagnostic_posterior(posterior, ...)

## Default S3 method:
diagnostic_posterior(posterior, diagnostic = c("ESS", "Rhat"), ...)

## S3 method for class 'stanreg'
diagnostic_posterior(
  posterior,
  diagnostic = "all",
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  ...
)

## S3 method for class 'brmsfit'
diagnostic_posterior(
  posterior,
  diagnostic = "all",
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  parameters = NULL,
  ...
)

Arguments

Details

Effective Sample (ESS) should be as large as possible, although for most applications, an effective sample size greater than 1000 is sufficient for stable estimates (Bürkner, 2017). The ESS corresponds to the number of independent samples with the same estimation power as the N autocorrelated samples. It is is a measure of "how much independent information there is in autocorrelated chains" (Kruschke 2015, p182-3).

Rhat should be the closest to 1. It should not be larger than 1.1 (Gelman and Rubin, 1992) or 1.01 (Vehtari et al., 2019). The split Rhat statistic quantifies the consistency of an ensemble of Markov chains.

Monte Carlo Standard Error (MCSE) is another measure of accuracy of the chains. It is defined as standard deviation of the chains divided by their effective sample size (the formula for mcse() is from Kruschke 2015, p. 187). The MCSE "provides a quantitative suggestion of how big the estimation noise is".

References

• Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical science, 7(4), 457-472.

• Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P. C. (2019). Rank-normalization, folding, and localization: An improved Rhat for assessing convergence of MCMC. arXiv preprint arXiv:1903.08008.

• Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press.

Examples

# rstanarm models
# -----------------------------------------------
model <- suppressWarnings(
  rstanarm::stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
)
diagnostic_posterior(model)

# brms models
# -----------------------------------------------
model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
diagnostic_posterior(model)

Description

A sample (simulated) dataset, used in tests and some examples.

Format

A data frame with 500 rows and 5 variables:

score

Score on the questionnaire, which ranges from 0 to 50 with higher scores representing harsher moral judgment

condition

one of three conditions, differing by the odor present in the room: a pleasant scent associated with cleanliness (lemon), a disgusting scent (sulfur), and a control condition in which no unusual odor is present

data("disgust")
head(disgust, n = 5)
#>   score condition
#> 1    13   control
#> 2    26   control
#> 3    30   control
#> 4    23   control
#> 5    34   control

Author(s)

Richard D. Morey

Description

Generate a sequence of n-quantiles, i.e., a sample of size n with a near-perfect distribution.

Usage

distribution(type = "normal", ...)

distribution_custom(n, type = "norm", ..., random = FALSE)

distribution_beta(n, shape1, shape2, ncp = 0, random = FALSE, ...)

distribution_binomial(n, size = 1, prob = 0.5, random = FALSE, ...)

distribution_binom(n, size = 1, prob = 0.5, random = FALSE, ...)

distribution_cauchy(n, location = 0, scale = 1, random = FALSE, ...)

distribution_chisquared(n, df, ncp = 0, random = FALSE, ...)

distribution_chisq(n, df, ncp = 0, random = FALSE, ...)

distribution_gamma(n, shape, scale = 1, random = FALSE, ...)

distribution_mixture_normal(n, mean = c(-3, 3), sd = 1, random = FALSE, ...)

distribution_normal(n, mean = 0, sd = 1, random = FALSE, ...)

distribution_gaussian(n, mean = 0, sd = 1, random = FALSE, ...)

distribution_nbinom(n, size, prob, mu, phi, random = FALSE, ...)

distribution_poisson(n, lambda = 1, random = FALSE, ...)

distribution_student(n, df, ncp, random = FALSE, ...)

distribution_t(n, df, ncp, random = FALSE, ...)

distribution_student_t(n, df, ncp, random = FALSE, ...)

distribution_tweedie(n, xi = NULL, mu, phi, power = NULL, random = FALSE, ...)

distribution_uniform(n, min = 0, max = 1, random = FALSE, ...)

rnorm_perfect(n, mean = 0, sd = 1)

Arguments

Details

When random = FALSE, these function return ⁠q*(ppoints(n), ...)⁠.

Examples

library(bayestestR)
x <- distribution(n = 10)
plot(density(x))

x <- distribution(type = "gamma", n = 100, shape = 2)
plot(density(x))

Description

This function returns the effective sample size (ESS).

Usage

effective_sample(model, ...)

## S3 method for class 'brmsfit'
effective_sample(
  model,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  parameters = NULL,
  ...
)

## S3 method for class 'stanreg'
effective_sample(
  model,
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  ...
)

Arguments

Details

Effective Sample (ESS) should be as large as possible, altough for most applications, an effective sample size greater than 1,000 is sufficient for stable estimates (Bürkner, 2017). The ESS corresponds to the number of independent samples with the same estimation power as the N autocorrelated samples. It is is a measure of “how much independent information there is in autocorrelated chains” (Kruschke 2015, p182-3).

Value

A data frame with two columns: Parameter name and effective sample size (ESS).

References

• Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press.

• Bürkner, P. C. (2017). brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1-28

Examples

library(rstanarm)
model <- suppressWarnings(
  stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
)
effective_sample(model)

Description

Perform a Test for Practical Equivalence for Bayesian and frequentist models.

Usage

equivalence_test(x, ...)

## Default S3 method:
equivalence_test(x, ...)

## S3 method for class 'data.frame'
equivalence_test(
  x,
  range = "default",
  ci = 0.95,
  rvar_col = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'stanreg'
equivalence_test(
  x,
  range = "default",
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'brmsfit'
equivalence_test(
  x,
  range = "default",
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

Arguments

Details

Documentation is accessible for:

• Bayesian models

• Frequentist models

For Bayesian models, the Test for Practical Equivalence is based on the "HDI+ROPE decision rule" (Kruschke, 2014, 2018) to check whether parameter values should be accepted or rejected against an explicitly formulated "null hypothesis" (i.e., a ROPE). In other words, it checks the percentage of the ⁠89%⁠ HDI that is the null region (the ROPE). If this percentage is sufficiently low, the null hypothesis is rejected. If this percentage is sufficiently high, the null hypothesis is accepted.

Using the ROPE and the HDI, Kruschke (2018) suggests using the percentage of the ⁠95%⁠ (or ⁠89%⁠, considered more stable) HDI that falls within the ROPE as a decision rule. If the HDI is completely outside the ROPE, the "null hypothesis" for this parameter is "rejected". If the ROPE completely covers the HDI, i.e., all most credible values of a parameter are inside the region of practical equivalence, the null hypothesis is accepted. Else, it’s undecided whether to accept or reject the null hypothesis. If the full ROPE is used (i.e., ⁠100%⁠ of the HDI), then the null hypothesis is rejected or accepted if the percentage of the posterior within the ROPE is smaller than to ⁠2.5%⁠ or greater than ⁠97.5%⁠. Desirable results are low proportions inside the ROPE (the closer to zero the better).

Some attention is required for finding suitable values for the ROPE limits (argument range). See 'Details' in rope_range() for further information.

Multicollinearity: Non-independent covariates

When parameters show strong correlations, i.e. when covariates are not independent, the joint parameter distributions may shift towards or away from the ROPE. In such cases, the test for practical equivalence may have inappropriate results. Collinearity invalidates ROPE and hypothesis testing based on univariate marginals, as the probabilities are conditional on independence. Most problematic are the results of the "undecided" parameters, which may either move further towards "rejection" or away from it (Kruschke 2014, 340f).

equivalence_test() performs a simple check for pairwise correlations between parameters, but as there can be collinearity between more than two variables, a first step to check the assumptions of this hypothesis testing is to look at different pair plots. An even more sophisticated check is the projection predictive variable selection (Piironen and Vehtari 2017).

Value

A data frame with following columns:

• Parameter The model parameter(s), if x is a model-object. If x is a vector, this column is missing.

• CI The probability of the HDI.

• ROPE_low, ROPE_high The limits of the ROPE. These values are identical for all parameters.

• ROPE_Percentage The proportion of the HDI that lies inside the ROPE.

• ROPE_Equivalence The "test result", as character. Either "rejected", "accepted" or "undecided".

• HDI_low , HDI_high The lower and upper HDI limits for the parameters.

Note

There is a print()-method with a digits-argument to control the amount of digits in the output, and there is a plot()-method to visualize the results from the equivalence-test (for models only).

References

• Kruschke, J. K. (2018). Rejecting or accepting parameter values in Bayesian estimation. Advances in Methods and Practices in Psychological Science, 1(2), 270-280. doi:10.1177/2515245918771304

• Kruschke, J. K. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press

• Piironen, J., & Vehtari, A. (2017). Comparison of Bayesian predictive methods for model selection. Statistics and Computing, 27(3), 711–735. doi:10.1007/s11222-016-9649-y

Examples

library(bayestestR)

equivalence_test(x = rnorm(1000, 0, 0.01), range = c(-0.1, 0.1))
equivalence_test(x = rnorm(1000, 0, 1), range = c(-0.1, 0.1))
equivalence_test(x = rnorm(1000, 1, 0.01), range = c(-0.1, 0.1))
equivalence_test(x = rnorm(1000, 1, 1), ci = c(.50, .99))

# print more digits
test <- equivalence_test(x = rnorm(1000, 1, 1), ci = c(.50, .99))
print(test, digits = 4)

model <- rstanarm::stan_glm(mpg ~ wt + cyl, data = mtcars)
equivalence_test(model)
# multiple ROPE ranges - asymmetric, symmetric, default
equivalence_test(model, range = list(c(10, 40), c(-5, -4), "default"))
# named ROPE ranges
equivalence_test(model, range = list(wt = c(-5, -4), `(Intercept)` = c(10, 40)))

# plot result
test <- equivalence_test(model)
plot(test)

equivalence_test(emmeans::emtrends(model, ~1, "wt", data = mtcars))

model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
equivalence_test(model)

bf <- BayesFactor::ttestBF(x = rnorm(100, 1, 1))
# equivalence_test(bf)

Description

This function is a wrapper over different methods of density estimation. By default, it uses the base R density with by default uses a different smoothing bandwidth ("SJ") from the legacy default implemented the base R density function ("nrd0"). However, Deng and Wickham suggest that method = "KernSmooth" is the fastest and the most accurate.

Usage

estimate_density(x, ...)

## S3 method for class 'data.frame'
estimate_density(
  x,
  method = "kernel",
  precision = 2^10,
  extend = FALSE,
  extend_scale = 0.1,
  bw = "SJ",
  ci = NULL,
  select = NULL,
  by = NULL,
  at = NULL,
  rvar_col = NULL,
  ...
)

Arguments

Note

There is also a plot()-method implemented in the see-package.

References

Deng, H., & Wickham, H. (2011). Density estimation in R. Electronic publication.

Examples

library(bayestestR)

set.seed(1)
x <- rnorm(250, mean = 1)

# Basic usage
density_kernel <- estimate_density(x) # default method is "kernel"

hist(x, prob = TRUE)
lines(density_kernel$x, density_kernel$y, col = "black", lwd = 2)
lines(density_kernel$x, density_kernel$CI_low, col = "gray", lty = 2)
lines(density_kernel$x, density_kernel$CI_high, col = "gray", lty = 2)
legend("topright",
  legend = c("Estimate", "95% CI"),
  col = c("black", "gray"), lwd = 2, lty = c(1, 2)
)

# Other Methods
density_logspline <- estimate_density(x, method = "logspline")
density_KernSmooth <- estimate_density(x, method = "KernSmooth")
density_mixture <- estimate_density(x, method = "mixture")

hist(x, prob = TRUE)
lines(density_kernel$x, density_kernel$y, col = "black", lwd = 2)
lines(density_logspline$x, density_logspline$y, col = "red", lwd = 2)
lines(density_KernSmooth$x, density_KernSmooth$y, col = "blue", lwd = 2)
lines(density_mixture$x, density_mixture$y, col = "green", lwd = 2)

# Extension
density_extended <- estimate_density(x, extend = TRUE)
density_default <- estimate_density(x, extend = FALSE)

hist(x, prob = TRUE)
lines(density_extended$x, density_extended$y, col = "red", lwd = 3)
lines(density_default$x, density_default$y, col = "black", lwd = 3)

# Multiple columns
head(estimate_density(iris))
head(estimate_density(iris, select = "Sepal.Width"))

# Grouped data
head(estimate_density(iris, by = "Species"))
head(estimate_density(iris$Petal.Width, by = iris$Species))

# rstanarm models
# -----------------------------------------------
library(rstanarm)
model <- suppressWarnings(
  stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
)
head(estimate_density(model))

library(emmeans)
head(estimate_density(emtrends(model, ~1, "wt", data = mtcars)))

# brms models
# -----------------------------------------------
library(brms)
model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
estimate_density(model)

Description

Compute the Equal-Tailed Interval (ETI) of posterior distributions using the quantiles method. The probability of being below this interval is equal to the probability of being above it. The ETI can be used in the context of uncertainty characterisation of posterior distributions as Credible Interval (CI).

Usage

eti(x, ...)

## S3 method for class 'numeric'
eti(x, ci = 0.95, verbose = TRUE, ...)

## S3 method for class 'data.frame'
eti(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...)

## S3 method for class 'stanreg'
eti(
  x,
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'brmsfit'
eti(
  x,
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'get_predicted'
eti(x, ci = 0.95, use_iterations = FALSE, verbose = TRUE, ...)

Arguments

Details

Unlike equal-tailed intervals (see eti()) that typically exclude ⁠2.5%⁠ from each tail of the distribution and always include the median, the HDI is not equal-tailed and therefore always includes the mode(s) of posterior distributions. While this can be useful to better represent the credibility mass of a distribution, the HDI also has some limitations. See spi() for details.

The ⁠95%⁠ or ⁠89%⁠ Credible Intervals (CI) are two reasonable ranges to characterize the uncertainty related to the estimation (see here for a discussion about the differences between these two values).

The ⁠89%⁠ intervals (ci = 0.89) are deemed to be more stable than, for instance, ⁠95%⁠ intervals (Kruschke, 2014). An effective sample size of at least 10.000 is recommended if one wants to estimate ⁠95%⁠ intervals with high precision (Kruschke, 2014, p. 183ff). Unfortunately, the default number of posterior samples for most Bayes packages (e.g., rstanarm or brms) is only 4.000 (thus, you might want to increase it when fitting your model). Moreover, 89 indicates the arbitrariness of interval limits - its only remarkable property is being the highest prime number that does not exceed the already unstable ⁠95%⁠ threshold (McElreath, 2015).

However, ⁠95%⁠ has some advantages too. For instance, it shares (in the case of a normal posterior distribution) an intuitive relationship with the standard deviation and it conveys a more accurate image of the (artificial) bounds of the distribution. Also, because it is wider, it makes analyses more conservative (i.e., the probability of covering 0 is larger for the ⁠95%⁠ CI than for lower ranges such as ⁠89%⁠), which is a good thing in the context of the reproducibility crisis.

A ⁠95%⁠ equal-tailed interval (ETI) has ⁠2.5%⁠ of the distribution on either side of its limits. It indicates the 2.5th percentile and the 97.5h percentile. In symmetric distributions, the two methods of computing credible intervals, the ETI and the HDI, return similar results.

This is not the case for skewed distributions. Indeed, it is possible that parameter values in the ETI have lower credibility (are less probable) than parameter values outside the ETI. This property seems undesirable as a summary of the credible values in a distribution.

On the other hand, the ETI range does change when transformations are applied to the distribution (for instance, for a log odds scale to probabilities): the lower and higher bounds of the transformed distribution will correspond to the transformed lower and higher bounds of the original distribution. On the contrary, applying transformations to the distribution will change the resulting HDI.

Value

A data frame with following columns:

• Parameter The model parameter(s), if x is a model-object. If x is a vector, this column is missing.

• CI The probability of the credible interval.

• CI_low, CI_high The lower and upper credible interval limits for the parameters.

See Also

Other ci: bci(), ci(), cwi(), hdi(), si(), spi()

Examples

library(bayestestR)

posterior <- rnorm(1000)
eti(posterior)
eti(posterior, ci = c(0.80, 0.89, 0.95))

df <- data.frame(replicate(4, rnorm(100)))
eti(df)
eti(df, ci = c(0.80, 0.89, 0.95))

model <- suppressWarnings(
  rstanarm::stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
)
eti(model)
eti(model, ci = c(0.80, 0.89, 0.95))

eti(emmeans::emtrends(model, ~1, "wt", data = mtcars))

model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
eti(model)
eti(model, ci = c(0.80, 0.89, 0.95))

bf <- BayesFactor::ttestBF(x = rnorm(100, 1, 1))
eti(bf)
eti(bf, ci = c(0.80, 0.89, 0.95))

Description

Compute the Highest Density Interval (HDI) of posterior distributions. All points within this interval have a higher probability density than points outside the interval. The HDI can be used in the context of uncertainty characterisation of posterior distributions as Credible Interval (CI).

Usage

hdi(x, ...)

## S3 method for class 'numeric'
hdi(x, ci = 0.95, verbose = TRUE, ...)

## S3 method for class 'data.frame'
hdi(x, ci = 0.95, rvar_col = NULL, verbose = TRUE, ...)

## S3 method for class 'stanreg'
hdi(
  x,
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'brmsfit'
hdi(
  x,
  ci = 0.95,
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  parameters = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'get_predicted'
hdi(x, ci = 0.95, use_iterations = FALSE, verbose = TRUE, ...)

Arguments

Details

Unlike equal-tailed intervals (see eti()) that typically exclude ⁠2.5%⁠ from each tail of the distribution and always include the median, the HDI is not equal-tailed and therefore always includes the mode(s) of posterior distributions. While this can be useful to better represent the credibility mass of a distribution, the HDI also has some limitations. See spi() for details.

The ⁠95%⁠ or ⁠89%⁠ Credible Intervals (CI) are two reasonable ranges to characterize the uncertainty related to the estimation (see here for a discussion about the differences between these two values).

The ⁠89%⁠ intervals (ci = 0.89) are deemed to be more stable than, for instance, ⁠95%⁠ intervals (Kruschke, 2014). An effective sample size of at least 10.000 is recommended if one wants to estimate ⁠95%⁠ intervals with high precision (Kruschke, 2014, p. 183ff). Unfortunately, the default number of posterior samples for most Bayes packages (e.g., rstanarm or brms) is only 4.000 (thus, you might want to increase it when fitting your model). Moreover, 89 indicates the arbitrariness of interval limits - its only remarkable property is being the highest prime number that does not exceed the already unstable ⁠95%⁠ threshold (McElreath, 2015).

However, ⁠95%⁠ has some advantages too. For instance, it shares (in the case of a normal posterior distribution) an intuitive relationship with the standard deviation and it conveys a more accurate image of the (artificial) bounds of the distribution. Also, because it is wider, it makes analyses more conservative (i.e., the probability of covering 0 is larger for the ⁠95%⁠ CI than for lower ranges such as ⁠89%⁠), which is a good thing in the context of the reproducibility crisis.

A ⁠95%⁠ equal-tailed interval (ETI) has ⁠2.5%⁠ of the distribution on either side of its limits. It indicates the 2.5th percentile and the 97.5h percentile. In symmetric distributions, the two methods of computing credible intervals, the ETI and the HDI, return similar results.

This is not the case for skewed distributions. Indeed, it is possible that parameter values in the ETI have lower credibility (are less probable) than parameter values outside the ETI. This property seems undesirable as a summary of the credible values in a distribution.

On the other hand, the ETI range does change when transformations are applied to the distribution (for instance, for a log odds scale to probabilities): the lower and higher bounds of the transformed distribution will correspond to the transformed lower and higher bounds of the original distribution. On the contrary, applying transformations to the distribution will change the resulting HDI.

Value

A data frame with following columns:

• Parameter The model parameter(s), if x is a model-object. If x is a vector, this column is missing.

• CI The probability of the credible interval.

• CI_low, CI_high The lower and upper credible interval limits for the parameters.

Note

There is also a plot()-method implemented in the see-package.

Author(s)

Credits go to ggdistribute and HDInterval.

References

• Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press.

• McElreath, R. (2015). Statistical rethinking: A Bayesian course with examples in R and Stan. Chapman and Hall/CRC.

See Also

Other interval functions, such as hdi(), eti(), bci(), spi(), si(), cwi().

Other ci: bci(), ci(), cwi(), eti(), si(), spi()

Examples

library(bayestestR)

posterior <- rnorm(1000)
hdi(posterior, ci = 0.89)
hdi(posterior, ci = c(0.80, 0.90, 0.95))

bayestestR::hdi(iris[1:4])
bayestestR::hdi(iris[1:4], ci = c(0.80, 0.90, 0.95))

model <- suppressWarnings(
  rstanarm::stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
)
bayestestR::hdi(model)
bayestestR::hdi(model, ci = c(0.80, 0.90, 0.95))

bayestestR::hdi(emmeans::emtrends(model, ~1, "wt", data = mtcars))

model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
bayestestR::hdi(model)
bayestestR::hdi(model, ci = c(0.80, 0.90, 0.95))

bf <- BayesFactor::ttestBF(x = rnorm(100, 1, 1))
bayestestR::hdi(bf)
bayestestR::hdi(bf, ci = c(0.80, 0.90, 0.95))

Description

Find the Highest Maximum A Posteriori probability estimate (MAP) of a posterior, i.e., the value associated with the highest probability density (the "peak" of the posterior distribution). In other words, it is an estimation of the mode for continuous parameters. Note that this function relies on estimate_density(), which by default uses a different smoothing bandwidth ("SJ") compared to the legacy default implemented the base R density() function ("nrd0").

Usage

map_estimate(x, ...)

## S3 method for class 'numeric'
map_estimate(x, precision = 2^10, method = "kernel", ...)

## S3 method for class 'stanreg'
map_estimate(
  x,
  precision = 2^10,
  method = "kernel",
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  ...
)

## S3 method for class 'brmsfit'
map_estimate(
  x,
  precision = 2^10,
  method = "kernel",
  effects = c("fixed", "random", "all"),
  component = c("conditional", "zi", "zero_inflated", "all"),
  parameters = NULL,
  ...
)

## S3 method for class 'data.frame'
map_estimate(x, precision = 2^10, method = "kernel", rvar_col = NULL, ...)

## S3 method for class 'get_predicted'
map_estimate(
  x,
  precision = 2^10,
  method = "kernel",
  use_iterations = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Value

A numeric value if x is a vector. If x is a model-object, returns a data frame with following columns:

• Parameter: The model parameter(s), if x is a model-object. If x is a vector, this column is missing.

• MAP_Estimate: The MAP estimate for the posterior or each model parameter.

Examples

library(bayestestR)

posterior <- rnorm(10000)
map_estimate(posterior)

plot(density(posterior))
abline(v = as.numeric(map_estimate(posterior)), col = "red")

model <- rstanarm::stan_glm(mpg ~ wt + cyl, data = mtcars)
map_estimate(model)

model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
map_estimate(model)

Description

This function returns the Monte Carlo Standard Error (MCSE).

Usage

mcse(model, ...)

## S3 method for class 'stanreg'
mcse(
  model,
  effects = c("fixed", "random", "all"),
  component = c("location", "all", "conditional", "smooth_terms", "sigma",
    "distributional", "auxiliary"),
  parameters = NULL,
  ...
)

Arguments

Details

Monte Carlo Standard Error (MCSE) is another measure of accuracy of the chains. It is defined as standard deviation of the chains divided by their effective sample size (the formula for mcse() is from Kruschke 2015, p. 187). The MCSE “provides a quantitative suggestion of how big the estimation noise is”.

References

Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press.

Examples

library(bayestestR)

model <- suppressWarnings(
  rstanarm::stan_glm(mpg ~ wt + am, data = mtcars, chains = 1, refresh = 0)
)
mcse(model)

Description

mediation() is a short summary for multivariate-response mediation-models, i.e. this function computes average direct and average causal mediation effects of multivariate response models.

Usage

mediation(model, ...)

## S3 method for class 'brmsfit'
mediation(
  model,
  treatment,
  mediator,
  response = NULL,
  centrality = "median",
  ci = 0.95,
  method = "ETI",
  ...
)

## S3 method for class 'stanmvreg'
mediation(
  model,
  treatment,
  mediator,
  response = NULL,
  centrality = "median",
  ci = 0.95,
  method = "ETI",
  ...
)

Arguments

Details

mediation() returns a data frame with information on the direct effect (mean value of posterior samples from treatment of the outcome model), mediator effect (mean value of posterior samples from mediator of the outcome model), indirect effect (mean value of the multiplication of the posterior samples from mediator of the outcome model and the posterior samples from treatment of the mediation model) and the total effect (mean value of sums of posterior samples used for the direct and indirect effect). The proportion mediated is the indirect effect divided by the total effect.