From b9fa98a9d477aecb469e17083b3e87c4bfd51e06 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 18:23:05 +0200 Subject: [PATCH] Flexible ROPE values for describe_posterior Fixes #643 --- DESCRIPTION | 2 +- R/describe_posterior.R | 8 ++- R/rope.R | 71 +++++++++++++++++------- man/describe_posterior.Rd | 8 ++- man/equivalence_test.Rd | 14 ++--- man/p_rope.Rd | 14 ++--- man/rope.Rd | 17 +++--- tests/testthat/test-describe_posterior.R | 5 ++ 8 files changed, 91 insertions(+), 48 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 99a51aec8..b19f9c897 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: bayestestR Title: Understand and Describe Bayesian Models and Posterior Distributions -Version: 0.14.0.8 +Version: 0.14.0.9 Authors@R: c(person(given = "Dominique", family = "Makowski", diff --git a/R/describe_posterior.R b/R/describe_posterior.R index 80085750f..715501f31 100644 --- a/R/describe_posterior.R +++ b/R/describe_posterior.R @@ -18,9 +18,10 @@ #' For each "test", the corresponding \pkg{bayestestR} function is called #' (e.g. [bayestestR::rope()] or [bayestestR::p_direction()]) and its results #' included in the summary output. -#' @param rope_range ROPE's lower and higher bounds. Should be a list of two -#' values (e.g., `c(-0.1, 0.1)`) or `"default"`. If `"default"`, -#' the bounds are set to `x +- 0.1*SD(response)`. +#' @param rope_range ROPE's lower and higher bounds. Should be a vector of two +#' values (e.g., `c(-0.1, 0.1)`), `"default"` or a list of numeric vectors of +#' the same length as numbers of parameters. If `"default"`, the bounds are +#' set to `x +- 0.1*SD(response)`. #' @param rope_ci The Credible Interval (CI) probability, corresponding to the #' proportion of HDI, to use for the percentage in ROPE. #' @param keep_iterations If `TRUE`, will keep all iterations (draws) of @@ -90,6 +91,7 @@ #' 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 #' # ----------------------------------------------- diff --git a/R/rope.R b/R/rope.R index 782967a8f..38f02f53c 100644 --- a/R/rope.R +++ b/R/rope.R @@ -6,13 +6,14 @@ #' @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` can be a vector of length two (e.g., `c(-0.1, +#' 0.1)`), or a list of numeric vector of the same length as numbers of +#' parameters (see 'Examples'). 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. #' @param ci_method The type of interval to use to quantify the percentage in @@ -110,6 +111,9 @@ #' rope(model) #' rope(model, ci = c(0.90, 0.95)) #' +#' # multiple ROPE ranges +#' rope(model, range = list(c(-10, 5), c(-0.2, 0.2), "default")) +#' #' library(emmeans) #' rope(emtrends(model, ~1, "wt"), ci = c(0.90, 0.95)) #' @@ -381,7 +385,7 @@ rope.stanreg <- function(x, range = "default", ci = 0.95, ci_method = "ETI", eff if (all(range == "default")) { range <- rope_range(x, verbose = verbose) - } else if (!all(is.numeric(range)) || length(range) != 2) { + } else if (!is.list(range) && (!all(is.numeric(range)) || length(range) != 2)) { insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") } @@ -442,7 +446,7 @@ rope.brmsfit <- function(x, "With a multivariate model, `range` should be 'default' or a list of named numeric vectors with length 2." ) } - } else if (!all(is.numeric(range)) || length(range) != 2) { + } else if (!is.list(range) && (!all(is.numeric(range)) || length(range) != 2)) { insight::format_error( "`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1))." ) @@ -514,7 +518,7 @@ rope.sim.merMod <- function(x, if (all(range == "default")) { range <- rope_range(x, verbose = verbose) - } else if (!all(is.numeric(range)) || length(range) != 2) { + } else if (!is.list(range) && (!all(is.numeric(range)) || length(range) != 2)) { insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") } @@ -574,7 +578,7 @@ rope.sim.merMod <- function(x, rope.sim <- function(x, range = "default", ci = 0.95, ci_method = "ETI", parameters = NULL, verbose = TRUE, ...) { if (all(range == "default")) { range <- rope_range(x, verbose = verbose) - } else if (!all(is.numeric(range)) || length(range) != 2) { + } else if (!is.list(range) && (!all(is.numeric(range)) || length(range) != 2)) { insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") } @@ -607,15 +611,42 @@ 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( - parms, - rope, - range = range, - ci = ci, - ci_method = ci_method, - verbose = verbose, - simplify = FALSE - ) + if (is.list(range)) { + if (length(range) != ncol(parms)) { + insight::format_error("Length of `range` (i.e. number of ROPE limits) should match the number of parameters.") + } + # check if list of values contains only valid values + checks <- vapply(range, function(r) { + !all(r == "default") || !all(is.numeric(r)) || length(r) != 2 + }, logical(1)) + if (!all(checks)) { + insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") + } + tmp <- mapply( + function(p, r) { + rope( + p, + range = r, + ci = ci, + ci_method = ci_method, + verbose = verbose + ) + }, + parms, + range, + SIMPLIFY = FALSE + ) + } else { + tmp <- sapply( + parms, + rope, + range = range, + ci = ci, + ci_method = ci_method, + verbose = verbose, + simplify = FALSE + ) + } HDI_area <- lapply(tmp, attr, which = "HDI_area") diff --git a/man/describe_posterior.Rd b/man/describe_posterior.Rd index bc3e0fcfe..a81b821fa 100644 --- a/man/describe_posterior.Rd +++ b/man/describe_posterior.Rd @@ -121,9 +121,10 @@ For each "test", the corresponding \pkg{bayestestR} function is called (e.g. \code{\link[=rope]{rope()}} or \code{\link[=p_direction]{p_direction()}}) and its results included in the summary output.} -\item{rope_range}{ROPE's lower and higher bounds. Should be a list of two -values (e.g., \code{c(-0.1, 0.1)}) or \code{"default"}. If \code{"default"}, -the bounds are set to \code{x +- 0.1*SD(response)}.} +\item{rope_range}{ROPE's lower and higher bounds. Should be a vector of two +values (e.g., \code{c(-0.1, 0.1)}), \code{"default"} or a list of numeric vectors of +the same length as numbers of parameters. If \code{"default"}, the bounds are +set to \code{x +- 0.1*SD(response)}.} \item{rope_ci}{The Credible Interval (CI) probability, corresponding to the proportion of HDI, to use for the percentage in ROPE.} @@ -210,6 +211,7 @@ if (require("rstanarm") && require("emmeans")) { 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 # ----------------------------------------------- diff --git a/man/equivalence_test.Rd b/man/equivalence_test.Rd index accbbcaae..d292fb09b 100644 --- a/man/equivalence_test.Rd +++ b/man/equivalence_test.Rd @@ -51,13 +51,13 @@ equivalence_test(x, ...) \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} can be a vector of length two (e.g., \code{c(-0.1, 0.1)}), or a list of numeric vector of the same length as numbers of +parameters (see 'Examples'). 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.} diff --git a/man/p_rope.Rd b/man/p_rope.Rd index 8791c0224..aa76e2d1e 100644 --- a/man/p_rope.Rd +++ b/man/p_rope.Rd @@ -42,13 +42,13 @@ p_rope(x, ...) \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} can be a vector of length two (e.g., \code{c(-0.1, 0.1)}), or a list of numeric vector of the same length as numbers of +parameters (see 'Examples'). 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{verbose}{Toggle off warnings.} diff --git a/man/rope.Rd b/man/rope.Rd index 87036d460..cf82a12e8 100644 --- a/man/rope.Rd +++ b/man/rope.Rd @@ -54,13 +54,13 @@ rope(x, ...) \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} can be a vector of length two (e.g., \code{c(-0.1, 0.1)}), or a list of numeric vector of the same length as numbers of +parameters (see 'Examples'). 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.} @@ -171,6 +171,9 @@ model <- suppressWarnings( rope(model) rope(model, ci = c(0.90, 0.95)) +# multiple ROPE ranges +rope(model, range = list(c(-10, 5), c(-0.2, 0.2), "default")) + library(emmeans) rope(emtrends(model, ~1, "wt"), ci = c(0.90, 0.95)) diff --git a/tests/testthat/test-describe_posterior.R b/tests/testthat/test-describe_posterior.R index f21a10787..d53955d46 100644 --- a/tests/testthat/test-describe_posterior.R +++ b/tests/testthat/test-describe_posterior.R @@ -146,6 +146,11 @@ test_that("describe_posterior", { ) expect_identical(dim(rez), c(4L, 21L)) + # allow multiple ropes + rez <- describe_posterior(x, rope_range = list(c(-1, 1), "default")) + expect_identical(rez$ROPE_low, c(-1, -0.1), tolerance = 1e-3) + expect_identical(rez$ROPE_high, c(1, 0.1), tolerance = 1e-3) + rez <- describe_posterior( x, centrality = NULL,