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/NEWS.md b/NEWS.md index db28c58f5..c3232190d 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()`, `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 - return results with appended grid-data. + return results with appended grid-data. * Usability improvements for `p_direction()`: 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/equivalence_test.R b/R/equivalence_test.R index 57d4b2058..638af6102 100644 --- a/R/equivalence_test.R +++ b/R/equivalence_test.R @@ -85,6 +85,10 @@ #' \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")) +#' # named ROPE ranges +#' equivalence_test(model, range = list(wt = c(-5, -4), `(Intercept)` = c(10, 40))) #' #' # plot result #' test <- equivalence_test(model) @@ -163,13 +167,33 @@ 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)) { + # check if list of values contains only valid values + range <- .check_list_range(range, x) + # apply thresholds to each column + 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( @@ -244,7 +268,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,24 +281,7 @@ 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 - ) - - 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/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] ) } diff --git a/R/p_rope.R b/R/p_rope.R index df8f3a9e4..5a93dda7b 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 46c8bd1c5..12a692c1b 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -15,7 +15,11 @@ #' (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 +#' - 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 #' @@ -53,6 +57,10 @@ #' chains = 2, refresh = 0 #' ) #' p_significance(model) +#' # multiple thresholds - asymmetric, symmetric, default +#' p_significance(model, threshold = list(c(-10, 5), 0.2, "default")) +#' # named thresholds +#' p_significance(model, threshold = list(wt = 0.2, `(Intercept)` = c(-10, 5))) #' } #' @export p_significance <- function(x, ...) { @@ -128,11 +136,26 @@ 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) { ps <- .p_significance(x[, 1], threshold = threshold, ...) + } else if (is.list(threshold)) { + # check if list of values contains only valid values + threshold <- .check_list_range(threshold, x, larger_two = TRUE) + # apply thresholds to each column + ps <- mapply( + function(p, thres) { + .p_significance( + p, + threshold = thres + ) + }, + x, + threshold, + SIMPLIFY = FALSE + ) } else { ps <- sapply(x, .p_significance, threshold = threshold, simplify = TRUE, ...) } @@ -268,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")) @@ -304,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) @@ -364,7 +393,42 @@ 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) { + out <- c(-1 * abs(out), abs(out)) + } + out + }) + } 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,13 +436,49 @@ 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)`)" ) } threshold } + +.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(named_range[is.na(pos)]) + )) + } + # now "fill" non-specified elements with "default" + out <- as.list(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.") + } + # 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)).") + } + range +} diff --git a/R/print.R b/R/print.R index fbe3a2fd6..dfc62eeb7 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,26 @@ 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$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_default( x = x, digits = digits, caption = caption, + ci_string = ci_string, ... ) } @@ -302,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.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/R/print_html.R b/R/print_html.R index a7a4be75c..fb487c405 100644 --- a/R/print_html.R +++ b/R/print_html.R @@ -31,22 +31,41 @@ 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", ...) } #' @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, ...) } @@ -86,7 +105,7 @@ print_html.bayesfactor_models <- function(x, log = log, show_names = show_names, caption = caption, - align = c("llr"), + align = "llr", ... ) } @@ -103,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 f8910a123..053ae0472 100644 --- a/R/print_md.R +++ b/R/print_md.R @@ -30,22 +30,41 @@ 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", ...) } #' @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, ...) } @@ -85,7 +104,7 @@ print_md.bayesfactor_models <- function(x, log = log, show_names = show_names, caption = caption, - align = c("llr"), + align = "llr", ... ) } @@ -102,7 +121,7 @@ print_md.bayesfactor_inclusion <- function(x, digits = digits, log = log, caption = caption, - align = c("lrrr"), + align = "lrrr", ... ) } diff --git a/R/rope.R b/R/rope.R index 782967a8f..24c3a9b4a 100644 --- a/R/rope.R +++ b/R/rope.R @@ -6,13 +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` 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. 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 @@ -33,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 @@ -110,6 +118,12 @@ #' 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")) +#' +#' # 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)) #' @@ -381,7 +395,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 +456,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 +528,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 +588,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 +621,35 @@ 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)) { + # check if list of values contains only valid values + range <- .check_list_range(range, parms) + # apply thresholds to each column + 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..27fc778f7 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} 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. 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.} @@ -164,6 +173,10 @@ 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, 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_rope.Rd b/man/p_rope.Rd index 8791c0224..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} 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. 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 30b2170e5..5d6d0fbaf 100644 --- a/man/p_significance.Rd +++ b/man/p_significance.Rd @@ -59,7 +59,11 @@ 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 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, @@ -131,6 +135,10 @@ model <- rstanarm::stan_glm(mpg ~ wt + cyl, chains = 2, refresh = 0 ) p_significance(model) +# multiple thresholds - asymmetric, symmetric, default +p_significance(model, threshold = list(c(-10, 5), 0.2, "default")) +# 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 87036d460..242edb880 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} 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. 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 @@ -171,6 +180,12 @@ 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")) + +# 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/_snaps/equivalence_test.md b/tests/testthat/_snaps/equivalence_test.md new file mode 100644 index 000000000..cf351f609 --- /dev/null +++ b/tests/testthat/_snaps/equivalence_test.md @@ -0,0 +1,72 @@ +# equivalence test, rstanarm + + 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] + + + +# 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-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), diff --git a/tests/testthat/test-describe_posterior.R b/tests/testthat/test-describe_posterior.R index f21a10787..db265a679 100644 --- a/tests/testthat/test-describe_posterior.R +++ b/tests/testthat/test-describe_posterior.R @@ -146,6 +146,20 @@ 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) + + 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-equivalence_test.R b/tests/testthat/test-equivalence_test.R new file mode 100644 index 000000000..e7ce74cac --- /dev/null +++ b/tests/testthat/test-equivalence_test.R @@ -0,0 +1,69 @@ +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") + + 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'" + ) +}) + + +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'" + ) +}) diff --git a/tests/testthat/test-p_rope.R b/tests/testthat/test-p_rope.R new file mode 100644 index 000000000..854722b76 --- /dev/null +++ b/tests/testthat/test-p_rope.R @@ -0,0 +1,19 @@ +test_that("p_rope", { + skip_if_offline() + skip_if_not_or_load_if_installed("rstanarm") + m <- insight::download_model("stanreg_merMod_5") + 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 6f542f3c0..d11ecbb25 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( @@ -49,7 +48,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 +57,50 @@ test_that("brms", { c(1.0000, 0.9985, 0.9785), tolerance = 0.01 ) + + out <- p_significance(m2, threshold = list(1, "default", 2), effects = "all") + expect_equal( + out$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 + ) + + 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" + ) +}) + +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 8367e3d80..0c9059854 100644 --- a/tests/testthat/test-rope.R +++ b/tests/testthat/test-rope.R @@ -49,6 +49,33 @@ test_that("rope", { rope(p, verbose = FALSE)$ROPE_Percentage, 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" + ) + expect_error( + 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" + ) })