From b9fa98a9d477aecb469e17083b3e87c4bfd51e06 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 18:23:05 +0200 Subject: [PATCH 01/22] 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, From 93de7a9f98ea64cef79be8617933303c844e2161 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 19:38:08 +0200 Subject: [PATCH 02/22] fix p_significance --- R/p_significance.R | 35 +++++++++++++++++++++++++++++++++-- R/print.R | 40 ++++++++++++++++++++++++++++------------ R/print_html.R | 18 +++++++++++++----- R/print_md.R | 18 +++++++++++++----- 4 files changed, 87 insertions(+), 24 deletions(-) diff --git a/R/p_significance.R b/R/p_significance.R index 46c8bd1c5..7c1e01898 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -133,6 +133,28 @@ p_significance.data.frame <- function(x, threshold = "default", rvar_col = NULL, if (ncol(x) == 1) { ps <- .p_significance(x[, 1], threshold = threshold, ...) + } else if (is.list(threshold)) { + if (length(threshold) != ncol(x)) { + insight::format_error("Length of `threshold` should match the number of parameters.") + } + # check if list of values contains only valid values + checks <- vapply(threshold, function(r) { + !all(r == "default") || !all(is.numeric(r)) || length(r) != 2 + }, logical(1)) + if (!all(checks)) { + insight::format_error("`threshold` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") + } + ps <- mapply( + function(p, thres) { + .p_significance( + p, + threshold = thres + ) + }, + x, + threshold, + SIMPLIFY = FALSE + ) } else { ps <- sapply(x, .p_significance, threshold = threshold, simplify = TRUE, ...) } @@ -365,6 +387,15 @@ as.double.p_significance <- as.numeric.p_significance #' @keywords internal .select_threshold_ps <- function(model = NULL, threshold = "default", verbose = TRUE) { + if (is.list(threshold)) { + lapply(threshold, .select_threshold_list, model = model, verbose = verbose) + } else { + .select_threshold_list(model = model, threshold = threshold, verbose = verbose) + } +} + +#' @keywords internal +.select_threshold_list <- function(model = NULL, threshold = "default", verbose = TRUE) { # If default if (all(threshold == "default")) { if (is.null(model)) { @@ -372,11 +403,11 @@ as.double.p_significance <- as.numeric.p_significance } else { threshold <- rope_range(model, verbose = verbose)[2] } - } else if (!all(is.numeric(threshold)) || length(threshold) > 2) { + } else if (!is.list(threshold) && (!all(is.numeric(threshold)) || length(threshold) > 2)) { insight::format_error( "`threshold` should be one of the following values:", "- \"default\", in which case the threshold is based on `rope_range()`", - "- a single numeric value (e.g., 0.1), which is used as range around zero (i.e. the threshold range is set to -0.1 and 0.1)", + "- a single numeric value (e.g., 0.1), which is used as range around zero (i.e. the threshold range is set to -0.1 and 0.1)", # nolint "- a numeric vector of length two (e.g., `c(-0.2, 0.1)`)" ) } diff --git a/R/print.R b/R/print.R index fbe3a2fd6..c4f464f58 100644 --- a/R/print.R +++ b/R/print.R @@ -70,14 +70,19 @@ print.map_estimate <- function(x, #' @export print.p_rope <- function(x, digits = 2, ...) { - caption <- sprintf( - "Proportion of samples inside the ROPE [%.*f, %.*f]", - digits, - x$ROPE_low[1], - digits, - x$ROPE_high[1] - ) - x$ROPE_low <- x$ROPE_high <- NULL + # check if we have multiple ROPE values + if (insight::n_unique(x$ROPE_low) > 1) { + caption <- "Proportion of samples inside the ROPE" + } else { + caption <- sprintf( + "Proportion of samples inside the ROPE [%.*f, %.*f]", + digits, + x$ROPE_low[1], + digits, + x$ROPE_high[1] + ) + x$ROPE_low <- x$ROPE_high <- NULL + } .print_default( x = x, digits = digits, @@ -90,14 +95,25 @@ print.p_rope <- function(x, digits = 2, ...) { #' @export print.p_significance <- function(x, digits = 2, ...) { - caption <- sprintf( - "Practical Significance (threshold: %s)", - insight::format_value(attributes(x)$threshold, digits = digits) - ) + threshold <- attributes(x)$threshold + if (is.list(threshold)) { + caption <- "Practical Significance" + out <- as.data.frame(do.call(rbind, threshold)) + colnames(out) <- c("ROPE_low", "ROPE_high") + x <- cbind(x, out) + ci_string <- "ROPE" + } else { + caption <- sprintf( + "Practical Significance (threshold: %s)", + insight::format_value(attributes(x)$threshold, digits = digits) + ) + ci_string <- NULL + } .print_default( x = x, digits = digits, caption = caption, + ci_string = ci_string, ... ) } diff --git a/R/print_html.R b/R/print_html.R index a7a4be75c..b9f958814 100644 --- a/R/print_html.R +++ b/R/print_html.R @@ -31,11 +31,19 @@ print_html.p_map <- function(x, digits = 2, caption = "MAP-based p-value", ...) #' @export print_html.p_rope <- function(x, digits = 2, ...) { - caption <- sprintf( - "Proportion of samples inside the ROPE [%.*f, %.*f]", - digits, x$ROPE_low[1], digits, x$ROPE_high[1] - ) - x$ROPE_low <- x$ROPE_high <- NULL + # check if we have multiple ROPE values + if (insight::n_unique(x$ROPE_low) > 1) { + caption <- "Proportion of samples inside the ROPE" + } else { + caption <- sprintf( + "Proportion of samples inside the ROPE [%.*f, %.*f]", + digits, + x$ROPE_low[1], + digits, + x$ROPE_high[1] + ) + x$ROPE_low <- x$ROPE_high <- NULL + } .print_html_default(x = x, digits = digits, caption = caption, ci_string = "ROPE", ...) } diff --git a/R/print_md.R b/R/print_md.R index f8910a123..643c7d50a 100644 --- a/R/print_md.R +++ b/R/print_md.R @@ -30,11 +30,19 @@ print_md.p_map <- function(x, digits = 2, caption = "MAP-based p-value", ...) { #' @export print_md.p_rope <- function(x, digits = 2, ...) { - caption <- sprintf( - "Proportion of samples inside the ROPE [%.*f, %.*f]", - digits, x$ROPE_low[1], digits, x$ROPE_high[1] - ) - x$ROPE_low <- x$ROPE_high <- NULL + # check if we have multiple ROPE values + if (insight::n_unique(x$ROPE_low) > 1) { + caption <- "Proportion of samples inside the ROPE" + } else { + caption <- sprintf( + "Proportion of samples inside the ROPE [%.*f, %.*f]", + digits, + x$ROPE_low[1], + digits, + x$ROPE_high[1] + ) + x$ROPE_low <- x$ROPE_high <- NULL + } .print_md_default(x = x, digits = digits, caption = caption, ci_string = "ROPE", ...) } From 81cbdbcf593daf1455ac36e4b670d028e0ee7177 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 20:16:07 +0200 Subject: [PATCH 03/22] fix print for p_significance --- R/print.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/print.R b/R/print.R index c4f464f58..6fbae4c9c 100644 --- a/R/print.R +++ b/R/print.R @@ -100,7 +100,8 @@ print.p_significance <- function(x, digits = 2, ...) { caption <- "Practical Significance" out <- as.data.frame(do.call(rbind, threshold)) colnames(out) <- c("ROPE_low", "ROPE_high") - x <- cbind(x, out) + x$ROPE_low <- out$ROPE_low + x$ROPE_high <- out$ROPE_high ci_string <- "ROPE" } else { caption <- sprintf( From 0789ab28200e5f03c7867f03830a9f23504ff9a6 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 20:17:25 +0200 Subject: [PATCH 04/22] docs --- R/p_significance.R | 5 ++++- man/p_significance.Rd | 5 ++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/R/p_significance.R b/R/p_significance.R index 7c1e01898..e9027d420 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -15,7 +15,8 @@ #' (i.e. the threshold range is set to -0.1 and 0.1, i.e. reflects a symmetric #' interval) #' - a numeric vector of length two (e.g., `c(-0.2, 0.1)`), useful for -#' asymmetric intervals. +#' asymmetric intervals +#' - a list of numeric vectors, where each vector corresponds to a parameter. #' @inheritParams rope #' @inheritParams hdi #' @@ -53,6 +54,8 @@ #' chains = 2, refresh = 0 #' ) #' p_significance(model) +#' # multiple thresholds +#' p_significance(model, threshold = list(c(-10, 5), c(-0.2, 0.2), "default")) #' } #' @export p_significance <- function(x, ...) { diff --git a/man/p_significance.Rd b/man/p_significance.Rd index 30b2170e5..1e4ad290c 100644 --- a/man/p_significance.Rd +++ b/man/p_significance.Rd @@ -59,7 +59,8 @@ and based on \code{\link[=rope_range]{rope_range()}} if a (Bayesian) model is pr (i.e. the threshold range is set to -0.1 and 0.1, i.e. reflects a symmetric interval) \item a numeric vector of length two (e.g., \code{c(-0.2, 0.1)}), useful for -asymmetric intervals. +asymmetric intervals +\item a list of numeric vectors, where each vector corresponds to a parameter. }} \item{use_iterations}{Logical, if \code{TRUE} and \code{x} is a \code{get_predicted} object, @@ -131,6 +132,8 @@ model <- rstanarm::stan_glm(mpg ~ wt + cyl, chains = 2, refresh = 0 ) p_significance(model) +# multiple thresholds +p_significance(model, threshold = list(c(-10, 5), c(-0.2, 0.2), "default")) } \dontshow{\}) # examplesIf} } From 5b45fe18cba8b02422c5aac5fce27e605e0b8576 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 20:18:42 +0200 Subject: [PATCH 05/22] lintr --- R/eti.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/eti.R b/R/eti.R index 4bcf75af1..404a9f036 100644 --- a/R/eti.R +++ b/R/eti.R @@ -316,8 +316,8 @@ eti.get_predicted <- function(x, ci = 0.95, use_iterations = FALSE, verbose = TR )) data.frame( - "CI" = ci, - "CI_low" = results[1], - "CI_high" = results[2] + CI = ci, + CI_low = results[1], + CI_high = results[2] ) } From ca77ee844d186a8234f631aaba684ec3ac45b4f4 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 20:21:04 +0200 Subject: [PATCH 06/22] print methods --- R/print_html.R | 21 ++++++++++++++++----- R/print_md.R | 21 ++++++++++++++++----- 2 files changed, 32 insertions(+), 10 deletions(-) diff --git a/R/print_html.R b/R/print_html.R index b9f958814..6e9f7883a 100644 --- a/R/print_html.R +++ b/R/print_html.R @@ -50,11 +50,22 @@ print_html.p_rope <- function(x, digits = 2, ...) { #' @export print_html.p_significance <- function(x, digits = 2, ...) { - caption <- sprintf( - "Practical Significance (threshold: %s)", - insight::format_value(attributes(x)$threshold, digits = digits) - ) - .print_html_default(x = x, digits = digits, caption = caption, ...) + threshold <- attributes(x)$threshold + if (is.list(threshold)) { + caption <- "Practical Significance" + out <- as.data.frame(do.call(rbind, threshold)) + colnames(out) <- c("ROPE_low", "ROPE_high") + x$ROPE_low <- out$ROPE_low + x$ROPE_high <- out$ROPE_high + ci_string <- "ROPE" + } else { + caption <- sprintf( + "Practical Significance (threshold: %s)", + insight::format_value(attributes(x)$threshold, digits = digits) + ) + ci_string <- NULL + } + .print_html_default(x = x, digits = digits, caption = caption, ci_string = ci_string, ...) } diff --git a/R/print_md.R b/R/print_md.R index 643c7d50a..c45cd27b2 100644 --- a/R/print_md.R +++ b/R/print_md.R @@ -49,11 +49,22 @@ print_md.p_rope <- function(x, digits = 2, ...) { #' @export print_md.p_significance <- function(x, digits = 2, ...) { - caption <- sprintf( - "Practical Significance (threshold: %s)", - insight::format_value(attributes(x)$threshold, digits = digits) - ) - .print_md_default(x = x, digits = digits, caption = caption, ...) + threshold <- attributes(x)$threshold + if (is.list(threshold)) { + caption <- "Practical Significance" + out <- as.data.frame(do.call(rbind, threshold)) + colnames(out) <- c("ROPE_low", "ROPE_high") + x$ROPE_low <- out$ROPE_low + x$ROPE_high <- out$ROPE_high + ci_string <- "ROPE" + } else { + caption <- sprintf( + "Practical Significance (threshold: %s)", + insight::format_value(attributes(x)$threshold, digits = digits) + ) + ci_string <- NULL + } + .print_md_default(x = x, digits = digits, caption = caption, ci_string = ci_string, ...) } From 22566e07f97920ce1e2a3aba00b12d1f0f35e844 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 21:24:58 +0200 Subject: [PATCH 07/22] add test --- R/p_rope.R | 2 +- R/p_significance.R | 14 ++++++++++---- man/p_significance.Rd | 4 ++-- tests/testthat/test-p_significance.R | 13 ++++++++++++- 4 files changed, 25 insertions(+), 8 deletions(-) diff --git a/R/p_rope.R b/R/p_rope.R index 3a8ae940f..d3531186e 100644 --- a/R/p_rope.R +++ b/R/p_rope.R @@ -229,7 +229,7 @@ p_rope.mcmc.list <- p_rope.mcmc #' @keywords internal .p_rope <- function(rope_rez) { cols <- c("Parameter", "ROPE_low", "ROPE_high", "ROPE_Percentage", "Effects", "Component") - out <- as.data.frame(rope_rez[cols[cols %in% names(rope_rez)]]) + out <- as.data.frame(rope_rez)[cols[cols %in% names(rope_rez)]] names(out)[names(out) == "ROPE_Percentage"] <- "p_ROPE" class(out) <- c("p_rope", "see_p_rope", "data.frame") diff --git a/R/p_significance.R b/R/p_significance.R index e9027d420..c00917612 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -54,8 +54,8 @@ #' chains = 2, refresh = 0 #' ) #' p_significance(model) -#' # multiple thresholds -#' p_significance(model, threshold = list(c(-10, 5), c(-0.2, 0.2), "default")) +#' # multiple thresholds - asymmetric, symmetric, default +#' p_significance(model, threshold = list(c(-10, 5), 0.2, "default")) #' } #' @export p_significance <- function(x, ...) { @@ -142,7 +142,7 @@ p_significance.data.frame <- function(x, threshold = "default", rvar_col = NULL, } # check if list of values contains only valid values checks <- vapply(threshold, function(r) { - !all(r == "default") || !all(is.numeric(r)) || length(r) != 2 + !all(r == "default") || !all(is.numeric(r)) || length(r) > 2 }, logical(1)) if (!all(checks)) { insight::format_error("`threshold` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") @@ -391,7 +391,13 @@ as.double.p_significance <- as.numeric.p_significance #' @keywords internal .select_threshold_ps <- function(model = NULL, threshold = "default", verbose = TRUE) { if (is.list(threshold)) { - lapply(threshold, .select_threshold_list, model = model, verbose = verbose) + lapply(threshold, function(i) { + out <- .select_threshold_list(model = model, threshold = i, verbose = verbose) + if (length(out) == 1) { + out <- c(-1 * abs(out), abs(out)) + } + out + }) } else { .select_threshold_list(model = model, threshold = threshold, verbose = verbose) } diff --git a/man/p_significance.Rd b/man/p_significance.Rd index 1e4ad290c..78e10864c 100644 --- a/man/p_significance.Rd +++ b/man/p_significance.Rd @@ -132,8 +132,8 @@ model <- rstanarm::stan_glm(mpg ~ wt + cyl, chains = 2, refresh = 0 ) p_significance(model) -# multiple thresholds -p_significance(model, threshold = list(c(-10, 5), c(-0.2, 0.2), "default")) +# multiple thresholds - asymmetric, symmetric, default +p_significance(model, threshold = list(c(-10, 5), 0.2, "default")) } \dontshow{\}) # examplesIf} } diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R index 6f542f3c0..cd4ccc95f 100644 --- a/tests/testthat/test-p_significance.R +++ b/tests/testthat/test-p_significance.R @@ -49,7 +49,7 @@ test_that("stanreg", { test_that("brms", { skip_if_offline() - skip_if_not_or_load_if_installed("rstanarm") + skip_if_not_or_load_if_installed("brms") m2 <- insight::download_model("brms_1") @@ -58,4 +58,15 @@ test_that("brms", { c(1.0000, 0.9985, 0.9785), tolerance = 0.01 ) + + expect_equal( + p_significance(m2, threshold = list(1, "default", 2), effects = "all")$ps, + c(1.00000, 0.99850, 0.12275), + tolerance = 0.01 + ) + expect_equal( + attributes(out)$threshold, + list(c(-1, 1), c(-0.60269480520891, 0.60269480520891), c(-2, 2)), + tolerance = 1e-4 + ) }) From b70ecbf9b71b1f6bb7c7f6c867ed64c0141310f0 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 22:20:03 +0200 Subject: [PATCH 08/22] Update test-p_significance.R --- tests/testthat/test-p_significance.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R index cd4ccc95f..0c4f2e488 100644 --- a/tests/testthat/test-p_significance.R +++ b/tests/testthat/test-p_significance.R @@ -59,8 +59,9 @@ test_that("brms", { tolerance = 0.01 ) + out <- p_significance(m2, threshold = list(1, "default", 2), effects = "all") expect_equal( - p_significance(m2, threshold = list(1, "default", 2), effects = "all")$ps, + out$ps, c(1.00000, 0.99850, 0.12275), tolerance = 0.01 ) From 9ed77e1d3f2b7d079538c18d31dfc726808bd0c6 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 08:42:04 +0200 Subject: [PATCH 09/22] news, tests --- NEWS.md | 10 +++++++--- tests/testthat/test-describe_posterior.R | 9 +++++++++ tests/testthat/test-p_rope.R | 20 ++++++++++++++++++++ tests/testthat/test-p_significance.R | 9 +++++++++ tests/testthat/test-rope.R | 9 +++++++++ 5 files changed, 54 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/test-p_rope.R diff --git a/NEWS.md b/NEWS.md index db28c58f5..f98458f3f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,14 +2,18 @@ ## Changes -* Support for `posterior::rvar`-type column in data frames. - For example, a data frame `df` with an `rvar` column `".pred"` can now be +* Support for `posterior::rvar`-type column in data frames. + For example, a data frame `df` with an `rvar` column `".pred"` can now be called directly via `p_direction(df, rvar_col = ".pred")`. * Added support for `{marginaleffects}` +* The ROPE or threshold ranges in `rope()`, `describe_posterior()` and + `p_significance()` can now be specified as a list. This allows for different + ranges for different parameters. + * Results from objects generated by `{emmeans}` (`emmGrid`/`emm_list`) now - return results with appended grid-data. + return results with appended grid-data. * Usability improvements for `p_direction()`: diff --git a/tests/testthat/test-describe_posterior.R b/tests/testthat/test-describe_posterior.R index d53955d46..db265a679 100644 --- a/tests/testthat/test-describe_posterior.R +++ b/tests/testthat/test-describe_posterior.R @@ -151,6 +151,15 @@ test_that("describe_posterior", { expect_identical(rez$ROPE_low, c(-1, -0.1), tolerance = 1e-3) expect_identical(rez$ROPE_high, c(1, 0.1), tolerance = 1e-3) + expect_error( + describe_posterior(x, rope_range = list(1, "default")), + regex = "should be 'default'" + ) + expect_error( + describe_posterior(x, rope_range = list(c(1, 1), c(2, 2), c(2, 3))), + regex = "Length of" + ) + rez <- describe_posterior( x, centrality = NULL, diff --git a/tests/testthat/test-p_rope.R b/tests/testthat/test-p_rope.R new file mode 100644 index 000000000..827614443 --- /dev/null +++ b/tests/testthat/test-p_rope.R @@ -0,0 +1,20 @@ +test_that("rope", { + skip_if_offline() + skip_if_not_or_load_if_installed("rstanarm") + m <- insight::download_model("stanreg_merMod_5") + p <- insight::get_parameters(m, effects = "all") + expect_equal( + p_rope(as.data.frame(m)[2:4], range = list(c(0, 40), "default", c(-1, 0.8)))$p_ROPE, + c(0.598, 0.002, 0.396), + tolerance = 1e-3 + ) + + expect_error( + p_rope(as.data.frame(m)[2:4], range = list(c(0, 40), c(-1, 0.8))), + regex = "Length of" + ) + expect_error( + p_rope(as.data.frame(m)[2:4], range = list(c(0, 40), "a", c(-1, 0.8))), + regex = "should be 'default'" + ) +}) diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R index 0c4f2e488..a03d42e0c 100644 --- a/tests/testthat/test-p_significance.R +++ b/tests/testthat/test-p_significance.R @@ -70,4 +70,13 @@ test_that("brms", { list(c(-1, 1), c(-0.60269480520891, 0.60269480520891), c(-2, 2)), tolerance = 1e-4 ) + + expect_error( + p_significance(m2, threshold = list(1, "a", 2), effects = "all"), + regex = "should be one of" + ) + expect_error( + p_significance(m2, threshold = list(1, 2, 3, 4), effects = "all"), + regex = "Length of" + ) }) diff --git a/tests/testthat/test-rope.R b/tests/testthat/test-rope.R index 8367e3d80..1c9012da7 100644 --- a/tests/testthat/test-rope.R +++ b/tests/testthat/test-rope.R @@ -49,6 +49,15 @@ test_that("rope", { rope(p, verbose = FALSE)$ROPE_Percentage, tolerance = 1e-3 ) + + expect_error( + rope(m, range = list(c(-0.1, 0.1), c(2, 2))), + regex = "Length of" + ) + expect_error( + rope(m, range = list(c(-0.1, 0.1), c(2, 2), "default", "a", c(1, 3))), + regex = "should be 'default'" + ) }) From 983a7190b2667963ba8c7ad7b4769c590ce35906 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 09:10:39 +0200 Subject: [PATCH 10/22] equi_test --- R/equivalence_test.R | 46 +++++++++++++++++++++----- R/print.equivalence_test.R | 20 ++++++----- man/equivalence_test.Rd | 2 ++ tests/testthat/test-equivalence_test.R | 32 ++++++++++++++++++ tests/testthat/test-p_rope.R | 3 +- 5 files changed, 83 insertions(+), 20 deletions(-) create mode 100644 tests/testthat/test-equivalence_test.R diff --git a/R/equivalence_test.R b/R/equivalence_test.R index 57d4b2058..05aa2c62d 100644 --- a/R/equivalence_test.R +++ b/R/equivalence_test.R @@ -85,6 +85,8 @@ #' \donttest{ #' 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")) #' #' # plot result #' test <- equivalence_test(model) @@ -244,7 +246,7 @@ equivalence_test.BFBayesFactor <- function(x, range = "default", ci = 0.95, verb verbose = TRUE) { if (all(range == "default")) { range <- rope_range(x, verbose = verbose) - } else if (!all(is.numeric(range)) || length(range) != 2L) { + } else if (!is.list(range) && (!all(is.numeric(range)) || length(range) != 2L)) { insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") } @@ -257,14 +259,40 @@ equivalence_test.BFBayesFactor <- function(x, range = "default", ci = 0.95, verb verbose = verbose ) - l <- sapply( - params, - equivalence_test, - range = range, - ci = ci, - verbose = verbose, - simplify = FALSE - ) + if (is.list(range)) { + if (length(range) != ncol(params)) { + 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)).") + } + l <- mapply( + function(p, r) { + equivalence_test( + p, + range = r, + ci = ci, + verbose = verbose + ) + }, + params, + range, + SIMPLIFY = FALSE + ) + } else { + l <- sapply( + params, + equivalence_test, + range = range, + ci = ci, + verbose = verbose, + simplify = FALSE + ) + } dat <- do.call(rbind, l) out <- data.frame( diff --git a/R/print.equivalence_test.R b/R/print.equivalence_test.R index 1b9d03cac..9bf61fdb6 100644 --- a/R/print.equivalence_test.R +++ b/R/print.equivalence_test.R @@ -2,7 +2,10 @@ print.equivalence_test <- function(x, digits = 2, ...) { orig_x <- x insight::print_color("# Test for Practical Equivalence\n\n", "blue") - cat(sprintf(" ROPE: [%.*f %.*f]\n\n", digits, x$ROPE_low[1], digits, x$ROPE_high[1])) + # print ROPE limits, if we just have one set of ROPE values + if (insight::n_unique(x$ROPE_low) == 1) { + cat(sprintf(" ROPE: [%.*f %.*f]\n\n", digits, x$ROPE_low[1], digits, x$ROPE_high[1])) + } # fix "sd" pattern model <- .retrieve_model(x) @@ -23,15 +26,8 @@ print.equivalence_test <- function(x, digits = 2, ...) { } } - # find the longest HDI-value, so we can align the brackets in the ouput - x$HDI_low <- sprintf("%.*f", digits, x$HDI_low) - x$HDI_high <- sprintf("%.*f", digits, x$HDI_high) - - maxlen_low <- max(nchar(x$HDI_low)) - maxlen_high <- max(nchar(x$HDI_high)) - x$ROPE_Percentage <- sprintf("%.*f %%", digits, x$ROPE_Percentage * 100) - x$HDI <- sprintf("[%*s %*s]", maxlen_low, x$HDI_low, maxlen_high, x$HDI_high) + x$HDI <- insight::format_ci(x$HDI_low, x$HDI_high, ci = NULL, digits = digits) ci <- unique(x$CI) keep.columns <- c( @@ -39,6 +35,12 @@ print.equivalence_test <- function(x, digits = 2, ...) { "ROPE_Equivalence", "ROPE_Percentage", "CI", "HDI" ) + # keep ROPE columns for multiple ROPE values + if (insight::n_unique(x$ROPE_low) > 1) { + keep.columns <- c(keep.columns, "ROPE") + x$ROPE <- insight::format_ci(x$ROPE_low, x$ROPE_high, ci = NULL, digits = digits) + } + x <- x[, intersect(keep.columns, colnames(x))] colnames(x)[which(colnames(x) == "ROPE_Equivalence")] <- "H0" diff --git a/man/equivalence_test.Rd b/man/equivalence_test.Rd index d292fb09b..ea0b3a182 100644 --- a/man/equivalence_test.Rd +++ b/man/equivalence_test.Rd @@ -164,6 +164,8 @@ print(test, digits = 4) \donttest{ 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, 5), c(-0.2, 0.2), "default")) # plot result test <- equivalence_test(model) diff --git a/tests/testthat/test-equivalence_test.R b/tests/testthat/test-equivalence_test.R new file mode 100644 index 000000000..632e58807 --- /dev/null +++ b/tests/testthat/test-equivalence_test.R @@ -0,0 +1,32 @@ +test_that("equivalence test", { + skip_if_offline() + skip_if_not_or_load_if_installed("rstanarm") + m <- insight::download_model("stanreg_merMod_5") + + out <- equivalence_test(m, verbose = FALSE) + expect_snapshot(print(out)) + + out <- equivalence_test( + m, + range = list(c(-1, 1), "default", c(0, 2), c(-2, 0), "default"), + verbose = FALSE + ) + expect_snapshot(print(out)) + + expect_error( + equivalence_test( + m, + range = list(c(-1, 1), "default", c(0, 2), c(-2, 0)), + verbose = FALSE + ), + regex = "Length of" + ) + expect_error( + equivalence_test( + m, + range = list(c(-1, 1), "default", c(0, 2), c(-2, 0), "a"), + verbose = FALSE + ), + regex = "should be 'default'" + ) +}) diff --git a/tests/testthat/test-p_rope.R b/tests/testthat/test-p_rope.R index 827614443..854722b76 100644 --- a/tests/testthat/test-p_rope.R +++ b/tests/testthat/test-p_rope.R @@ -1,8 +1,7 @@ -test_that("rope", { +test_that("p_rope", { skip_if_offline() skip_if_not_or_load_if_installed("rstanarm") m <- insight::download_model("stanreg_merMod_5") - p <- insight::get_parameters(m, effects = "all") expect_equal( p_rope(as.data.frame(m)[2:4], range = list(c(0, 40), "default", c(-1, 0.8)))$p_ROPE, c(0.598, 0.002, 0.396), From e171da2561286f3adbb9d43a732b33f3a5c7bce4 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 09:13:55 +0200 Subject: [PATCH 11/22] fix test --- tests/testthat/_snaps/equivalence_test.md | 36 +++++++++++++++++++++++ tests/testthat/test-equivalence_test.R | 2 +- 2 files changed, 37 insertions(+), 1 deletion(-) create mode 100644 tests/testthat/_snaps/equivalence_test.md diff --git a/tests/testthat/_snaps/equivalence_test.md b/tests/testthat/_snaps/equivalence_test.md new file mode 100644 index 000000000..8e54d8ac0 --- /dev/null +++ b/tests/testthat/_snaps/equivalence_test.md @@ -0,0 +1,36 @@ +# equivalence test + + Code + print(out) + Output + # Test for Practical Equivalence + + ROPE: [-0.18 0.18] + + Parameter | H0 | inside ROPE | 95% HDI + ----------------------------------------------------- + (Intercept) | Rejected | 0.00 % | [-2.68, -0.50] + size | Accepted | 100.00 % | [-0.04, 0.07] + period2 | Rejected | 0.00 % | [-1.61, -0.36] + period3 | Rejected | 0.00 % | [-1.77, -0.40] + period4 | Rejected | 0.00 % | [-2.52, -0.76] + + + +--- + + Code + print(out) + Output + # Test for Practical Equivalence + + Parameter | H0 | inside ROPE | 95% HDI | ROPE + ---------------------------------------------------------------------- + (Intercept) | Undecided | 15.82 % | [-2.68, -0.50] | [-1.00, 1.00] + size | Accepted | 100.00 % | [-0.04, 0.07] | [-0.10, 0.10] + period2 | Rejected | 0.00 % | [-1.61, -0.36] | [0.00, 2.00] + period3 | Accepted | 100.00 % | [-1.77, -0.40] | [-2.00, 0.00] + period4 | Rejected | 0.00 % | [-2.52, -0.76] | [-0.10, 0.10] + + + diff --git a/tests/testthat/test-equivalence_test.R b/tests/testthat/test-equivalence_test.R index 632e58807..bf0f83706 100644 --- a/tests/testthat/test-equivalence_test.R +++ b/tests/testthat/test-equivalence_test.R @@ -1,6 +1,6 @@ test_that("equivalence test", { skip_if_offline() - skip_if_not_or_load_if_installed("rstanarm") + # skip_if_not_or_load_if_installed("rstanarm") m <- insight::download_model("stanreg_merMod_5") out <- equivalence_test(m, verbose = FALSE) From 968d0eefcc642f106701627607616dec0a94b264 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 09:52:49 +0200 Subject: [PATCH 12/22] fix --- tests/testthat/test-equivalence_test.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-equivalence_test.R b/tests/testthat/test-equivalence_test.R index bf0f83706..632e58807 100644 --- a/tests/testthat/test-equivalence_test.R +++ b/tests/testthat/test-equivalence_test.R @@ -1,6 +1,6 @@ test_that("equivalence test", { skip_if_offline() - # skip_if_not_or_load_if_installed("rstanarm") + skip_if_not_or_load_if_installed("rstanarm") m <- insight::download_model("stanreg_merMod_5") out <- equivalence_test(m, verbose = FALSE) From 8de01b4f058572457feef9b5190fba87985c70f4 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 10:03:37 +0200 Subject: [PATCH 13/22] update --- NEWS.md | 4 +- R/equivalence_test.R | 86 +++++++++-------------- man/equivalence_test.Rd | 2 +- tests/testthat/_snaps/equivalence_test.md | 38 +++++++++- tests/testthat/test-equivalence_test.R | 39 +++++++++- 5 files changed, 113 insertions(+), 56 deletions(-) diff --git a/NEWS.md b/NEWS.md index f98458f3f..c3232190d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,8 +8,8 @@ * Added support for `{marginaleffects}` -* The ROPE or threshold ranges in `rope()`, `describe_posterior()` and - `p_significance()` can now be specified as a list. This allows for different +* The ROPE or threshold ranges in `rope()`, `describe_posterior()`, `p_significance()` + and `equivalence_test()` can now be specified as a list. This allows for different ranges for different parameters. * Results from objects generated by `{emmeans}` (`emmGrid`/`emm_list`) now diff --git a/R/equivalence_test.R b/R/equivalence_test.R index 05aa2c62d..589b8d7fc 100644 --- a/R/equivalence_test.R +++ b/R/equivalence_test.R @@ -165,13 +165,40 @@ equivalence_test.data.frame <- function(x, range = "default", ci = 0.95, rvar_co return(.append_datagrid(out, x)) } - l <- insight::compact_list(lapply( - x, - equivalence_test, - range = range, - ci = ci, - verbose = verbose - )) + # multiple ranges for the parameters - iterate over parameters and range + if (is.list(range)) { + if (length(range) != ncol(x)) { + 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)).") + } + l <- insight::compact_list(mapply( + function(p, r) { + equivalence_test( + p, + range = r, + ci = ci, + verbose = verbose + ) + }, + x, + range, + SIMPLIFY = FALSE + )) + } else { + l <- insight::compact_list(lapply( + x, + equivalence_test, + range = range, + ci = ci, + verbose = verbose + )) + } dat <- do.call(rbind, l) out <- data.frame( @@ -259,50 +286,7 @@ equivalence_test.BFBayesFactor <- function(x, range = "default", ci = 0.95, verb verbose = verbose ) - if (is.list(range)) { - if (length(range) != ncol(params)) { - 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)).") - } - l <- mapply( - function(p, r) { - equivalence_test( - p, - range = r, - ci = ci, - verbose = verbose - ) - }, - params, - range, - SIMPLIFY = FALSE - ) - } else { - l <- sapply( - params, - equivalence_test, - range = range, - ci = ci, - verbose = verbose, - simplify = FALSE - ) - } - - dat <- do.call(rbind, l) - out <- data.frame( - Parameter = rep(names(l), each = nrow(dat) / length(l)), - dat, - stringsAsFactors = FALSE - ) - - class(out) <- unique(c("equivalence_test", "see_equivalence_test", class(out))) - out + equivalence_test(params, range = range, ci = ci, verbose = verbose) } diff --git a/man/equivalence_test.Rd b/man/equivalence_test.Rd index ea0b3a182..81f5d42fd 100644 --- a/man/equivalence_test.Rd +++ b/man/equivalence_test.Rd @@ -165,7 +165,7 @@ 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, 5), c(-0.2, 0.2), "default")) +equivalence_test(model, range = list(c(10, 40), c(-5, -4), "default")) # plot result test <- equivalence_test(model) diff --git a/tests/testthat/_snaps/equivalence_test.md b/tests/testthat/_snaps/equivalence_test.md index 8e54d8ac0..cf351f609 100644 --- a/tests/testthat/_snaps/equivalence_test.md +++ b/tests/testthat/_snaps/equivalence_test.md @@ -1,4 +1,4 @@ -# equivalence test +# equivalence test, rstanarm Code print(out) @@ -34,3 +34,39 @@ +# equivalence test, df + + Code + print(out) + Output + # Test for Practical Equivalence + + ROPE: [-0.10 0.10] + + Parameter | H0 | inside ROPE | 95% HDI + ----------------------------------------------------- + (Intercept) | Rejected | 0.00 % | [-2.68, -0.50] + size | Accepted | 100.00 % | [-0.04, 0.07] + period2 | Rejected | 0.00 % | [-1.61, -0.36] + period3 | Rejected | 0.00 % | [-1.77, -0.40] + period4 | Rejected | 0.00 % | [-2.52, -0.76] + + + +--- + + Code + print(out) + Output + # Test for Practical Equivalence + + Parameter | H0 | inside ROPE | 95% HDI | ROPE + ---------------------------------------------------------------------- + (Intercept) | Undecided | 15.82 % | [-2.68, -0.50] | [-1.00, 1.00] + size | Accepted | 100.00 % | [-0.04, 0.07] | [-0.10, 0.10] + period2 | Rejected | 0.00 % | [-1.61, -0.36] | [0.00, 2.00] + period3 | Accepted | 100.00 % | [-1.77, -0.40] | [-2.00, 0.00] + period4 | Rejected | 0.00 % | [-2.52, -0.76] | [-0.10, 0.10] + + + diff --git a/tests/testthat/test-equivalence_test.R b/tests/testthat/test-equivalence_test.R index 632e58807..e7ce74cac 100644 --- a/tests/testthat/test-equivalence_test.R +++ b/tests/testthat/test-equivalence_test.R @@ -1,4 +1,6 @@ -test_that("equivalence test", { +skip_on_cran() + +test_that("equivalence test, rstanarm", { skip_if_offline() skip_if_not_or_load_if_installed("rstanarm") m <- insight::download_model("stanreg_merMod_5") @@ -30,3 +32,38 @@ test_that("equivalence test", { regex = "should be 'default'" ) }) + + +test_that("equivalence test, df", { + skip_if_offline() + skip_if_not_or_load_if_installed("rstanarm") + m <- insight::download_model("stanreg_merMod_5") + params <- as.data.frame(m)[1:5] + + out <- equivalence_test(params, verbose = FALSE) + expect_snapshot(print(out)) + + out <- equivalence_test( + params, + range = list(c(-1, 1), "default", c(0, 2), c(-2, 0), "default"), + verbose = FALSE + ) + expect_snapshot(print(out)) + + expect_error( + equivalence_test( + params, + range = list(c(-1, 1), "default", c(0, 2), c(-2, 0)), + verbose = FALSE + ), + regex = "Length of" + ) + expect_error( + equivalence_test( + params, + range = list(c(-1, 1), "default", c(0, 2), c(-2, 0), "a"), + verbose = FALSE + ), + regex = "should be 'default'" + ) +}) From 81fd43233cb60bcca3857c98e22a5d1e640d5431 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 10:12:17 +0200 Subject: [PATCH 14/22] DRY --- R/equivalence_test.R | 11 ++--------- R/p_significance.R | 28 +++++++++++++++++++--------- R/rope.R | 11 ++--------- 3 files changed, 23 insertions(+), 27 deletions(-) diff --git a/R/equivalence_test.R b/R/equivalence_test.R index 589b8d7fc..fa089ed75 100644 --- a/R/equivalence_test.R +++ b/R/equivalence_test.R @@ -167,16 +167,9 @@ equivalence_test.data.frame <- function(x, range = "default", ci = 0.95, rvar_co # multiple ranges for the parameters - iterate over parameters and range if (is.list(range)) { - if (length(range) != ncol(x)) { - 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)).") - } + .check_list_range(range, x) + # apply thresholds to each column l <- insight::compact_list(mapply( function(p, r) { equivalence_test( diff --git a/R/p_significance.R b/R/p_significance.R index c00917612..af075cef7 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -137,16 +137,9 @@ p_significance.data.frame <- function(x, threshold = "default", rvar_col = NULL, if (ncol(x) == 1) { ps <- .p_significance(x[, 1], threshold = threshold, ...) } else if (is.list(threshold)) { - if (length(threshold) != ncol(x)) { - insight::format_error("Length of `threshold` should match the number of parameters.") - } # check if list of values contains only valid values - checks <- vapply(threshold, function(r) { - !all(r == "default") || !all(is.numeric(r)) || length(r) > 2 - }, logical(1)) - if (!all(checks)) { - insight::format_error("`threshold` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") - } + .check_list_range(threshold, x, larger_two = TRUE) + # apply thresholds to each column ps <- mapply( function(p, thres) { .p_significance( @@ -422,3 +415,20 @@ as.double.p_significance <- as.numeric.p_significance } threshold } + +.check_list_range <- function(range, params, larger_two = FALSE) { + if (length(range) != ncol(params)) { + 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) { + if (larger_two) { + !all(r == "default") || !all(is.numeric(r)) || length(r) > 2 + } else { + !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)).") + } +} diff --git a/R/rope.R b/R/rope.R index 38f02f53c..a07726504 100644 --- a/R/rope.R +++ b/R/rope.R @@ -612,16 +612,9 @@ 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) { 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)).") - } + .check_list_range(range, params) + # apply thresholds to each column tmp <- mapply( function(p, r) { rope( From 8dd0b65d0d65af40fbb99785dabf516e6dc6707d Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 11:23:33 +0200 Subject: [PATCH 15/22] fix typo --- R/rope.R | 2 +- tests/testthat/test-bayesian_as_frequentist.R | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/rope.R b/R/rope.R index a07726504..8158e6091 100644 --- a/R/rope.R +++ b/R/rope.R @@ -613,7 +613,7 @@ rope.sim <- function(x, range = "default", ci = 0.95, ci_method = "ETI", paramet .prepare_rope_df <- function(parms, range, ci, ci_method, verbose) { if (is.list(range)) { # check if list of values contains only valid values - .check_list_range(range, params) + .check_list_range(range, parms) # apply thresholds to each column tmp <- mapply( function(p, r) { diff --git a/tests/testthat/test-bayesian_as_frequentist.R b/tests/testthat/test-bayesian_as_frequentist.R index 8066f3483..0fafdbcb9 100644 --- a/tests/testthat/test-bayesian_as_frequentist.R +++ b/tests/testthat/test-bayesian_as_frequentist.R @@ -34,7 +34,7 @@ test_that("brms beta to freq", { skip_if_not_or_load_if_installed("betareg") set.seed(333) - m <- insight::download_model("brms_beta_1") + m <- suppressWarnings(insight::download_model("brms_beta_1")) data(FoodExpenditure, package = "betareg") m1 <- glmmTMB::glmmTMB( I(food / income) ~ income + (1 | persons), @@ -56,7 +56,7 @@ test_that("ordbetareg to freq", { set.seed(333) data(sleepstudy, package = "lme4") - m <- insight::download_model("ordbetareg_1") + m <- suppressWarnings(insight::download_model("ordbetareg_1")) sleepstudy$y <- datawizard::normalize(sleepstudy$Reaction) m1 <- glmmTMB::glmmTMB( y ~ Days + (Days | Subject), From 2f17e457d10a5e9c8397f9688e4be6b6a067527f Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 11:38:24 +0200 Subject: [PATCH 16/22] fix --- R/equivalence_test.R | 2 +- R/p_significance.R | 21 ++++++++++++++++++++- R/rope.R | 2 +- tests/testthat/test-rope.R | 14 ++++++++++++++ 4 files changed, 36 insertions(+), 3 deletions(-) diff --git a/R/equivalence_test.R b/R/equivalence_test.R index fa089ed75..daec8dc79 100644 --- a/R/equivalence_test.R +++ b/R/equivalence_test.R @@ -168,7 +168,7 @@ equivalence_test.data.frame <- function(x, range = "default", ci = 0.95, rvar_co # multiple ranges for the parameters - iterate over parameters and range if (is.list(range)) { # check if list of values contains only valid values - .check_list_range(range, x) + range <- .check_list_range(range, x) # apply thresholds to each column l <- insight::compact_list(mapply( function(p, r) { diff --git a/R/p_significance.R b/R/p_significance.R index af075cef7..f2f9a34e2 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -138,7 +138,7 @@ p_significance.data.frame <- function(x, threshold = "default", rvar_col = NULL, ps <- .p_significance(x[, 1], threshold = threshold, ...) } else if (is.list(threshold)) { # check if list of values contains only valid values - .check_list_range(threshold, x, larger_two = TRUE) + threshold <- .check_list_range(threshold, x, larger_two = TRUE) # apply thresholds to each column ps <- mapply( function(p, thres) { @@ -417,6 +417,24 @@ as.double.p_significance <- as.numeric.p_significance } .check_list_range <- function(range, params, larger_two = FALSE) { + # if we have named elements, complete list + named_range <- names(range) + if (!is.null(named_range)) { + # find out which name belongs to which parameter + pos <- match(named_range, colnames(params)) + # if not all element names were found, error + if (anyNA(pos)) { + insight::format_error(paste( + "Not all elements of `range` were found in the parameters. Please check following names:", + toString(names_range[is.na]) + )) + } + # now "fill" non-specified elements with "default" + out <- rep("default", ncol(params)) + out[pos] <- range + # overwrite former range + range <- out + } if (length(range) != ncol(params)) { insight::format_error("Length of `range` (i.e. number of ROPE limits) should match the number of parameters.") } @@ -431,4 +449,5 @@ as.double.p_significance <- as.numeric.p_significance if (!all(checks)) { insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") } + range } diff --git a/R/rope.R b/R/rope.R index 8158e6091..61eec4d32 100644 --- a/R/rope.R +++ b/R/rope.R @@ -613,7 +613,7 @@ rope.sim <- function(x, range = "default", ci = 0.95, ci_method = "ETI", paramet .prepare_rope_df <- function(parms, range, ci, ci_method, verbose) { if (is.list(range)) { # check if list of values contains only valid values - .check_list_range(range, parms) + range <- .check_list_range(range, parms) # apply thresholds to each column tmp <- mapply( function(p, r) { diff --git a/tests/testthat/test-rope.R b/tests/testthat/test-rope.R index 1c9012da7..48e71e7aa 100644 --- a/tests/testthat/test-rope.R +++ b/tests/testthat/test-rope.R @@ -50,6 +50,20 @@ test_that("rope", { tolerance = 1e-3 ) + # list range + expect_equal( + rope(m, range = list(c(-1, 0.1), "default", "default", c(-1, 1), c(-1.5, -1)))$ROPE_Percentage, + c(0.15823, 1, 0, 0.3903, 0.38186) + tolerance = 1e-3 + ) + + # named elements, chooses "default" for unnamed + expect_equal( + rope(m, range = list(c(-1, 0.1), "default", "default", c(-1, 1), c(-1.5, -1)))$ROPE_Percentage, + rope(m, range = list("(Intercept)" = c(-1, 0.1), period4 = c(-1.5, -1), period3 = c(-1, 1)))$ROPE_Percentage, + tolerance = 1e-3 + ) + expect_error( rope(m, range = list(c(-0.1, 0.1), c(2, 2))), regex = "Length of" From d18afa44b5ab228c1ed50dc1d9632d2651048541 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 12:07:01 +0200 Subject: [PATCH 17/22] allow named vectors --- R/p_significance.R | 57 +++++++++++++++++++++------- R/rope.R | 25 +++++++----- man/equivalence_test.Rd | 23 +++++++---- man/p_rope.Rd | 23 +++++++---- man/p_significance.Rd | 7 +++- man/rope.Rd | 25 ++++++++---- tests/testthat/test-p_significance.R | 27 ++++++++++++- tests/testthat/test-rope.R | 4 ++ 8 files changed, 145 insertions(+), 46 deletions(-) diff --git a/R/p_significance.R b/R/p_significance.R index f2f9a34e2..6edd70c8d 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -16,7 +16,10 @@ #' interval) #' - a numeric vector of length two (e.g., `c(-0.2, 0.1)`), useful for #' asymmetric intervals -#' - a list of numeric vectors, where each vector corresponds to a parameter. +#' - a list of numeric vectors, where each vector corresponds to a parameter +#' - a list of *named* numeric vectors, where names correspond to parameter +#' names. In this case, all parameters that have no matching name in `threshold` +#' will be set to `"default"`. #' @inheritParams rope #' @inheritParams hdi #' @@ -56,6 +59,8 @@ #' p_significance(model) #' # multiple thresholds - asymmetric, symmetric, default #' p_significance(model, threshold = list(c(-10, 5), 0.2, "default")) +#' # names thresholds +#' p_significance(model, threshold = list(wt = 0.2, `(Intercept)` = c(-10, 5))) #' } #' @export p_significance <- function(x, ...) { @@ -131,7 +136,7 @@ p_significance.data.frame <- function(x, threshold = "default", rvar_col = NULL, } - threshold <- .select_threshold_ps(threshold = threshold) + threshold <- .select_threshold_ps(threshold = threshold, params = x) x <- .select_nums(x) if (ncol(x) == 1) { @@ -286,12 +291,15 @@ p_significance.stanreg <- function(x, ...) { effects <- match.arg(effects) component <- match.arg(component) - threshold <- .select_threshold_ps(model = x, threshold = threshold, verbose = verbose) + params <- insight::get_parameters(x, effects = effects, component = component, parameters = parameters) - result <- p_significance( - insight::get_parameters(x, effects = effects, component = component, parameters = parameters), - threshold = threshold + threshold <- .select_threshold_ps( + model = x, + threshold = threshold, + params = params, + verbose = verbose ) + result <- p_significance(params, threshold = threshold) cleaned_parameters <- insight::clean_parameters(x) out <- .prepare_output(result, cleaned_parameters, inherits(x, "stanmvreg")) @@ -322,12 +330,15 @@ p_significance.brmsfit <- function(x, ...) { effects <- match.arg(effects) component <- match.arg(component) - threshold <- .select_threshold_ps(model = x, threshold = threshold, verbose = verbose) + params <- insight::get_parameters(x, effects = effects, component = component, parameters = parameters) - result <- p_significance( - insight::get_parameters(x, effects = effects, component = component, parameters = parameters), - threshold = threshold + threshold <- .select_threshold_ps( + model = x, + threshold = threshold, + params = params, + verbose = verbose ) + result <- p_significance(params, threshold = threshold) cleaned_parameters <- insight::clean_parameters(x) out <- .prepare_output(result, cleaned_parameters) @@ -382,8 +393,28 @@ as.double.p_significance <- as.numeric.p_significance # helpers -------------------------- #' @keywords internal -.select_threshold_ps <- function(model = NULL, threshold = "default", verbose = TRUE) { +.select_threshold_ps <- function(model = NULL, threshold = "default", params = NULL, verbose = TRUE) { if (is.list(threshold)) { + # if we have named elements, complete list + if (!is.null(params)) { + named_threshold <- names(threshold) + if (!is.null(named_threshold)) { + # find out which name belongs to which parameter + pos <- match(named_threshold, colnames(params)) + # if not all element names were found, error + if (anyNA(pos)) { + insight::format_error(paste( + "Not all elements of `threshold` were found in the parameters. Please check following names:", + toString(named_threshold[is.na(pos)]) + )) + } + # now "fill" non-specified elements with "default" + out <- as.list(rep("default", ncol(params))) + out[pos] <- threshold + # overwrite former threshold + threshold <- out + } + } lapply(threshold, function(i) { out <- .select_threshold_list(model = model, threshold = i, verbose = verbose) if (length(out) == 1) { @@ -426,11 +457,11 @@ as.double.p_significance <- as.numeric.p_significance if (anyNA(pos)) { insight::format_error(paste( "Not all elements of `range` were found in the parameters. Please check following names:", - toString(names_range[is.na]) + toString(named_range[is.na(pos)]) )) } # now "fill" non-specified elements with "default" - out <- rep("default", ncol(params)) + out <- as.list(rep("default", ncol(params))) out[pos] <- range # overwrite former range range <- out diff --git a/R/rope.R b/R/rope.R index 61eec4d32..b678c24eb 100644 --- a/R/rope.R +++ b/R/rope.R @@ -6,14 +6,21 @@ #' @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` 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. +#' depending on the number of outcome variables a vector or a list. For models +#' with one response, `range` can be: +#' +#' - a vector of length two (e.g., `c(-0.1, 0.1)`), +#' - a list of numeric vector of the same length as numbers of parameters (see +#' 'Examples'). +#' - a list of *named* numeric vectors, where names correspond to parameter +#' names. In this case, all parameters that have no matching name in `range` +#' will be set to `"default"`. +#' +#' 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()`] 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 @@ -34,7 +41,7 @@ #' 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. +#' [`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 diff --git a/man/equivalence_test.Rd b/man/equivalence_test.Rd index 81f5d42fd..adc86a51e 100644 --- a/man/equivalence_test.Rd +++ b/man/equivalence_test.Rd @@ -51,13 +51,22 @@ 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} 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.} +depending on the number of outcome variables a vector or a list. For models +with one response, \code{range} can be: +\itemize{ +\item a vector of length two (e.g., \code{c(-0.1, 0.1)}), +\item a list of numeric vector of the same length as numbers of parameters (see +'Examples'). +\item a list of \emph{named} numeric vectors, where names correspond to parameter +names. In this case, all parameters that have no matching name in \code{range} +will be set to \code{"default"}. +} + +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 aa76e2d1e..4d27fe375 100644 --- a/man/p_rope.Rd +++ b/man/p_rope.Rd @@ -42,13 +42,22 @@ 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} 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.} +depending on the number of outcome variables a vector or a list. For models +with one response, \code{range} can be: +\itemize{ +\item a vector of length two (e.g., \code{c(-0.1, 0.1)}), +\item a list of numeric vector of the same length as numbers of parameters (see +'Examples'). +\item a list of \emph{named} numeric vectors, where names correspond to parameter +names. In this case, all parameters that have no matching name in \code{range} +will be set to \code{"default"}. +} + +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/p_significance.Rd b/man/p_significance.Rd index 78e10864c..38abc8d2e 100644 --- a/man/p_significance.Rd +++ b/man/p_significance.Rd @@ -60,7 +60,10 @@ and based on \code{\link[=rope_range]{rope_range()}} if a (Bayesian) model is pr interval) \item a numeric vector of length two (e.g., \code{c(-0.2, 0.1)}), useful for asymmetric intervals -\item a list of numeric vectors, where each vector corresponds to a parameter. +\item a list of numeric vectors, where each vector corresponds to a parameter +\item a list of \emph{named} numeric vectors, where names correspond to parameter +names. In this case, all parameters that have no matching name in \code{threshold} +will be set to \code{"default"}. }} \item{use_iterations}{Logical, if \code{TRUE} and \code{x} is a \code{get_predicted} object, @@ -134,6 +137,8 @@ model <- rstanarm::stan_glm(mpg ~ wt + cyl, p_significance(model) # multiple thresholds - asymmetric, symmetric, default p_significance(model, threshold = list(c(-10, 5), 0.2, "default")) +# names thresholds +p_significance(model, threshold = list(wt = 0.2, `(Intercept)` = c(-10, 5))) } \dontshow{\}) # examplesIf} } diff --git a/man/rope.Rd b/man/rope.Rd index cf82a12e8..21f960185 100644 --- a/man/rope.Rd +++ b/man/rope.Rd @@ -54,13 +54,22 @@ 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} 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.} +depending on the number of outcome variables a vector or a list. For models +with one response, \code{range} can be: +\itemize{ +\item a vector of length two (e.g., \code{c(-0.1, 0.1)}), +\item a list of numeric vector of the same length as numbers of parameters (see +'Examples'). +\item a list of \emph{named} numeric vectors, where names correspond to parameter +names. In this case, all parameters that have no matching name in \code{range} +will be set to \code{"default"}. +} + +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.} @@ -107,7 +116,7 @@ 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. +\code{\link[=rope_range]{rope_range()}} function. 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 diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R index a03d42e0c..dbfae00a5 100644 --- a/tests/testthat/test-p_significance.R +++ b/tests/testthat/test-p_significance.R @@ -37,7 +37,6 @@ test_that("p_significance", { test_that("stanreg", { skip_if_offline() skip_if_not_or_load_if_installed("rstanarm") - m <- insight::download_model("stanreg_merMod_5") expect_equal( @@ -80,3 +79,29 @@ test_that("brms", { regex = "Length of" ) }) + +test_that("stan", { + skip_if_offline() + skip_if_not_or_load_if_installed("rstanarm") + m <- insight::download_model("stanreg_merMod_5") + + expect_equal( + p_significance(m, threshold = list("(Intercept)" = 1, period4 = 1.5, period3 = 0.5))$ps, + p_significance(m, threshold = list(1, "default", "default", 0.5, 1.5))$ps, + tolerance = 1e-4 + ) + + expect_error( + p_significance(m, threshold = list("(Intercept)" = 1, point = 1.5, period3 = 0.5)), + regex = "Not all elements" + ) + expect_error( + p_significance(m, threshold = list(1, "a", 2), effects = "all"), + regex = "should be one of" + ) + expect_error( + p_significance(m, threshold = list(1, 2, 3, 4), effects = "all"), + regex = "Length of" + ) + +}) diff --git a/tests/testthat/test-rope.R b/tests/testthat/test-rope.R index 48e71e7aa..c7c8eead5 100644 --- a/tests/testthat/test-rope.R +++ b/tests/testthat/test-rope.R @@ -72,6 +72,10 @@ test_that("rope", { rope(m, range = list(c(-0.1, 0.1), c(2, 2), "default", "a", c(1, 3))), regex = "should be 'default'" ) + expect_error( + rope(m, range = list("(Intercept)" = c(-1, 0.1), pointout = c(-1.5, -1), period3 = c(-1, 1))), + regex = "Not all elements" + ) }) From ccfa2ef9295c7dbde17cbe9f22a218abcb681bd8 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 12:18:07 +0200 Subject: [PATCH 18/22] fix --- R/equivalence_test.R | 2 ++ R/p_significance.R | 2 +- R/rope.R | 3 +++ man/equivalence_test.Rd | 2 ++ man/p_significance.Rd | 2 +- man/rope.Rd | 3 +++ tests/testthat/test-rope.R | 2 +- 7 files changed, 13 insertions(+), 3 deletions(-) diff --git a/R/equivalence_test.R b/R/equivalence_test.R index daec8dc79..638af6102 100644 --- a/R/equivalence_test.R +++ b/R/equivalence_test.R @@ -87,6 +87,8 @@ #' 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) diff --git a/R/p_significance.R b/R/p_significance.R index 6edd70c8d..12a692c1b 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -59,7 +59,7 @@ #' p_significance(model) #' # multiple thresholds - asymmetric, symmetric, default #' p_significance(model, threshold = list(c(-10, 5), 0.2, "default")) -#' # names thresholds +#' # named thresholds #' p_significance(model, threshold = list(wt = 0.2, `(Intercept)` = c(-10, 5))) #' } #' @export diff --git a/R/rope.R b/R/rope.R index b678c24eb..24c3a9b4a 100644 --- a/R/rope.R +++ b/R/rope.R @@ -121,6 +121,9 @@ #' # multiple ROPE ranges #' rope(model, range = list(c(-10, 5), c(-0.2, 0.2), "default")) #' +#' # named ROPE ranges +#' rope(model, range = list(gear = c(-3, 2), wt = c(-0.2, 0.2))) +#' #' library(emmeans) #' rope(emtrends(model, ~1, "wt"), ci = c(0.90, 0.95)) #' diff --git a/man/equivalence_test.Rd b/man/equivalence_test.Rd index adc86a51e..27fc778f7 100644 --- a/man/equivalence_test.Rd +++ b/man/equivalence_test.Rd @@ -175,6 +175,8 @@ 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) diff --git a/man/p_significance.Rd b/man/p_significance.Rd index 38abc8d2e..5d6d0fbaf 100644 --- a/man/p_significance.Rd +++ b/man/p_significance.Rd @@ -137,7 +137,7 @@ model <- rstanarm::stan_glm(mpg ~ wt + cyl, p_significance(model) # multiple thresholds - asymmetric, symmetric, default p_significance(model, threshold = list(c(-10, 5), 0.2, "default")) -# names thresholds +# named thresholds p_significance(model, threshold = list(wt = 0.2, `(Intercept)` = c(-10, 5))) } \dontshow{\}) # examplesIf} diff --git a/man/rope.Rd b/man/rope.Rd index 21f960185..242edb880 100644 --- a/man/rope.Rd +++ b/man/rope.Rd @@ -183,6 +183,9 @@ rope(model, ci = c(0.90, 0.95)) # multiple ROPE ranges rope(model, range = list(c(-10, 5), c(-0.2, 0.2), "default")) +# named ROPE ranges +rope(model, range = list(gear = c(-3, 2), wt = c(-0.2, 0.2))) + library(emmeans) rope(emtrends(model, ~1, "wt"), ci = c(0.90, 0.95)) diff --git a/tests/testthat/test-rope.R b/tests/testthat/test-rope.R index c7c8eead5..0c9059854 100644 --- a/tests/testthat/test-rope.R +++ b/tests/testthat/test-rope.R @@ -53,7 +53,7 @@ test_that("rope", { # list range expect_equal( rope(m, range = list(c(-1, 0.1), "default", "default", c(-1, 1), c(-1.5, -1)))$ROPE_Percentage, - c(0.15823, 1, 0, 0.3903, 0.38186) + c(0.15823, 1, 0, 0.3903, 0.38186), tolerance = 1e-3 ) From a02342b49e01a64558ee4e8f8879ff98b63b2a06 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 12:22:31 +0200 Subject: [PATCH 19/22] Update test-p_significance.R --- tests/testthat/test-p_significance.R | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R index dbfae00a5..d11ecbb25 100644 --- a/tests/testthat/test-p_significance.R +++ b/tests/testthat/test-p_significance.R @@ -103,5 +103,4 @@ test_that("stan", { p_significance(m, threshold = list(1, 2, 3, 4), effects = "all"), regex = "Length of" ) - }) From cb967ba1f0ac0f0e1e502063681b4363dc3ccc39 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 13:48:01 +0200 Subject: [PATCH 20/22] fix test warning --- tests/testthat/test-check_prior.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-check_prior.R b/tests/testthat/test-check_prior.R index 50028c4e2..4b2512e8c 100644 --- a/tests/testthat/test-check_prior.R +++ b/tests/testthat/test-check_prior.R @@ -92,7 +92,7 @@ test_that("check_prior - brms (linux)", { expect_warning(expect_identical( check_prior(model2, method = "lakeland")$Prior_Quality, c( - "informative", "misinformative", "informative", "informative", + "informative", "informative", "informative", "informative", "informative", "not determinable", "not determinable", "not determinable" ) )) From bf2ce499f022c7234f51166bcd42aa242d9b2804 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 13:51:35 +0200 Subject: [PATCH 21/22] lintr --- R/print.R | 2 +- R/print_html.R | 4 ++-- R/print_md.R | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/R/print.R b/R/print.R index 6fbae4c9c..dfc62eeb7 100644 --- a/R/print.R +++ b/R/print.R @@ -319,7 +319,7 @@ print.bayesfactor_parameters <- function(x, digits = 3, log = FALSE, ...) { sep = " ", header = NULL, format = "text", - align = align, + align = align )) invisible(x) diff --git a/R/print_html.R b/R/print_html.R index 6e9f7883a..fb487c405 100644 --- a/R/print_html.R +++ b/R/print_html.R @@ -105,7 +105,7 @@ print_html.bayesfactor_models <- function(x, log = log, show_names = show_names, caption = caption, - align = c("llr"), + align = "llr", ... ) } @@ -122,7 +122,7 @@ print_html.bayesfactor_inclusion <- function(x, digits = digits, log = log, caption = caption, - align = c("lrrr"), + align = "lrrr", ... ) } diff --git a/R/print_md.R b/R/print_md.R index c45cd27b2..053ae0472 100644 --- a/R/print_md.R +++ b/R/print_md.R @@ -104,7 +104,7 @@ print_md.bayesfactor_models <- function(x, log = log, show_names = show_names, caption = caption, - align = c("llr"), + align = "llr", ... ) } @@ -121,7 +121,7 @@ print_md.bayesfactor_inclusion <- function(x, digits = digits, log = log, caption = caption, - align = c("lrrr"), + align = "lrrr", ... ) } From ab60d467fbd109bf42e13742f7247b170f6c8670 Mon Sep 17 00:00:00 2001 From: Daniel Date: Tue, 17 Sep 2024 14:42:44 +0200 Subject: [PATCH 22/22] revert --- tests/testthat/test-check_prior.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-check_prior.R b/tests/testthat/test-check_prior.R index 4b2512e8c..50028c4e2 100644 --- a/tests/testthat/test-check_prior.R +++ b/tests/testthat/test-check_prior.R @@ -92,7 +92,7 @@ test_that("check_prior - brms (linux)", { expect_warning(expect_identical( check_prior(model2, method = "lakeland")$Prior_Quality, c( - "informative", "informative", "informative", "informative", + "informative", "misinformative", "informative", "informative", "informative", "not determinable", "not determinable", "not determinable" ) ))