diff --git a/R/bayesfactor_models.R b/R/bayesfactor_models.R index 168cab5a3..742e414b5 100644 --- a/R/bayesfactor_models.R +++ b/R/bayesfactor_models.R @@ -9,13 +9,13 @@ #' @param ... Fitted models (see details), all fit on the same data, or a single #' `BFBayesFactor` object (see 'Details'). Ignored in `as.matrix()`, #' `update()`. If the following named arguments are present, they are passed -#' to [insight::get_loglikelihood] (see details): +#' to [`insight::get_loglikelihood()`] (see details): #' - `estimator` (defaults to `"ML"`) #' - `check_response` (defaults to `FALSE`) #' @param denominator Either an integer indicating which of the models to use as #' the denominator, or a model to be used as a denominator. Ignored for #' `BFBayesFactor`. -#' @param object,x A [bayesfactor_models()] object. +#' @param object,x A [`bayesfactor_models()`] object. #' @param subset Vector of model indices to keep or remove. #' @param reference Index of model to reference to, or `"top"` to #' reference to the best model, or `"bottom"` to reference to the worst @@ -26,9 +26,9 @@ #' implemented in the \href{https://easystats.github.io/see/}{\pkg{see}-package}. #' #' @details -#' If the passed models are supported by \pkg{insight} the DV of all models will be tested for equality -#' (else this is assumed to be true), and the models' terms will be extracted (allowing for follow-up -#' analysis with `bayesfactor_inclusion`). +#' 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 \CRANpkg{bridgesampling} package. #' - `brmsfit` models must have been fitted with `save_pars = save_pars(all = TRUE)`. @@ -44,9 +44,9 @@ #' 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 -#' (\cite{Gronau, Singmann, & Wagenmakers, 2017}). If less than 40,000 samples +#' (_Gronau, Singmann, & Wagenmakers, 2017_). If less than 40,000 samples #' are detected, `bayesfactor_models()` gives a warning. -#' \cr \cr +#' #' See also [the Bayes factors vignette](https://easystats.github.io/bayestestR/articles/bayes_factors.html). #' #' @inheritSection bayesfactor_parameters Interpreting Bayes Factors diff --git a/R/ci.R b/R/ci.R index e838dc215..cdc5d9eb9 100644 --- a/R/ci.R +++ b/R/ci.R @@ -4,36 +4,37 @@ #' (SI) for Bayesian and frequentist models. The Documentation is accessible #' for: #' -#' \itemize{ -#' \item [Bayesian models](https://easystats.github.io/bayestestR/articles/credible_interval.html) -#' \item [Frequentist models](https://easystats.github.io/parameters/reference/ci.default.html) -#' } +#' - [Bayesian models](https://easystats.github.io/bayestestR/articles/credible_interval.html) +#' - [Frequentist models](https://easystats.github.io/parameters/reference/ci.default.html) #' -#' @param x A `stanreg` or `brmsfit` model, or a vector representing a posterior distribution. -#' @param method Can be ['ETI'][eti] (default), ['HDI'][hdi], ['BCI'][bci], ['SPI'][spi] or ['SI'][si]. +#' @param x A `stanreg` or `brmsfit` model, or a vector representing a posterior +#' distribution. +#' @param method Can be ["ETI"][eti] (default), ["HDI"][hdi], ["BCI"][bci], +#' ["SPI"][spi] or ["SI"][si]. #' @param ci Value or vector of probability of the CI (between 0 and 1) -#' to be estimated. Default to `.95` (`95%`). +#' to be estimated. Default to `0.95` (`95%`). #' @inheritParams hdi #' @inheritParams si #' @inherit hdi seealso #' @family ci #' #' @return A data frame with following columns: -#' \itemize{ -#' \item `Parameter` The model parameter(s), if `x` is a model-object. If `x` is a vector, this column is missing. -#' \item `CI` The probability of the credible interval. -#' \item `CI_low`, `CI_high` The lower and upper credible interval limits for the parameters. -#' } +#' +#' - `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 -#' \dQuote{Given any value in the interval and the background assumptions, -#' the data should not seem very surprising} (\cite{Gelman & Greenland 2019}). -#' \cr \cr -#' There is also a [`plot()`-method](https://easystats.github.io/see/articles/bayestestR.html) implemented in the \href{https://easystats.github.io/see/}{\pkg{see}-package}. +#' "Given any value in the interval and the background assumptions, +#' the data should not seem very surprising" (_Gelman & Greenland 2019_). #' -#' @references Gelman A, Greenland S. Are confidence intervals better termed "uncertainty intervals"? BMJ 2019;l5381. 10.1136/bmj.l5381 +#' There is also a [`plot()`-method](https://easystats.github.io/see/articles/bayestestR.html) implemented in the \href{https://easystats.github.io/see/}{\pkg{see}-package}. #' +#' @references +#' Gelman A, Greenland S. Are confidence intervals better termed "uncertainty +#' intervals"? BMJ 2019;l5381. 10.1136/bmj.l5381 #' #' @examplesIf require("rstanarm", quietly = TRUE) #' library(bayestestR) @@ -138,9 +139,12 @@ ci <- function(x, ...) { ) ) } else { - insight::format_error( - "`method` should be 'ETI' (for equal-tailed interval),'HDI' (for highest density interval), 'BCI' (for bias corrected and accelerated bootstrap intervals), 'SPI' (for shortest probability interval) or 'SI' (for support interval)." - ) + insight::format_error(paste0( + "`method` should be 'ETI' (for equal-tailed interval), ", + "'HDI' (for highest density interval), 'BCI' (for bias corrected and ", + "accelerated bootstrap intervals), 'SPI' (for shortest probability ", + "interval) or 'SI' (for support interval)." + )) } } @@ -157,7 +161,6 @@ ci.numeric <- function(x, ci = 0.95, method = "ETI", verbose = TRUE, BF = 1, ... ci.data.frame <- ci.numeric - #' @export ci.draws <- function(x, ci = 0.95, method = "ETI", verbose = TRUE, BF = 1, ...) { .ci_bayesian(.posterior_draws_to_df(x), ci = ci, method = method, verbose = verbose, BF = BF, ...) @@ -185,7 +188,6 @@ ci.emmGrid <- function(x, ci = NULL, ...) { ci.emm_list <- ci.emmGrid - #' @rdname ci #' @export ci.sim.merMod <- function(x, @@ -207,7 +209,6 @@ ci.sim.merMod <- function(x, } - #' @rdname ci #' @export ci.sim <- function(x, @@ -227,7 +228,6 @@ ci.sim <- function(x, } - #' @rdname ci #' @export ci.stanreg <- function(x, @@ -285,22 +285,17 @@ ci.brmsfit <- function(x, ) } - #' @export ci.stanfit <- ci.stanreg - #' @export ci.blavaan <- ci.stanreg - #' @rdname ci #' @export ci.BFBayesFactor <- ci.numeric - - #' @rdname ci #' @export ci.MCMCglmm <- function(x, ci = 0.95, method = "ETI", verbose = TRUE, ...) { @@ -338,22 +333,17 @@ ci.bcplm <- function(x, ci = 0.95, method = "ETI", verbose = TRUE, ...) { ci(insight::get_parameters(x), ci = ci, method = method, verbose = verbose, ...) } - #' @export ci.blrm <- ci.bcplm - #' @export ci.mcmc <- ci.bcplm - #' @export ci.mcmc.list <- ci.bcplm - #' @export ci.BGGM <- ci.bcplm - #' @export ci.get_predicted <- ci.data.frame diff --git a/R/describe_posterior.R b/R/describe_posterior.R index e62b055ad..b292cfaee 100644 --- a/R/describe_posterior.R +++ b/R/describe_posterior.R @@ -92,7 +92,7 @@ #' #' # emmeans estimates #' # ----------------------------------------------- -#' describe_posterior(suppressWarnings(emtrends(model, ~1, "wt"))) +#' describe_posterior(emtrends(model, ~1, "wt")) #' } #' #' # BayesFactor objects diff --git a/R/equivalence_test.R b/R/equivalence_test.R index 1703e8f88..a0da152f6 100644 --- a/R/equivalence_test.R +++ b/R/equivalence_test.R @@ -89,9 +89,7 @@ #' plot(test) #' #' library(emmeans) -#' equivalence_test( -#' suppressWarnings(emtrends(model, ~1, "wt", data = mtcars)) -#' ) +#' equivalence_test(emtrends(model, ~1, "wt", data = mtcars)) #' #' library(brms) #' model <- brms::brm(mpg ~ wt + cyl, data = mtcars) diff --git a/R/estimate_density.R b/R/estimate_density.R index f65184749..554290adf 100644 --- a/R/estimate_density.R +++ b/R/estimate_density.R @@ -82,9 +82,7 @@ #' head(estimate_density(model)) #' #' library(emmeans) -#' head(estimate_density(suppressWarnings( -#' emtrends(model, ~1, "wt", data = mtcars) -#' ))) +#' head(estimate_density(emtrends(model, ~1, "wt", data = mtcars))) #' #' # brms models #' # ----------------------------------------------- diff --git a/R/eti.R b/R/eti.R index 9ec8cb12e..04d47498d 100644 --- a/R/eti.R +++ b/R/eti.R @@ -31,7 +31,7 @@ #' eti(model, ci = c(0.80, 0.89, 0.95)) #' #' library(emmeans) -#' eti(suppressWarnings(emtrends(model, ~1, "wt", data = mtcars))) +#' eti(emtrends(model, ~1, "wt", data = mtcars)) #' #' library(brms) #' model <- brms::brm(mpg ~ wt + cyl, data = mtcars) diff --git a/R/hdi.R b/R/hdi.R index 9acea5f1c..0b365de90 100644 --- a/R/hdi.R +++ b/R/hdi.R @@ -31,22 +31,24 @@ #' 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 +#' mass of a distribution, the HDI also has some limitations. See [`spi()`] for #' details. -#' \cr \cr +#' #' The [`95%` or `89%` Credible Intervals (CI)](https://easystats.github.io/bayestestR/articles/credible_interval.html) -#' are two reasonable ranges to characterize the uncertainty related to the estimation (see [here](https://easystats.github.io/bayestestR/articles/credible_interval.html) for a discussion about the differences between these two values). -#' \cr +#' are two reasonable ranges to characterize the uncertainty related to the +#' estimation (see [here](https://easystats.github.io/bayestestR/articles/credible_interval.html) +#' 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 (\cite{Kruschke, 2014}). An effective sample size +#' 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 (\cite{Kruschke, 2014, p. 183ff}). Unfortunately, the +#' 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 (\cite{McElreath, 2015}). -#' \cr +#' exceed the already unstable `95%` threshold (_McElreath, 2015_). +#' #' However, `95%` has some [advantages #' too](https://easystats.github.io/blog/posts/bayestestr_95/). For instance, it #' shares (in the case of a normal posterior distribution) an intuitive @@ -55,17 +57,17 @@ #' 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. -#' \cr \cr +#' #' 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][hdi], return similar results. -#' \cr \cr +#' #' 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. -#' \cr \cr +#' #' 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 @@ -96,9 +98,7 @@ #' bayestestR::hdi(model, ci = c(0.80, 0.90, 0.95)) #' #' library(emmeans) -#' bayestestR::hdi(suppressWarnings( -#' emtrends(model, ~1, "wt", data = mtcars) -#' )) +#' bayestestR::hdi(emtrends(model, ~1, "wt", data = mtcars)) #' #' library(brms) #' model <- brms::brm(mpg ~ wt + cyl, data = mtcars) diff --git a/R/p_direction.R b/R/p_direction.R index 938386534..b743ce9d2 100644 --- a/R/p_direction.R +++ b/R/p_direction.R @@ -106,9 +106,7 @@ #' # emmeans #' # ----------------------------------------------- #' if (require("emmeans")) { -#' p_direction(suppressWarnings( -#' emtrends(model, ~1, "wt", data = mtcars) -#' )) +#' p_direction(emtrends(model, ~1, "wt", data = mtcars)) #' } #' #' # brms models diff --git a/R/point_estimate.R b/R/point_estimate.R index e139f4014..a461c3cba 100644 --- a/R/point_estimate.R +++ b/R/point_estimate.R @@ -44,7 +44,7 @@ #' # ----------------------------------------------- #' library(emmeans) #' point_estimate( -#' suppressWarnings(emtrends(model, ~1, "wt", data = mtcars)), +#' emtrends(model, ~1, "wt", data = mtcars), #' centrality = c("median", "MAP") #' ) #' diff --git a/R/rope.R b/R/rope.R index 925caba3a..cbc79cefc 100644 --- a/R/rope.R +++ b/R/rope.R @@ -1,87 +1,99 @@ #' Region of Practical Equivalence (ROPE) #' -#' Compute the proportion of the HDI (default to the `89%` HDI) of a posterior distribution that lies within a region of practical equivalence. +#' Compute the proportion of the HDI (default to the `89%` HDI) of a posterior +#' distribution that lies within a region of practical equivalence. #' -#' @param x Vector representing a posterior distribution. Can also be a `stanreg` or `brmsfit` model. +#' @param x Vector representing a posterior distribution. Can also be a +#' `stanreg` or `brmsfit` model. #' @param range ROPE's lower and higher bounds. Should be `"default"` or -#' depending on the number of outcome variables a vector or a list. In models with one response, -#' `range` should be a vector of length two (e.g., `c(-0.1, 0.1)`). In -#' multivariate models, `range` should be a list with a numeric vectors for -#' each response variable. Vector names should correspond to the name of the response -#' variables. If `"default"` and input is a vector, the range is set to `c(-0.1, -#' 0.1)`. If `"default"` and input is a Bayesian model, -#' [`rope_range()`][rope_range] is used. +#' depending on the number of outcome variables a vector or a list. In +#' models with one response, `range` should be a vector of length two (e.g., +#' `c(-0.1, 0.1)`). In multivariate models, `range` should be a list with a +#' numeric vectors for each response variable. Vector names should correspond +#' to the name of the response variables. If `"default"` and input is a vector, +#' the range is set to `c(-0.1, 0.1)`. If `"default"` and input is a Bayesian +#' model, [`rope_range()`][rope_range] is used. #' @param ci The Credible Interval (CI) probability, corresponding to the -#' proportion of HDI, to use for the percentage in ROPE. +#' proportion of HDI, to use for the percentage in ROPE. #' @param ci_method The type of interval to use to quantify the percentage in -#' ROPE. Can be 'HDI' (default) or 'ETI'. See [ci()]. +#' ROPE. Can be 'HDI' (default) or 'ETI'. See [ci()]. #' #' @inheritParams hdi #' -#' @details -#' \subsection{ROPE}{ -#' Statistically, the probability of a posterior distribution of being -#' different from 0 does not make much sense (the probability of a single value -#' null hypothesis in a continuous distribution is 0). Therefore, the idea -#' underlining ROPE is to let the user define an area around the null value -#' enclosing values that are *equivalent to the null* value for practical -#' purposes (\cite{Kruschke 2010, 2011, 2014}). -#' \cr \cr -#' Kruschke (2018) suggests that such null value could be set, by default, -#' to the -0.1 to 0.1 range of a standardized parameter (negligible effect -#' size according to Cohen, 1988). This could be generalized: For instance, -#' for linear models, the ROPE could be set as `0 +/- .1 * sd(y)`. -#' This ROPE range can be automatically computed for models using the -#' [rope_range] function. -#' \cr \cr -#' Kruschke (2010, 2011, 2014) suggests using the proportion of the `95%` -#' (or `89%`, considered more stable) [HDI][hdi] that falls within the -#' ROPE as an index for "null-hypothesis" testing (as understood under the -#' Bayesian framework, see [`equivalence_test()`][equivalence_test]). -#' } -#' \subsection{Sensitivity to parameter's scale}{ -#' It is important to consider the unit (i.e., the scale) of the predictors -#' when using an index based on the ROPE, as the correct interpretation of the -#' ROPE as representing a region of practical equivalence to zero is dependent -#' on the scale of the predictors. Indeed, the percentage in ROPE depend on -#' the unit of its parameter. In other words, as the ROPE represents a fixed -#' portion of the response's scale, its proximity with a coefficient depends -#' on the scale of the coefficient itself. -#' } -#' \subsection{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. Collinearity invalidates ROPE and hypothesis -#' testing based on univariate marginals, as the probabilities are conditional -#' on independence. Most problematic are parameters that only have partial -#' overlap with the ROPE region. In case of collinearity, the (joint) distributions -#' of these parameters may either get an increased or decreased ROPE, which -#' means that inferences based on `rope()` are inappropriate -#' (\cite{Kruschke 2014, 340f}). -#' \cr \cr -#' `rope()` 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 (\cite{Piironen and Vehtari 2017}). -#' } -#' \subsection{Strengths and Limitations}{ -#' **Strengths:** Provides information related to the practical relevance of the effects. -#' \cr \cr -#' **Limitations:** A ROPE range needs to be arbitrarily defined. Sensitive to the scale (the unit) of the predictors. Not sensitive to highly significant effects. -#' } +#' @section ROPE: +#' Statistically, the probability of a posterior distribution of being +#' different from 0 does not make much sense (the probability of a single value +#' null hypothesis in a continuous distribution is 0). Therefore, the idea +#' underlining ROPE is to let the user define an area around the null value +#' enclosing values that are *equivalent to the null* value for practical +#' purposes (_Kruschke 2010, 2011, 2014_). +#' +#' Kruschke (2018) suggests that such null value could be set, by default, +#' to the -0.1 to 0.1 range of a standardized parameter (negligible effect +#' size according to Cohen, 1988). This could be generalized: For instance, +#' for linear models, the ROPE could be set as `0 +/- .1 * sd(y)`. +#' This ROPE range can be automatically computed for models using the +#' [rope_range] function. +#' +#' Kruschke (2010, 2011, 2014) suggests using the proportion of the `95%` +#' (or `89%`, considered more stable) [HDI][hdi] that falls within the +#' ROPE as an index for "null-hypothesis" testing (as understood under the +#' Bayesian framework, see [`equivalence_test()`][equivalence_test]). +#' +#' @section Sensitivity to parameter's scale: +#' It is important to consider the unit (i.e., the scale) of the predictors +#' when using an index based on the ROPE, as the correct interpretation of the +#' ROPE as representing a region of practical equivalence to zero is dependent +#' on the scale of the predictors. Indeed, the percentage in ROPE depend on +#' the unit of its parameter. In other words, as the ROPE represents a fixed +#' portion of the response's scale, its proximity with a coefficient depends +#' on the scale of the coefficient itself. +#' +#' @section 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. Collinearity invalidates ROPE and hypothesis +#' testing based on univariate marginals, as the probabilities are conditional +#' on independence. Most problematic are parameters that only have partial +#' overlap with the ROPE region. In case of collinearity, the (joint) distributions +#' of these parameters may either get an increased or decreased ROPE, which +#' means that inferences based on `rope()` are inappropriate +#' (_Kruschke 2014, 340f_). +#' +#' `rope()` 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_). +#' +#' @section Strengths and Limitations: +#' **Strengths:** Provides information related to the practical relevance of +#' the effects. +#' +#' **Limitations:** A ROPE range needs to be arbitrarily defined. Sensitive to +#' the scale (the unit) of the predictors. Not sensitive to highly significant +#' effects. #' #' @note There is also a [`plot()`-method](https://easystats.github.io/see/articles/bayestestR.html) implemented in the \href{https://easystats.github.io/see/}{\pkg{see}-package}. #' -#' @references \itemize{ -#' \item Cohen, J. (1988). Statistical power analysis for the behavioural sciences. -#' \item Kruschke, J. K. (2010). What to believe: Bayesian methods for data analysis. Trends in cognitive sciences, 14(7), 293-300. \doi{10.1016/j.tics.2010.05.001}. -#' \item Kruschke, J. K. (2011). Bayesian assessment of null values via parameter estimation and model comparison. Perspectives on Psychological Science, 6(3), 299-312. \doi{10.1177/1745691611406925}. -#' \item Kruschke, J. K. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press. \doi{10.1177/2515245918771304}. -#' \item 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}. -#' \item Makowski D, Ben-Shachar MS, Chen SHA, 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} -#' \item 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} -#' } +#' @references +#' - Cohen, J. (1988). Statistical power analysis for the behavioural sciences. +#' - Kruschke, J. K. (2010). What to believe: Bayesian methods for data analysis. +#' Trends in cognitive sciences, 14(7), 293-300. \doi{10.1016/j.tics.2010.05.001}. +#' - Kruschke, J. K. (2011). Bayesian assessment of null values via parameter +#' estimation and model comparison. Perspectives on Psychological Science, +#' 6(3), 299-312. \doi{10.1177/1745691611406925}. +#' - Kruschke, J. K. (2014). Doing Bayesian data analysis: A tutorial with R, +#' JAGS, and Stan. Academic Press. \doi{10.1177/2515245918771304}. +#' - 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}. +#' - Makowski D, Ben-Shachar MS, Chen SHA, 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} +#' - 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) @@ -89,32 +101,36 @@ #' rope(x = rnorm(1000, 0, 0.01), range = c(-0.1, 0.1)) #' rope(x = rnorm(1000, 0, 1), range = c(-0.1, 0.1)) #' rope(x = rnorm(1000, 1, 0.01), range = c(-0.1, 0.1)) -#' rope(x = rnorm(1000, 1, 1), ci = c(.90, .95)) +#' rope(x = rnorm(1000, 1, 1), ci = c(0.90, 0.95)) #' \dontrun{ #' library(rstanarm) #' model <- suppressWarnings( #' stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0) #' ) #' rope(model) -#' rope(model, ci = c(.90, .95)) +#' rope(model, ci = c(0.90, 0.95)) #' #' library(emmeans) -#' rope(emtrends(model, ~1, "wt"), ci = c(.90, .95)) +#' rope(emtrends(model, ~1, "wt"), ci = c(0.90, 0.95)) #' #' library(brms) -#' model <- brms::brm(mpg ~ wt + cyl, data = mtcars) +#' model <- brm(mpg ~ wt + cyl, data = mtcars) #' rope(model) -#' rope(model, ci = c(.90, .95)) +#' rope(model, ci = c(0.90, 0.95)) #' #' library(brms) -#' model <- brms::brm(brms::mvbind(mpg, disp) ~ wt + cyl, data = mtcars) +#' model <- brm( +#' bf(mvbind(mpg, disp) ~ wt + cyl) + set_rescor(rescor = TRUE), +#' data = mtcars, +#' set_rescore +#' ) #' rope(model) -#' rope(model, ci = c(.90, .95)) +#' rope(model, ci = c(0.90, 0.95)) #' #' library(BayesFactor) #' bf <- ttestBF(x = rnorm(100, 1, 1)) #' rope(bf) -#' rope(bf, ci = c(.90, .95)) +#' rope(bf, ci = c(0.90, 0.95)) #' } #' @export rope <- function(x, ...) { @@ -151,7 +167,6 @@ rope.numeric <- function(x, range = "default", ci = 0.95, ci_method = "ETI", ver # "do.call(rbind)" does not bind attribute values together # so we need to capture the information about HDI separately - out <- do.call(rbind, rope_values) if (nrow(out) > 1) { out$ROPE_Percentage <- as.numeric(out$ROPE_Percentage) @@ -190,7 +205,6 @@ rope.data.frame <- function(x, range = "default", ci = 0.95, ci_method = "ETI", } - #' @export rope.draws <- function(x, range = "default", ci = 0.95, ci_method = "ETI", verbose = TRUE, ...) { rope(.posterior_draws_to_df(x), range = range, ci = ci, ci_method = ci_method, verbose = verbose, ...) @@ -200,7 +214,6 @@ rope.draws <- function(x, range = "default", ci = 0.95, ci_method = "ETI", verbo rope.rvar <- rope.draws - #' @export rope.emmGrid <- function(x, range = "default", ci = 0.95, ci_method = "ETI", verbose = TRUE, ...) { xdf <- insight::get_parameters(x) @@ -214,7 +227,6 @@ rope.emmGrid <- function(x, range = "default", ci = 0.95, ci_method = "ETI", ver rope.emm_list <- rope.emmGrid - #' @export rope.BFBayesFactor <- function(x, range = "default", ci = 0.95, ci_method = "ETI", verbose = TRUE, ...) { if (all(range == "default")) { @@ -225,7 +237,6 @@ rope.BFBayesFactor <- function(x, range = "default", ci = 0.95, ci_method = "ETI out } - #' @export rope.bamlss <- rope.BFBayesFactor @@ -276,7 +287,6 @@ rope.BGGM <- rope.bcplm rope.mcmc.list <- rope.bcplm - #' @keywords internal .rope <- function(x, range = c(-0.1, 0.1), ci = 0.95, ci_method = "ETI", verbose = TRUE) { ci_bounds <- ci(x, ci = ci, method = ci_method, verbose = verbose) @@ -304,7 +314,6 @@ rope.mcmc.list <- rope.bcplm } - #' @rdname rope #' @export rope.stanreg <- function(x, range = "default", ci = 0.95, ci_method = "ETI", effects = c("fixed", "random", "all"), component = c("location", "all", "conditional", "smooth_terms", "sigma", "distributional", "auxiliary"), parameters = NULL, verbose = TRUE, ...) { @@ -433,7 +442,6 @@ rope.brmsfit <- function(x, } - #' @export rope.sim.merMod <- function(x, range = "default", @@ -503,7 +511,6 @@ rope.sim.merMod <- function(x, } - #' @export rope.sim <- function(x, range = "default", ci = 0.95, ci_method = "ETI", parameters = NULL, verbose = TRUE, ...) { if (all(range == "default")) { @@ -539,8 +546,6 @@ rope.sim <- function(x, range = "default", ci = 0.95, ci_method = "ETI", paramet } - - #' @keywords internal .prepare_rope_df <- function(parms, range, ci, ci_method, verbose) { tmp <- sapply( diff --git a/R/si.R b/R/si.R index 6d1e41e22..ad67b787f 100644 --- a/R/si.R +++ b/R/si.R @@ -79,7 +79,8 @@ #' si(brms_model) #' } #' @references -#' Wagenmakers, E., Gronau, Q. F., Dablander, F., & Etz, A. (2018, November 22). The Support Interval. \doi{10.31234/osf.io/zwnxb} +#' Wagenmakers, E., Gronau, Q. F., Dablander, F., & Etz, A. (2018, November 22). +#' The Support Interval. \doi{10.31234/osf.io/zwnxb} #' #' @export si <- function(posterior, prior = NULL, BF = 1, verbose = TRUE, ...) { diff --git a/man/bayesfactor_models.Rd b/man/bayesfactor_models.Rd index 2a7b69060..89e83db95 100644 --- a/man/bayesfactor_models.Rd +++ b/man/bayesfactor_models.Rd @@ -22,7 +22,7 @@ bf_models(..., denominator = 1, verbose = TRUE) \item{...}{Fitted models (see details), all fit on the same data, or a single \code{BFBayesFactor} object (see 'Details'). Ignored in \code{as.matrix()}, \code{update()}. If the following named arguments are present, they are passed -to \link[insight:get_loglikelihood]{insight::get_loglikelihood} (see details): +to \code{\link[insight:get_loglikelihood]{insight::get_loglikelihood()}} (see details): \itemize{ \item \code{estimator} (defaults to \code{"ML"}) \item \code{check_response} (defaults to \code{FALSE}) @@ -53,9 +53,9 @@ This function computes or extracts Bayes factors from fitted models. The \verb{bf_*} function is an alias of the main function. } \details{ -If the passed models are supported by \pkg{insight} the DV of all models will be tested for equality -(else this is assumed to be true), and the models' terms will be extracted (allowing for follow-up -analysis with \code{bayesfactor_inclusion}). +If the passed models are supported by \strong{insight} the DV of all models will +be tested for equality (else this is assumed to be true), and the models' +terms will be extracted (allowing for follow-up analysis with \code{bayesfactor_inclusion}). \itemize{ \item For \code{brmsfit} or \code{stanreg} models, Bayes factors are computed using the \CRANpkg{bridgesampling} package. \itemize{ @@ -74,9 +74,9 @@ are the 4 P's: \strong{P}roper \strong{P}riors and \strong{P}lentiful 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 -(\cite{Gronau, Singmann, & Wagenmakers, 2017}). If less than 40,000 samples +(\emph{Gronau, Singmann, & Wagenmakers, 2017}). If less than 40,000 samples are detected, \code{bayesfactor_models()} gives a warning. -\cr \cr + See also \href{https://easystats.github.io/bayestestR/articles/bayes_factors.html}{the Bayes factors vignette}. } \note{ diff --git a/man/bci.Rd b/man/bci.Rd index 2bace9b0f..48ac91f10 100644 --- a/man/bci.Rd +++ b/man/bci.Rd @@ -90,7 +90,8 @@ be abbreviated. Only applies to \pkg{brms}-models.} \value{ A data frame with following columns: \itemize{ -\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a vector, this column is missing. +\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a +vector, this column is missing. \item \code{CI} The probability of the credible interval. \item \code{CI_low}, \code{CI_high} The lower and upper credible interval limits for the parameters. } @@ -106,20 +107,22 @@ from each tail of the distribution and always include the median, the HDI is distributions. While this can be useful to better represent the credibility mass of a distribution, the HDI also has some limitations. See \code{\link[=spi]{spi()}} for details. -\cr \cr + The \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{\verb{95\%} or \verb{89\%} Credible Intervals (CI)} -are two reasonable ranges to characterize the uncertainty related to the estimation (see \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{here} for a discussion about the differences between these two values). -\cr +are two reasonable ranges to characterize the uncertainty related to the +estimation (see \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{here} +for a discussion about the differences between these two values). + The \verb{89\%} intervals (\code{ci = 0.89}) are deemed to be more stable than, for -instance, \verb{95\%} intervals (\cite{Kruschke, 2014}). An effective sample size +instance, \verb{95\%} intervals (\emph{Kruschke, 2014}). An effective sample size of at least 10.000 is recommended if one wants to estimate \verb{95\%} intervals -with high precision (\cite{Kruschke, 2014, p. 183ff}). Unfortunately, the +with high precision (\emph{Kruschke, 2014, p. 183ff}). Unfortunately, the default number of posterior samples for most Bayes packages (e.g., \code{rstanarm} or \code{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 \verb{95\%} threshold (\cite{McElreath, 2015}). -\cr +exceed the already unstable \verb{95\%} threshold (\emph{McElreath, 2015}). + However, \verb{95\%} has some \href{https://easystats.github.io/blog/posts/bayestestr_95/}{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 @@ -127,17 +130,17 @@ 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 \verb{95\%} CI than for lower ranges such as \verb{89\%}), which is a good thing in the context of the reproducibility crisis. -\cr \cr + A \verb{95\%} equal-tailed interval (ETI) has \verb{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 \link[=hdi]{HDI}, return similar results. -\cr \cr + 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. -\cr \cr + 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 diff --git a/man/ci.Rd b/man/ci.Rd index 2cc9e1b80..a80f688d8 100644 --- a/man/ci.Rd +++ b/man/ci.Rd @@ -60,14 +60,16 @@ ci(x, ...) \method{ci}{MCMCglmm}(x, ci = 0.95, method = "ETI", verbose = TRUE, ...) } \arguments{ -\item{x}{A \code{stanreg} or \code{brmsfit} model, or a vector representing a posterior distribution.} +\item{x}{A \code{stanreg} or \code{brmsfit} model, or a vector representing a posterior +distribution.} \item{...}{Currently not used.} \item{ci}{Value or vector of probability of the CI (between 0 and 1) -to be estimated. Default to \code{.95} (\verb{95\%}).} +to be estimated. Default to \code{0.95} (\verb{95\%}).} -\item{method}{Can be \link[=eti]{'ETI'} (default), \link[=hdi]{'HDI'}, \link[=bci]{'BCI'}, \link[=spi]{'SPI'} or \link[=si]{'SI'}.} +\item{method}{Can be \link[=eti]{"ETI"} (default), \link[=hdi]{"HDI"}, \link[=bci]{"BCI"}, +\link[=spi]{"SPI"} or \link[=si]{"SI"}.} \item{verbose}{Toggle off warnings.} @@ -89,7 +91,8 @@ be abbreviated. Only applies to \pkg{brms}-models.} \value{ A data frame with following columns: \itemize{ -\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a vector, this column is missing. +\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a +vector, this column is missing. \item \code{CI} The probability of the credible interval. \item \code{CI_low}, \code{CI_high} The lower and upper credible interval limits for the parameters. } @@ -108,9 +111,9 @@ for: \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 -\dQuote{Given any value in the interval and the background assumptions, -the data should not seem very surprising} (\cite{Gelman & Greenland 2019}). -\cr \cr +"Given any value in the interval and the background assumptions, +the data should not seem very surprising" (\emph{Gelman & Greenland 2019}). + 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}. } \examples{ @@ -143,7 +146,8 @@ ci(model, method = "HDI") \dontshow{\}) # examplesIf} } \references{ -Gelman A, Greenland S. Are confidence intervals better termed "uncertainty intervals"? BMJ 2019;l5381. 10.1136/bmj.l5381 +Gelman A, Greenland S. Are confidence intervals better termed "uncertainty +intervals"? BMJ 2019;l5381. 10.1136/bmj.l5381 } \seealso{ Other ci: diff --git a/man/cwi.Rd b/man/cwi.Rd index 253ad9166..6b3db0ae0 100644 --- a/man/cwi.Rd +++ b/man/cwi.Rd @@ -24,7 +24,8 @@ resemble the arguments of the \code{.numeric} or \code{.data.frame}methods.} \value{ A data frame with following columns: \itemize{ -\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a vector, this column is missing. +\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a +vector, this column is missing. \item \code{CI} The probability of the credible interval. \item \code{CI_low}, \code{CI_high} The lower and upper credible interval limits for the parameters. } diff --git a/man/describe_posterior.Rd b/man/describe_posterior.Rd index 31015e25f..0f1131e34 100644 --- a/man/describe_posterior.Rd +++ b/man/describe_posterior.Rd @@ -84,7 +84,7 @@ to the estimate(s) (\code{SD} and \code{MAD} for \code{mean} and \code{median}, Dispersion is not available for \code{"MAP"} or \code{"mode"} centrality indices.} \item{ci}{Value or vector of probability of the CI (between 0 and 1) -to be estimated. Default to \code{.95} (\verb{95\%}).} +to be estimated. Default to \code{0.95} (\verb{95\%}).} \item{ci_method}{The type of index used for Credible Interval. Can be \code{"ETI"} (default, see \code{\link[=eti]{eti()}}), \code{"HDI"} @@ -187,7 +187,7 @@ if (require("rstanarm") && require("emmeans")) { # emmeans estimates # ----------------------------------------------- - describe_posterior(suppressWarnings(emtrends(model, ~1, "wt"))) + describe_posterior(emtrends(model, ~1, "wt")) } # BayesFactor objects diff --git a/man/equivalence_test.Rd b/man/equivalence_test.Rd index 9a2b2fb00..36aabd050 100644 --- a/man/equivalence_test.Rd +++ b/man/equivalence_test.Rd @@ -38,17 +38,19 @@ equivalence_test(x, ...) ) } \arguments{ -\item{x}{Vector representing a posterior distribution. Can also be a \code{stanreg} or \code{brmsfit} model.} +\item{x}{Vector representing a posterior distribution. Can also be a +\code{stanreg} or \code{brmsfit} model.} \item{...}{Currently not used.} \item{range}{ROPE's lower and higher bounds. Should be \code{"default"} or -depending on the number of outcome variables a vector or a list. In models with one response, -\code{range} should be a vector of length two (e.g., \code{c(-0.1, 0.1)}). In -multivariate models, \code{range} should be a list with a numeric vectors for -each response variable. Vector names should correspond to the name of the response -variables. If \code{"default"} and input is a vector, the range is set to \code{c(-0.1, 0.1)}. If \code{"default"} and input is a Bayesian model, -\code{\link[=rope_range]{rope_range()}} is used.} +depending on the number of outcome variables a vector or a list. In +models with one response, \code{range} should be a vector of length two (e.g., +\code{c(-0.1, 0.1)}). In multivariate models, \code{range} should be a list with a +numeric vectors for each response variable. Vector names should correspond +to the name of the response variables. If \code{"default"} and input is a vector, +the range is set to \code{c(-0.1, 0.1)}. If \code{"default"} and input is a Bayesian +model, \code{\link[=rope_range]{rope_range()}} is used.} \item{ci}{The Credible Interval (CI) probability, corresponding to the proportion of HDI, to use for the percentage in ROPE.} @@ -152,9 +154,7 @@ test <- equivalence_test(model) plot(test) library(emmeans) -equivalence_test( - suppressWarnings(emtrends(model, ~1, "wt", data = mtcars)) -) +equivalence_test(emtrends(model, ~1, "wt", data = mtcars)) library(brms) model <- brms::brm(mpg ~ wt + cyl, data = mtcars) diff --git a/man/estimate_density.Rd b/man/estimate_density.Rd index 001a2e712..c6460eeaf 100644 --- a/man/estimate_density.Rd +++ b/man/estimate_density.Rd @@ -121,9 +121,7 @@ model <- suppressWarnings( head(estimate_density(model)) library(emmeans) -head(estimate_density(suppressWarnings( - emtrends(model, ~1, "wt", data = mtcars) -))) +head(estimate_density(emtrends(model, ~1, "wt", data = mtcars))) # brms models # ----------------------------------------------- diff --git a/man/eti.Rd b/man/eti.Rd index 3161a8a9c..9e7e7ce62 100644 --- a/man/eti.Rd +++ b/man/eti.Rd @@ -62,7 +62,8 @@ for the output.} \value{ A data frame with following columns: \itemize{ -\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a vector, this column is missing. +\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a +vector, this column is missing. \item \code{CI} The probability of the credible interval. \item \code{CI_low}, \code{CI_high} The lower and upper credible interval limits for the parameters. } @@ -81,20 +82,22 @@ from each tail of the distribution and always include the median, the HDI is distributions. While this can be useful to better represent the credibility mass of a distribution, the HDI also has some limitations. See \code{\link[=spi]{spi()}} for details. -\cr \cr + The \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{\verb{95\%} or \verb{89\%} Credible Intervals (CI)} -are two reasonable ranges to characterize the uncertainty related to the estimation (see \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{here} for a discussion about the differences between these two values). -\cr +are two reasonable ranges to characterize the uncertainty related to the +estimation (see \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{here} +for a discussion about the differences between these two values). + The \verb{89\%} intervals (\code{ci = 0.89}) are deemed to be more stable than, for -instance, \verb{95\%} intervals (\cite{Kruschke, 2014}). An effective sample size +instance, \verb{95\%} intervals (\emph{Kruschke, 2014}). An effective sample size of at least 10.000 is recommended if one wants to estimate \verb{95\%} intervals -with high precision (\cite{Kruschke, 2014, p. 183ff}). Unfortunately, the +with high precision (\emph{Kruschke, 2014, p. 183ff}). Unfortunately, the default number of posterior samples for most Bayes packages (e.g., \code{rstanarm} or \code{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 \verb{95\%} threshold (\cite{McElreath, 2015}). -\cr +exceed the already unstable \verb{95\%} threshold (\emph{McElreath, 2015}). + However, \verb{95\%} has some \href{https://easystats.github.io/blog/posts/bayestestr_95/}{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 @@ -102,17 +105,17 @@ 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 \verb{95\%} CI than for lower ranges such as \verb{89\%}), which is a good thing in the context of the reproducibility crisis. -\cr \cr + A \verb{95\%} equal-tailed interval (ETI) has \verb{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 \link[=hdi]{HDI}, return similar results. -\cr \cr + 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. -\cr \cr + 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 @@ -139,7 +142,7 @@ eti(model) eti(model, ci = c(0.80, 0.89, 0.95)) library(emmeans) -eti(suppressWarnings(emtrends(model, ~1, "wt", data = mtcars))) +eti(emtrends(model, ~1, "wt", data = mtcars)) library(brms) model <- brms::brm(mpg ~ wt + cyl, data = mtcars) diff --git a/man/hdi.Rd b/man/hdi.Rd index 2e5f9f95c..59c3919e9 100644 --- a/man/hdi.Rd +++ b/man/hdi.Rd @@ -65,7 +65,8 @@ for the output.} \value{ A data frame with following columns: \itemize{ -\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a vector, this column is missing. +\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a +vector, this column is missing. \item \code{CI} The probability of the credible interval. \item \code{CI_low}, \code{CI_high} The lower and upper credible interval limits for the parameters. } @@ -83,20 +84,22 @@ from each tail of the distribution and always include the median, the HDI is distributions. While this can be useful to better represent the credibility mass of a distribution, the HDI also has some limitations. See \code{\link[=spi]{spi()}} for details. -\cr \cr + The \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{\verb{95\%} or \verb{89\%} Credible Intervals (CI)} -are two reasonable ranges to characterize the uncertainty related to the estimation (see \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{here} for a discussion about the differences between these two values). -\cr +are two reasonable ranges to characterize the uncertainty related to the +estimation (see \href{https://easystats.github.io/bayestestR/articles/credible_interval.html}{here} +for a discussion about the differences between these two values). + The \verb{89\%} intervals (\code{ci = 0.89}) are deemed to be more stable than, for -instance, \verb{95\%} intervals (\cite{Kruschke, 2014}). An effective sample size +instance, \verb{95\%} intervals (\emph{Kruschke, 2014}). An effective sample size of at least 10.000 is recommended if one wants to estimate \verb{95\%} intervals -with high precision (\cite{Kruschke, 2014, p. 183ff}). Unfortunately, the +with high precision (\emph{Kruschke, 2014, p. 183ff}). Unfortunately, the default number of posterior samples for most Bayes packages (e.g., \code{rstanarm} or \code{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 \verb{95\%} threshold (\cite{McElreath, 2015}). -\cr +exceed the already unstable \verb{95\%} threshold (\emph{McElreath, 2015}). + However, \verb{95\%} has some \href{https://easystats.github.io/blog/posts/bayestestr_95/}{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 @@ -104,17 +107,17 @@ 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 \verb{95\%} CI than for lower ranges such as \verb{89\%}), which is a good thing in the context of the reproducibility crisis. -\cr \cr + A \verb{95\%} equal-tailed interval (ETI) has \verb{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 \link[=hdi]{HDI}, return similar results. -\cr \cr + 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. -\cr \cr + 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 @@ -143,9 +146,7 @@ bayestestR::hdi(model) bayestestR::hdi(model, ci = c(0.80, 0.90, 0.95)) library(emmeans) -bayestestR::hdi(suppressWarnings( - emtrends(model, ~1, "wt", data = mtcars) -)) +bayestestR::hdi(emtrends(model, ~1, "wt", data = mtcars)) library(brms) model <- brms::brm(mpg ~ wt + cyl, data = mtcars) diff --git a/man/mediation.Rd b/man/mediation.Rd index 7b749260d..de4d9643f 100644 --- a/man/mediation.Rd +++ b/man/mediation.Rd @@ -63,9 +63,10 @@ additional predictor), \code{y}, is not transformed. Then we could use \code{"mode"} or \code{"all"}.} \item{ci}{Value or vector of probability of the CI (between 0 and 1) -to be estimated. Default to \code{.95} (\verb{95\%}).} +to be estimated. Default to \code{0.95} (\verb{95\%}).} -\item{method}{Can be \link[=eti]{'ETI'} (default), \link[=hdi]{'HDI'}, \link[=bci]{'BCI'}, \link[=spi]{'SPI'} or \link[=si]{'SI'}.} +\item{method}{Can be \link[=eti]{"ETI"} (default), \link[=hdi]{"HDI"}, \link[=bci]{"BCI"}, +\link[=spi]{"SPI"} or \link[=si]{"SI"}.} } \value{ A data frame with direct, indirect, mediator and diff --git a/man/p_direction.Rd b/man/p_direction.Rd index f6282310c..cf90a3531 100644 --- a/man/p_direction.Rd +++ b/man/p_direction.Rd @@ -163,9 +163,7 @@ if (require("rstanarm")) { # emmeans # ----------------------------------------------- if (require("emmeans")) { - p_direction(suppressWarnings( - emtrends(model, ~1, "wt", data = mtcars) - )) + p_direction(emtrends(model, ~1, "wt", data = mtcars)) } # brms models diff --git a/man/p_rope.Rd b/man/p_rope.Rd index 169ad4e91..99515923c 100644 --- a/man/p_rope.Rd +++ b/man/p_rope.Rd @@ -31,17 +31,19 @@ p_rope(x, ...) ) } \arguments{ -\item{x}{Vector representing a posterior distribution. Can also be a \code{stanreg} or \code{brmsfit} model.} +\item{x}{Vector representing a posterior distribution. Can also be a +\code{stanreg} or \code{brmsfit} model.} \item{...}{Currently not used.} \item{range}{ROPE's lower and higher bounds. Should be \code{"default"} or -depending on the number of outcome variables a vector or a list. In models with one response, -\code{range} should be a vector of length two (e.g., \code{c(-0.1, 0.1)}). In -multivariate models, \code{range} should be a list with a numeric vectors for -each response variable. Vector names should correspond to the name of the response -variables. If \code{"default"} and input is a vector, the range is set to \code{c(-0.1, 0.1)}. If \code{"default"} and input is a Bayesian model, -\code{\link[=rope_range]{rope_range()}} is used.} +depending on the number of outcome variables a vector or a list. In +models with one response, \code{range} should be a vector of length two (e.g., +\code{c(-0.1, 0.1)}). In multivariate models, \code{range} should be a list with a +numeric vectors for each response variable. Vector names should correspond +to the name of the response variables. If \code{"default"} and input is a vector, +the range is set to \code{c(-0.1, 0.1)}. If \code{"default"} and input is a Bayesian +model, \code{\link[=rope_range]{rope_range()}} is used.} \item{effects}{Should results for fixed effects, random effects or both be returned? Only applies to mixed models. May be abbreviated.} diff --git a/man/p_significance.Rd b/man/p_significance.Rd index ce59b5bcc..2564e744b 100644 --- a/man/p_significance.Rd +++ b/man/p_significance.Rd @@ -33,7 +33,8 @@ p_significance(x, ...) ) } \arguments{ -\item{x}{Vector representing a posterior distribution. Can also be a \code{stanreg} or \code{brmsfit} model.} +\item{x}{Vector representing a posterior distribution. Can also be a +\code{stanreg} or \code{brmsfit} model.} \item{...}{Currently not used.} diff --git a/man/point_estimate.Rd b/man/point_estimate.Rd index 333b08359..635c3a627 100644 --- a/man/point_estimate.Rd +++ b/man/point_estimate.Rd @@ -99,7 +99,7 @@ point_estimate(model, centrality = c("median", "MAP")) # ----------------------------------------------- library(emmeans) point_estimate( - suppressWarnings(emtrends(model, ~1, "wt", data = mtcars)), + emtrends(model, ~1, "wt", data = mtcars), centrality = c("median", "MAP") ) diff --git a/man/rope.Rd b/man/rope.Rd index 86bfe31bc..5957e4f99 100644 --- a/man/rope.Rd +++ b/man/rope.Rd @@ -37,17 +37,19 @@ rope(x, ...) ) } \arguments{ -\item{x}{Vector representing a posterior distribution. Can also be a \code{stanreg} or \code{brmsfit} model.} +\item{x}{Vector representing a posterior distribution. Can also be a +\code{stanreg} or \code{brmsfit} model.} \item{...}{Currently not used.} \item{range}{ROPE's lower and higher bounds. Should be \code{"default"} or -depending on the number of outcome variables a vector or a list. In models with one response, -\code{range} should be a vector of length two (e.g., \code{c(-0.1, 0.1)}). In -multivariate models, \code{range} should be a list with a numeric vectors for -each response variable. Vector names should correspond to the name of the response -variables. If \code{"default"} and input is a vector, the range is set to \code{c(-0.1, 0.1)}. If \code{"default"} and input is a Bayesian model, -\code{\link[=rope_range]{rope_range()}} is used.} +depending on the number of outcome variables a vector or a list. In +models with one response, \code{range} should be a vector of length two (e.g., +\code{c(-0.1, 0.1)}). In multivariate models, \code{range} should be a list with a +numeric vectors for each response variable. Vector names should correspond +to the name of the response variables. If \code{"default"} and input is a vector, +the range is set to \code{c(-0.1, 0.1)}. If \code{"default"} and input is a Bayesian +model, \code{\link[=rope_range]{rope_range()}} is used.} \item{ci}{The Credible Interval (CI) probability, corresponding to the proportion of HDI, to use for the percentage in ROPE.} @@ -71,30 +73,36 @@ filtered by default, so only parameters that typically appear in the for the output.} } \description{ -Compute the proportion of the HDI (default to the \verb{89\%} HDI) of a posterior distribution that lies within a region of practical equivalence. +Compute the proportion of the HDI (default to the \verb{89\%} HDI) of a posterior +distribution that lies within a region of practical equivalence. } -\details{ -\subsection{ROPE}{ +\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}. +} +\section{ROPE}{ + Statistically, the probability of a posterior distribution of being different from 0 does not make much sense (the probability of a single value null hypothesis in a continuous distribution is 0). Therefore, the idea underlining ROPE is to let the user define an area around the null value enclosing values that are \emph{equivalent to the null} value for practical -purposes (\cite{Kruschke 2010, 2011, 2014}). -\cr \cr +purposes (\emph{Kruschke 2010, 2011, 2014}). + Kruschke (2018) suggests that such null value could be set, by default, to the -0.1 to 0.1 range of a standardized parameter (negligible effect size according to Cohen, 1988). This could be generalized: For instance, for linear models, the ROPE could be set as \verb{0 +/- .1 * sd(y)}. This ROPE range can be automatically computed for models using the \link{rope_range} function. -\cr \cr + Kruschke (2010, 2011, 2014) suggests using the proportion of the \verb{95\%} (or \verb{89\%}, considered more stable) \link[=hdi]{HDI} that falls within the ROPE as an index for "null-hypothesis" testing (as understood under the Bayesian framework, see \code{\link[=equivalence_test]{equivalence_test()}}). } -\subsection{Sensitivity to parameter's scale}{ + +\section{Sensitivity to parameter's scale}{ + It is important to consider the unit (i.e., the scale) of the predictors when using an index based on the ROPE, as the correct interpretation of the ROPE as representing a region of practical equivalence to zero is dependent @@ -103,7 +111,9 @@ the unit of its parameter. In other words, as the ROPE represents a fixed portion of the response's scale, its proximity with a coefficient depends on the scale of the coefficient itself. } -\subsection{Multicollinearity: Non-independent covariates}{ + +\section{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. Collinearity invalidates ROPE and hypothesis @@ -112,65 +122,81 @@ on independence. Most problematic are parameters that only have partial overlap with the ROPE region. In case of collinearity, the (joint) distributions of these parameters may either get an increased or decreased ROPE, which means that inferences based on \code{rope()} are inappropriate -(\cite{Kruschke 2014, 340f}). -\cr \cr +(\emph{Kruschke 2014, 340f}). + \code{rope()} 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 (\cite{Piironen and Vehtari 2017}). -} -\subsection{Strengths and Limitations}{ -\strong{Strengths:} Provides information related to the practical relevance of the effects. -\cr \cr -\strong{Limitations:} A ROPE range needs to be arbitrarily defined. Sensitive to the scale (the unit) of the predictors. Not sensitive to highly significant effects. -} +predictive variable selection (\emph{Piironen and Vehtari 2017}). } -\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}. + +\section{Strengths and Limitations}{ + +\strong{Strengths:} Provides information related to the practical relevance of +the effects. + +\strong{Limitations:} A ROPE range needs to be arbitrarily defined. Sensitive to +the scale (the unit) of the predictors. Not sensitive to highly significant +effects. } + \examples{ library(bayestestR) rope(x = rnorm(1000, 0, 0.01), range = c(-0.1, 0.1)) rope(x = rnorm(1000, 0, 1), range = c(-0.1, 0.1)) rope(x = rnorm(1000, 1, 0.01), range = c(-0.1, 0.1)) -rope(x = rnorm(1000, 1, 1), ci = c(.90, .95)) +rope(x = rnorm(1000, 1, 1), ci = c(0.90, 0.95)) \dontrun{ library(rstanarm) model <- suppressWarnings( stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0) ) rope(model) -rope(model, ci = c(.90, .95)) +rope(model, ci = c(0.90, 0.95)) library(emmeans) -rope(emtrends(model, ~1, "wt"), ci = c(.90, .95)) +rope(emtrends(model, ~1, "wt"), ci = c(0.90, 0.95)) library(brms) -model <- brms::brm(mpg ~ wt + cyl, data = mtcars) +model <- brm(mpg ~ wt + cyl, data = mtcars) rope(model) -rope(model, ci = c(.90, .95)) +rope(model, ci = c(0.90, 0.95)) library(brms) -model <- brms::brm(brms::mvbind(mpg, disp) ~ wt + cyl, data = mtcars) +model <- brm( + bf(mvbind(mpg, disp) ~ wt + cyl) + set_rescor(rescor = TRUE), + data = mtcars, + set_rescore +) rope(model) -rope(model, ci = c(.90, .95)) +rope(model, ci = c(0.90, 0.95)) library(BayesFactor) bf <- ttestBF(x = rnorm(100, 1, 1)) rope(bf) -rope(bf, ci = c(.90, .95)) +rope(bf, ci = c(0.90, 0.95)) } } \references{ \itemize{ \item Cohen, J. (1988). Statistical power analysis for the behavioural sciences. -\item Kruschke, J. K. (2010). What to believe: Bayesian methods for data analysis. Trends in cognitive sciences, 14(7), 293-300. \doi{10.1016/j.tics.2010.05.001}. -\item Kruschke, J. K. (2011). Bayesian assessment of null values via parameter estimation and model comparison. Perspectives on Psychological Science, 6(3), 299-312. \doi{10.1177/1745691611406925}. -\item Kruschke, J. K. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press. \doi{10.1177/2515245918771304}. -\item 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}. -\item Makowski D, Ben-Shachar MS, Chen SHA, 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} -\item 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} +\item Kruschke, J. K. (2010). What to believe: Bayesian methods for data analysis. +Trends in cognitive sciences, 14(7), 293-300. \doi{10.1016/j.tics.2010.05.001}. +\item Kruschke, J. K. (2011). Bayesian assessment of null values via parameter +estimation and model comparison. Perspectives on Psychological Science, +6(3), 299-312. \doi{10.1177/1745691611406925}. +\item Kruschke, J. K. (2014). Doing Bayesian data analysis: A tutorial with R, +JAGS, and Stan. Academic Press. \doi{10.1177/2515245918771304}. +\item 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}. +\item Makowski D, Ben-Shachar MS, Chen SHA, 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} +\item 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} } } diff --git a/man/sexit_thresholds.Rd b/man/sexit_thresholds.Rd index c4f2cfe17..08ff415ae 100644 --- a/man/sexit_thresholds.Rd +++ b/man/sexit_thresholds.Rd @@ -7,7 +7,8 @@ sexit_thresholds(x, ...) } \arguments{ -\item{x}{Vector representing a posterior distribution. Can also be a \code{stanreg} or \code{brmsfit} model.} +\item{x}{Vector representing a posterior distribution. Can also be a +\code{stanreg} or \code{brmsfit} model.} \item{...}{Currently not used.} } diff --git a/man/si.Rd b/man/si.Rd index b7d42c5eb..436470e5c 100644 --- a/man/si.Rd +++ b/man/si.Rd @@ -189,7 +189,8 @@ si(brms_model) \dontshow{\}) # examplesIf} } \references{ -Wagenmakers, E., Gronau, Q. F., Dablander, F., & Etz, A. (2018, November 22). The Support Interval. \doi{10.31234/osf.io/zwnxb} +Wagenmakers, E., Gronau, Q. F., Dablander, F., & Etz, A. (2018, November 22). +The Support Interval. \doi{10.31234/osf.io/zwnxb} } \seealso{ Other ci: diff --git a/man/spi.Rd b/man/spi.Rd index 085be3a30..f809085a6 100644 --- a/man/spi.Rd +++ b/man/spi.Rd @@ -62,7 +62,8 @@ for the output.} \value{ A data frame with following columns: \itemize{ -\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a vector, this column is missing. +\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a +vector, this column is missing. \item \code{CI} The probability of the credible interval. \item \code{CI_low}, \code{CI_high} The lower and upper credible interval limits for the parameters. }