diff --git a/R/bayesfactor_models.R b/R/bayesfactor_models.R index db3bef13b..c8340868e 100644 --- a/R/bayesfactor_models.R +++ b/R/bayesfactor_models.R @@ -611,7 +611,7 @@ as.matrix.bayesfactor_models <- function(x, ...) { # Else... Get marginal likelihood if (verbose) { - message("Computation of Marginal Likelihood: estimating marginal likelihood, please wait...") + insight::format_alert("Computation of Marginal Likelihood: estimating marginal likelihood, please wait...") } # Should probably allow additional arguments such as reps or cores to for bridge_sampler bridgesampling::bridge_sampler(mod, silent = TRUE) diff --git a/R/bayesfactor_parameters.R b/R/bayesfactor_parameters.R index 93b993b29..e13f4367c 100644 --- a/R/bayesfactor_parameters.R +++ b/R/bayesfactor_parameters.R @@ -199,7 +199,7 @@ bayesfactor_pointnull <- function(posterior, verbose = TRUE, ...) { if (length(null) > 1 && verbose) { - message("'null' is a range - computing a ROPE based Bayes factor.") + insight::format_alert("`null` is a range - computing a ROPE based Bayes factor.") } bayesfactor_parameters( diff --git a/R/describe_posterior.R b/R/describe_posterior.R index 910e43421..063949ac1 100644 --- a/R/describe_posterior.R +++ b/R/describe_posterior.R @@ -499,8 +499,6 @@ describe_posterior.default <- function(posteriors, ...) { row.names(out) <- NULL } - - # Prepare output attr(out, "ci_method") <- ci_method out diff --git a/R/map_estimate.R b/R/map_estimate.R index 5c8da0be6..ac06708ec 100644 --- a/R/map_estimate.R +++ b/R/map_estimate.R @@ -1,16 +1,22 @@ #' Maximum A Posteriori probability estimate (MAP) #' -#' Find the **Highest Maximum A Posteriori probability estimate (MAP)** of a posterior, i.e., the value associated with the highest probability density (the "peak" of the posterior distribution). In other words, it is an estimation of the *mode* for continuous parameters. Note that this function relies on [estimate_density], which by default uses a different smoothing bandwidth (`"SJ"`) compared to the legacy default implemented the base R [density] function (`"nrd0"`). +#' Find the **Highest Maximum A Posteriori probability estimate (MAP)** of a +#' posterior, i.e., the value associated with the highest probability density +#' (the "peak" of the posterior distribution). In other words, it is an estimation +#' of the *mode* for continuous parameters. Note that this function relies on +#' [`estimate_density()`], which by default uses a different smoothing bandwidth +#' (`"SJ"`) compared to the legacy default implemented the base R [`density()`] +#' function (`"nrd0"`). #' #' @inheritParams hdi #' @inheritParams estimate_density #' #' @return A numeric value if `x` is a vector. If `x` is a model-object, #' returns a data frame with following columns: -#' \itemize{ -#' \item `Parameter` The model parameter(s), if `x` is a model-object. If `x` is a vector, this column is missing. -#' \item `MAP_Estimate` The MAP estimate for the posterior or each model parameter. -#' } +#' +#' - `Parameter`: The model parameter(s), if `x` is a model-object. If `x` is a +#' vector, this column is missing. +#' - `MAP_Estimate`: The MAP estimate for the posterior or each model parameter. #' #' @examples #' \dontrun{ diff --git a/R/mediation.R b/R/mediation.R index 6efd87d84..2f3dcd9eb 100644 --- a/R/mediation.R +++ b/R/mediation.R @@ -179,7 +179,7 @@ mediation.stanmvreg <- function(model, treatment, mediator, response = NULL, cen # check for binary response. In this case, user should rescale variables modelinfo <- insight::model_info(model) if (any(sapply(modelinfo, function(i) i$is_binomial, simplify = TRUE))) { - message("One of moderator or outcome is binary, so direct and indirect effects may be on different scales. Consider rescaling model predictors, e.g. with `effectsize::standardize()`.") + insight::format_alert("One of moderator or outcome is binary, so direct and indirect effects may be on different scales. Consider rescaling model predictors, e.g. with `effectsize::standardize()`.") } # model responses @@ -355,7 +355,7 @@ print.bayestestR_mediation <- function(x, digits = 3, ...) { ) if (any(prop_mediated_ori$Estimate < 0)) { - message("\nDirect and indirect effects have opposite directions. The proportion mediated is not meaningful.") + insight::format_alert("\nDirect and indirect effects have opposite directions. The proportion mediated is not meaningful.") } } diff --git a/R/p_direction.R b/R/p_direction.R index 310b50516..eeb1952d5 100644 --- a/R/p_direction.R +++ b/R/p_direction.R @@ -459,7 +459,7 @@ p_direction.get_predicted <- function(x, ...) { if ("iterations" %in% names(attributes(x))) { out <- p_direction(as.data.frame(t(attributes(x)$iterations)), ...) } else { - stop("No iterations present in the output.", call. = FALSE) + insight::format_error("No iterations present in the output.") } attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x)) out diff --git a/R/plot.R b/R/plot.R index bb1de42f5..03177a874 100644 --- a/R/plot.R +++ b/R/plot.R @@ -46,42 +46,49 @@ plot.bayestestR_eti <- function(x, ...) { NextMethod() } + #' @export plot.bayestestR_si <- function(x, ...) { insight::check_if_installed("see", "to plot support intervals") NextMethod() } + #' @export plot.bayesfactor_parameters <- function(x, ...) { insight::check_if_installed("see", "to plot Savage-Dickey Bayes factor") NextMethod() } + #' @export plot.bayesfactor_models <- function(x, ...) { insight::check_if_installed("see", "to plot models' Bayes factors") NextMethod() } + #' @export plot.estimate_density <- function(x, ...) { insight::check_if_installed("see", "to plot densities") NextMethod() } + #' @export plot.estimate_density_df <- function(x, ...) { insight::check_if_installed("see", "to plot models' densities") NextMethod() } + #' @export plot.p_significance <- function(x, ...) { insight::check_if_installed("see", "to plot practical significance") NextMethod() } + #' @export plot.describe_posterior <- function(x, stack = FALSE, ...) { insight::check_if_installed("see", "to plot posterior samples") @@ -91,6 +98,6 @@ plot.describe_posterior <- function(x, stack = FALSE, ...) { graphics::plot(estimate_density(model), stack = stack, ...) + ggplot2::labs(title = "Posterior Samples", x = NULL, y = NULL) } else { - warning(insight::format_message("Could not find model-object. Try ' plot(estimate_density(model))' instead."), call. = FALSE) + insight::format_warning("Could not find model-object. Try `plot(estimate_density(model))` instead.") } } diff --git a/R/point_estimate.R b/R/point_estimate.R index f23c14d6d..1b6dcf5ef 100644 --- a/R/point_estimate.R +++ b/R/point_estimate.R @@ -2,13 +2,22 @@ #' #' Compute various point-estimates, such as the mean, the median or the MAP, to describe posterior distributions. #' -#' @param centrality The point-estimates (centrality indices) to compute. Character (vector) or list with one or more of these options: `"median"`, `"mean"`, `"MAP"` or `"all"`. -#' @param dispersion Logical, if `TRUE`, computes indices of dispersion related to the estimate(s) (`SD` and `MAD` for `mean` and `median`, respectively). -#' @param threshold For `centrality = "trimmed"` (i.e. trimmed mean), indicates the fraction (0 to 0.5) of observations to be trimmed from each end of the vector before the mean is computed. +#' @param centrality The point-estimates (centrality indices) to compute. Character +#' (vector) or list with one or more of these options: `"median"`, `"mean"`, `"MAP"` +#' (see [`map_estimate()`]), `"trimmed"` (which is just `mean(x, trim = threshold)`), +#' `"mode"` or `"all"`. +#' @param dispersion Logical, if `TRUE`, computes indices of dispersion related +#' to the estimate(s) (`SD` and `MAD` for `mean` and `median`, respectively). +#' Dispersion is not available for `"MAP"` or `"mode"` centrality indices. +#' @param threshold For `centrality = "trimmed"` (i.e. trimmed mean), indicates +#' the fraction (0 to 0.5) of observations to be trimmed from each end of the +#' vector before the mean is computed. #' @param ... Additional arguments to be passed to or from methods. #' @inheritParams hdi #' -#' @references Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. (2019). *Indices of Effect Existence and Significance in the Bayesian Framework*. Frontiers in Psychology 2019;10:2767. \doi{10.3389/fpsyg.2019.02767} +#' @references Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. +#' (2019). *Indices of Effect Existence and Significance in the Bayesian Framework*. +#' Frontiers in Psychology 2019;10:2767. \doi{10.3389/fpsyg.2019.02767} #' #' @note There is also a [`plot()`-method](https://easystats.github.io/see/articles/bayestestR.html) implemented in the \href{https://easystats.github.io/see/}{\pkg{see}-package}. #' @@ -59,13 +68,15 @@ point_estimate <- function(x, ...) { #' @export point_estimate.default <- function(x, ...) { - stop(insight::format_message(paste0("'point_estimate()' is not yet implemented for objects of class '", class(x)[1], "'.")), call. = FALSE) + insight::format_error( + paste0("'point_estimate()' is not yet implemented for objects of class '", class(x)[1], "'.") + ) } #' @rdname point_estimate #' @export -point_estimate.numeric <- function(x, centrality = "all", dispersion = FALSE, threshold = .1, ...) { +point_estimate.numeric <- function(x, centrality = "all", dispersion = FALSE, threshold = 0.1, ...) { centrality <- match.arg(tolower(centrality), c("median", "mean", "map", "trimmed", "mode", "all"), several.ok = TRUE) if ("all" %in% centrality) { estimate_list <- c("median", "mean", "map") @@ -119,15 +130,8 @@ point_estimate.numeric <- function(x, centrality = "all", dispersion = FALSE, th } -.mode_estimate <- function(x) { - ux <- unique(x) - ux[which.max(tabulate(match(x, ux)))] -} - - - #' @export -point_estimate.data.frame <- function(x, centrality = "all", dispersion = FALSE, threshold = .1, ...) { +point_estimate.data.frame <- function(x, centrality = "all", dispersion = FALSE, threshold = 0.1, ...) { x <- .select_nums(x) if (ncol(x) == 1) { @@ -148,8 +152,14 @@ point_estimate.data.frame <- function(x, centrality = "all", dispersion = FALSE, #' @export -point_estimate.draws <- function(x, centrality = "all", dispersion = FALSE, threshold = .1, ...) { - point_estimate(.posterior_draws_to_df(x), centrality = centrality, dispersion = dispersion, threshold = threshold, ...) +point_estimate.draws <- function(x, centrality = "all", dispersion = FALSE, threshold = 0.1, ...) { + point_estimate( + .posterior_draws_to_df(x), + centrality = centrality, + dispersion = dispersion, + threshold = threshold, + ... + ) } #' @export @@ -186,16 +196,25 @@ point_estimate.BGGM <- point_estimate.bcplm #' @export point_estimate.bamlss <- function(x, centrality = "all", dispersion = FALSE, component = c("conditional", "location", "all"), ...) { component <- match.arg(component) - out <- point_estimate(insight::get_parameters(x, component = component), centrality = centrality, dispersion = dispersion, ...) - out <- .add_clean_parameters_attribute(out, x) - out + out <- point_estimate( + insight::get_parameters(x, component = component), + centrality = centrality, + dispersion = dispersion, + ... + ) + .add_clean_parameters_attribute(out, x) } #' @export point_estimate.MCMCglmm <- function(x, centrality = "all", dispersion = FALSE, ...) { nF <- x$Fixed$nfl - point_estimate(as.data.frame(x$Sol[, 1:nF, drop = FALSE]), centrality = centrality, dispersion = dispersion, ...) + point_estimate( + as.data.frame(x$Sol[, 1:nF, drop = FALSE]), + centrality = centrality, + dispersion = dispersion, + ... + ) } @@ -212,17 +231,6 @@ point_estimate.emmGrid <- function(x, centrality = "all", dispersion = FALSE, .. point_estimate.emm_list <- point_estimate.emmGrid -# Helper ------------------------------------------------------------------ - - - -#' @keywords internal -.point_estimate_models <- function(x, effects, component, parameters, centrality = "all", dispersion = FALSE, ...) { - out <- point_estimate(insight::get_parameters(x, effects = effects, component = component, parameters = parameters), centrality = centrality, dispersion = dispersion, ...) - out -} - - #' @rdname point_estimate #' @export point_estimate.stanreg <- function(x, centrality = "all", dispersion = FALSE, effects = c("fixed", "random", "all"), component = c("location", "all", "conditional", "smooth_terms", "sigma", "distributional", "auxiliary"), parameters = NULL, ...) { @@ -294,7 +302,6 @@ point_estimate.sim.merMod <- function(x, centrality = "all", dispersion = FALSE, } - #' @export point_estimate.sim <- function(x, centrality = "all", dispersion = FALSE, parameters = NULL, ...) { out <- .point_estimate_models( @@ -314,7 +321,6 @@ point_estimate.sim <- function(x, centrality = "all", dispersion = FALSE, parame } - #' @rdname point_estimate #' @export point_estimate.BFBayesFactor <- function(x, centrality = "all", dispersion = FALSE, ...) { @@ -341,3 +347,22 @@ point_estimate.get_predicted <- function(x, ...) { as.numeric(x) } } + + +# Helper ------------------------------------------------------------------ + +#' @keywords internal +.point_estimate_models <- function(x, effects, component, parameters, centrality = "all", dispersion = FALSE, ...) { + point_estimate( + insight::get_parameters(x, effects = effects, component = component, parameters = parameters), + centrality = centrality, + dispersion = dispersion, + ... + ) +} + +#' @keywords internal +.mode_estimate <- function(x) { + ux <- unique(x) + ux[which.max(tabulate(match(x, ux)))] +} diff --git a/R/reshape_iterations.R b/R/reshape_iterations.R index 1a07632df..ad7b13f2f 100644 --- a/R/reshape_iterations.R +++ b/R/reshape_iterations.R @@ -32,9 +32,9 @@ reshape_iterations <- function(x, prefix = c("draw", "iter", "iteration", "sim") prefix <- prefix[min(which(sapply(tolower(prefix), function(prefix) sum(grepl(prefix, tolower(names(x)), fixed = TRUE)) > 1)))] if (is.na(prefix) || is.null(prefix)) { - stop(insight::format_message( + insight::format_error( "Couldn't find columns corresponding to iterations in your dataframe, please specify the correct prefix." - ), call. = FALSE) + ) } # Get column names @@ -58,7 +58,7 @@ reshape_iterations <- function(x, prefix = c("draw", "iter", "iteration", "sim") ) row.names(long) <- NULL - class(long) <- class(long)[which(class(long) == "data.frame"):length(class(long))] + class(long) <- class(long)[which(inherits(long, "data.frame")):length(class(long))] long } diff --git a/R/rope.R b/R/rope.R index 146171c63..6072febd5 100644 --- a/R/rope.R +++ b/R/rope.R @@ -139,7 +139,7 @@ rope.numeric <- function(x, range = "default", ci = 0.95, ci_method = "ETI", ver if (all(range == "default")) { range <- c(-0.1, 0.1) } else if (!all(is.numeric(range)) || length(range) != 2) { - stop("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).", call. = FALSE) + insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") } rope_values <- lapply(ci, function(i) { @@ -231,7 +231,14 @@ rope.bamlss <- rope.BFBayesFactor #' @export rope.MCMCglmm <- function(x, range = "default", ci = 0.95, ci_method = "ETI", verbose = TRUE, ...) { nF <- x$Fixed$nfl - out <- rope(as.data.frame(x$Sol[, 1:nF, drop = FALSE]), range = range, ci = ci, ci_method = ci_method, verbose = verbose, ...) + out <- rope( + as.data.frame(x$Sol[, 1:nF, drop = FALSE]), + range = range, + ci = ci, + ci_method = ci_method, + verbose = verbose, + ... + ) attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x)) out } @@ -361,10 +368,14 @@ rope.brmsfit <- function(x, length(range) < length(insight::find_response(x)) || !all(names(range) %in% insight::find_response(x)) ) { - stop("With a multivariate model, `range` should be 'default' or a list of named numeric vectors with length 2.", call. = FALSE) + insight::format_error( + "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) { - stop("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).", call. = FALSE) + insight::format_error( + "`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1))." + ) } # check for possible collinearity that might bias ROPE and print a warning @@ -422,13 +433,20 @@ rope.brmsfit <- function(x, #' @export -rope.sim.merMod <- function(x, range = "default", ci = 0.95, ci_method = "ETI", effects = c("fixed", "random", "all"), parameters = NULL, verbose = TRUE, ...) { +rope.sim.merMod <- function(x, + range = "default", + ci = 0.95, + ci_method = "ETI", + effects = c("fixed", "random", "all"), + parameters = NULL, + verbose = TRUE, + ...) { effects <- match.arg(effects) if (all(range == "default")) { range <- rope_range(x, verbose = verbose) } else if (!all(is.numeric(range)) || length(range) != 2) { - stop("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).", call. = FALSE) + insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") } list <- lapply(c("fixed", "random"), function(.x) { @@ -489,7 +507,7 @@ rope.sim <- function(x, range = "default", ci = 0.95, ci_method = "ETI", paramet if (all(range == "default")) { range <- rope_range(x, verbose = verbose) } else if (!all(is.numeric(range)) || length(range) != 2) { - stop("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).", call. = FALSE) + insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).") } parms <- insight::get_parameters(x, parameters = parameters) diff --git a/R/simulate_data.R b/R/simulate_data.R index b59fb67c6..ba14aaf12 100644 --- a/R/simulate_data.R +++ b/R/simulate_data.R @@ -68,21 +68,21 @@ simulate_correlation <- function(n = 100, if (is.matrix(r)) { if (isSymmetric(r)) { if (any(r > 1)) { - stop("'r' should only contain values between -1 and 1.", call. = FALSE) + insight::format_error("`r` should only contain values between -1 and 1.") } else { sigma <- r } } else { - stop("'r' should be a symetric matrix (relative to the diagonal).", call. = FALSE) + insight::format_error("`r` should be a symetric matrix (relative to the diagonal).") } } else if (length(r) == 1L) { if (abs(r) > 1) { - stop("'r' should only contain values between -1 and 1.", call. = FALSE) + insight::format_error("`r` should only contain values between -1 and 1.") } else { sigma <- matrix(c(1, r, r, 1), nrow = 2) } } else { - stop("'r' should be a value (e.g., r = 0.5) or a square matrix.", call. = FALSE) + insight::format_error("`r` should be a value (e.g., r = 0.5) or a square matrix.") } @@ -107,10 +107,8 @@ simulate_correlation <- function(n = 100, data <- as.data.frame(data) # Rename - if (!is.null(names)) { - if (length(names) == ncol(data)) { - names(data) <- names - } + if (!is.null(names) && length(names) == ncol(data)) { + names(data) <- names } data } @@ -128,10 +126,8 @@ simulate_ttest <- function(n = 100, d = 0.5, names = NULL, ...) { data <- data.frame(y = as.factor(y), x = x) names(data) <- paste0("V", 0:(ncol(data) - 1)) - if (!is.null(names)) { - if (length(names) == ncol(data)) { - names(data) <- names - } + if (!is.null(names) && length(names) == ncol(data)) { + names(data) <- names } data } @@ -149,10 +145,8 @@ simulate_difference <- function(n = 100, d = 0.5, names = NULL, ...) { ) names(data) <- paste0("V", 0:(ncol(data) - 1)) - if (!is.null(names)) { - if (length(names) == ncol(data)) { - names(data) <- names - } + if (!is.null(names) && length(names) == ncol(data)) { + names(data) <- names } data } diff --git a/R/simulate_priors.R b/R/simulate_priors.R index 70ff0bf1a..1bb5780eb 100644 --- a/R/simulate_priors.R +++ b/R/simulate_priors.R @@ -127,7 +127,7 @@ simulate_prior.bcplm <- function(model, n = 1000, verbose = TRUE, ...) { } if (sim_error_msg && verbose) { - warning(paste0("Can't simulate priors from a ", prior$Distribution, " distribution."), call. = FALSE) + insight::format_warning(paste0("Can't simulate priors from a ", prior$Distribution, " distribution.")) } simulated$.bamboozled <- NULL diff --git a/R/simulate_simpson.R b/R/simulate_simpson.R index 01dc59ba7..63861a4a5 100644 --- a/R/simulate_simpson.R +++ b/R/simulate_simpson.R @@ -30,7 +30,7 @@ simulate_simpson <- function(n = 100, difference = 1, group_prefix = "G_") { if (n <= 3) { - stop("The number of observation `n` should be higher than 3", call. = FALSE) + insight::format_error("The number of observations `n` should be larger than 3.") } data <- data.frame() diff --git a/R/spi.R b/R/spi.R index 6a31fac5e..e5196c6da 100644 --- a/R/spi.R +++ b/R/spi.R @@ -18,8 +18,7 @@ #' than the HDI, because, the _"HDI can be noisy (that is, have a high Monte Carlo error)"_ #' (Liu et al. 2015). Furthermore, the HDI is sensitive to additional assumptions, #' in particular assumptions related to the different estimation methods, which -#' can make the HDI less accurate or reliable (see also discussion -#' [here](https://twitter.com/betanalpha/status/1479107186030624771)). +#' can make the HDI less accurate or reliable. #' #' @references #' Liu, Y., Gelman, A., & Zheng, T. (2015). Simulation-efficient shortest probability intervals. Statistics and Computing, 25(4), 809–819. https://doi.org/10.1007/s11222-015-9563-8 @@ -48,7 +47,7 @@ spi <- function(x, ...) { #' @export spi.default <- function(x, ...) { - stop(insight::format_message(paste0("'spi()' is not yet implemented for objects of class '", class(x)[1], "'.")), call. = FALSE) + insight::format_error(paste0("'spi()' is not yet implemented for objects of class '", class(x)[1], "'.")) } @@ -227,7 +226,7 @@ spi.get_predicted <- function(x, ...) { if ("iterations" %in% names(attributes(x))) { out <- spi(as.data.frame(t(attributes(x)$iterations)), ...) } else { - stop("No iterations present in the output.", call. = FALSE) + insight::format_error("No iterations present in the output.") } attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x)) out @@ -272,7 +271,7 @@ spi.get_predicted <- function(x, ...) { # lower bound - if (all(!is.na(k)) && all(k == 1)) { + if (!anyNA(k) && all(k == 1)) { x.l <- l } else { x.l <- tryCatch( @@ -284,7 +283,7 @@ spi.get_predicted <- function(x, ...) { frac <- 1 while (is.null(x.l)) { - frac <- frac - .1 + frac <- frac - 0.1 x.l <- tryCatch( { .spi_lower(bw = frac * bw, n.sims = n.sims, k = k, l = l, dens = dens, x = x) @@ -292,17 +291,15 @@ spi.get_predicted <- function(x, ...) { error = function(e) NULL ) - if (frac <= .1) { - message(insight::color_text(insight::format_message( - "Could not find a solution for the SPI lower bound." - ), color = "red")) + if (frac <= 0.1) { + insight::format_alert("Could not find a solution for the SPI lower bound.") x.l <- NA } } } # upper bound - if (all(!is.na(ui)) && all(ui == n.sims)) { + if (!anyNA(ui) && all(ui == n.sims)) { x.u <- u } else { x.u <- tryCatch( @@ -314,7 +311,7 @@ spi.get_predicted <- function(x, ...) { frac <- 1 while (is.null(x.u)) { - frac <- frac - .1 + frac <- frac - 0.1 x.u <- tryCatch( { .spi_upper(bw = frac * bw, n.sims = n.sims, ui = ui, u = u, dens = dens, x = x) @@ -322,10 +319,8 @@ spi.get_predicted <- function(x, ...) { error = function(e) NULL ) - if (frac <= .1) { - message(insight::color_text(insight::format_message( - "Could not find a solution for the SPI upper bound." - ), color = "red")) + if (frac <= 0.1) { + insight::format_alert("Could not find a solution for the SPI upper bound.") x.u <- NA } } @@ -370,35 +365,33 @@ spi.get_predicted <- function(x, ...) { # create constraint matrix A.l <- matrix(0, nrow = range_ll_lu + 3, ncol = range_ll_lu + 1) A.l[1, ] <- 1 - if (bw > 1) { - if (k > 2) { - for (j in 1:(range_ll_k - 1)) { - if (x[l.l + j + 1] == x[l.l + j]) { - A.l[1 + j, j + 1] <- 1 - A.l[1 + j, j + 2] <- -1 - } else { - aa <- (x[l.l + j] - x[l.l + j - 1]) / (x[l.l + j + 1] - x[l.l + j]) - A.l[1 + j, j] <- 1 - A.l[1 + j, j + 1] <- -(aa + 1) - A.l[1 + j, j + 2] <- aa - } + if (bw > 1 && k > 2) { + for (j in 1:(range_ll_k - 1)) { + if (x[l.l + j + 1] == x[l.l + j]) { + A.l[1 + j, j + 1] <- 1 + A.l[1 + j, j + 2] <- -1 + } else { + aa <- (x[l.l + j] - x[l.l + j - 1]) / (x[l.l + j + 1] - x[l.l + j]) + A.l[1 + j, j] <- 1 + A.l[1 + j, j + 1] <- -(aa + 1) + A.l[1 + j, j + 2] <- aa } + } - for (j in 0:(l.u - k - 2)) { - if (x[k + j + 1] == x[k + j + 2]) { - A.l[range_ll_k + 1 + j, range_ll_k + 2 + j] <- 1 - A.l[range_ll_k + 1 + j, range_ll_k + 3 + j] <- -1 - } else { - aa <- (x[k + j] - x[k + j + 1]) / (x[k + j + 1] - x[k + j + 2]) - A.l[range_ll_k + 1 + j, range_ll_k + 1 + j] <- -1 - A.l[range_ll_k + 1 + j, range_ll_k + 2 + j] <- aa + 1 - A.l[range_ll_k + 1 + j, range_ll_k + 3 + j] <- -aa - } + for (j in 0:(l.u - k - 2)) { + if (x[k + j + 1] == x[k + j + 2]) { + A.l[range_ll_k + 1 + j, range_ll_k + 2 + j] <- 1 + A.l[range_ll_k + 1 + j, range_ll_k + 3 + j] <- -1 + } else { + aa <- (x[k + j] - x[k + j + 1]) / (x[k + j + 1] - x[k + j + 2]) + A.l[range_ll_k + 1 + j, range_ll_k + 1 + j] <- -1 + A.l[range_ll_k + 1 + j, range_ll_k + 2 + j] <- aa + 1 + A.l[range_ll_k + 1 + j, range_ll_k + 3 + j] <- -aa } } } if (x[k + 1] == x[k]) { - aa <- (x[k] - x[k - 1]) / (x[k + 1] - x[k] + .000001) + aa <- (x[k] - x[k - 1]) / (x[k + 1] - x[k] + 0.000001) } else { aa <- (x[k] - x[k - 1]) / (x[k + 1] - x[k]) } @@ -449,33 +442,31 @@ spi.get_predicted <- function(x, ...) { # create constraint matrix A.u <- matrix(0, nrow = range_ul_uu + 3, ncol = range_ul_uu + 1) A.u[1, ] <- 1 - if (bw > 1) { - if (range_ul_ui > 1) { - for (j in 1:(range_ul_ui - 1)) { - if (x[u.l + j + 1] == x[u.l + j]) { - A.u[1 + j, j + 1] <- 1 - A.u[1 + j, j + 2] <- -1 - } else { - aa <- (x[u.l + j] - x[u.l + j - 1]) / (x[u.l + j + 1] - x[u.l + j]) - A.u[1 + j, j] <- 1 - A.u[1 + j, j + 1] <- -(aa + 1) - A.u[1 + j, j + 2] <- aa - } + if (bw > 1 && range_ul_ui > 1) { + for (j in 1:(range_ul_ui - 1)) { + if (x[u.l + j + 1] == x[u.l + j]) { + A.u[1 + j, j + 1] <- 1 + A.u[1 + j, j + 2] <- -1 + } else { + aa <- (x[u.l + j] - x[u.l + j - 1]) / (x[u.l + j + 1] - x[u.l + j]) + A.u[1 + j, j] <- 1 + A.u[1 + j, j + 1] <- -(aa + 1) + A.u[1 + j, j + 2] <- aa } + } - i <- 0 - for (j in (range_ul_ui):(range_ul_uu - 2)) { - if (x[ui + i + 1] == x[ui + i + 2]) { - A.u[1 + j, j + 2] <- 1 - A.u[1 + j, j + 3] <- -1 - } else { - aa <- (x[ui + i] - x[ui + i + 1]) / (x[ui + i + 1] - x[ui + i + 2]) - A.u[1 + j, j + 1] <- -1 - A.u[1 + j, j + 2] <- aa + 1 - A.u[1 + j, j + 3] <- -aa - } - i <- i + 1 + i <- 0 + for (j in (range_ul_ui):(range_ul_uu - 2)) { + if (x[ui + i + 1] == x[ui + i + 2]) { + A.u[1 + j, j + 2] <- 1 + A.u[1 + j, j + 3] <- -1 + } else { + aa <- (x[ui + i] - x[ui + i + 1]) / (x[ui + i + 1] - x[ui + i + 2]) + A.u[1 + j, j + 1] <- -1 + A.u[1 + j, j + 2] <- aa + 1 + A.u[1 + j, j + 3] <- -aa } + i <- i + 1 } } if (x[ui + 1] == x[ui]) { diff --git a/R/unupdate.R b/R/unupdate.R index 69ce6b351..0e0adacf3 100644 --- a/R/unupdate.R +++ b/R/unupdate.R @@ -9,7 +9,8 @@ #' #' @param model A fitted Bayesian model. #' @param verbose Toggle warnings. -#' @param newdata List of `data.frames` to update the model with new data. Required even if the original data should be used. +#' @param newdata List of `data.frames` to update the model with new data. +#' Required even if the original data should be used. #' @param ... Not used #' #' @return A model un-fitted to the data, representing the prior model. @@ -33,15 +34,13 @@ unupdate.stanreg <- function(model, verbose = TRUE, ...) { } if (verbose) { - message("Sampling priors, please wait...") + insight::format_alert("Sampling priors, please wait...") } prior_dists <- sapply(rstanarm::prior_summary(model), `[[`, "dist") if (anyNA(prior_dists)) { - stop( - "Cannot sample from flat priors (such as when priors are ", - "set to 'NULL' in a 'stanreg' model).", - call. = FALSE + insight::format_error( + "Cannot sample from flat priors (such as when priors are set to 'NULL' in a 'stanreg' model)." ) } @@ -64,24 +63,22 @@ unupdate.brmsfit <- function(model, verbose = TRUE, ...) { } if (verbose) { - message("Sampling priors, please wait...") + insight::format_alert("Sampling priors, please wait...") } - utils::capture.output( + utils::capture.output({ model_prior <- try(suppressMessages(suppressWarnings( stats::update(model, sample_prior = "only", refresh = 0) )), silent = TRUE) - ) + }) if (methods::is(model_prior, "try-error")) { if (grepl("proper priors", model_prior, fixed = TRUE)) { - stop( - "Cannot sample from flat priors (such as the default ", - "priors for fixed-effects in a 'brmsfit' model).", - call. = FALSE + insight::format_error( + "Cannot sample from flat priors (such as the default priors for fixed-effects in a 'brmsfit' model)." ) } else { - stop(model_prior, call. = FALSE) + insight::format_error(model_prior) } } @@ -102,10 +99,10 @@ unupdate.brmsfit_multiple <- function(model, } if (verbose) { - message("Sampling priors, please wait...") + insight::format_alert("Sampling priors, please wait...") } - utils::capture.output(model_prior <- + utils::capture.output({model_prior <- try(suppressMessages(suppressWarnings( stats::update( model, @@ -113,7 +110,7 @@ unupdate.brmsfit_multiple <- function(model, newdata = newdata, refresh = 0 ) - )), silent = TRUE)) + )), silent = TRUE)}) if (methods::is(model_prior, "try-error")) { if (grepl("proper priors", model_prior, fixed = TRUE)) { @@ -141,12 +138,12 @@ unupdate.blavaan <- function(model, verbose = TRUE, ...) { } if (verbose) { - message("Sampling priors, please wait...") + insight::format_alert("Sampling priors, please wait...") } cl$prisamp <- TRUE suppressMessages(suppressWarnings( - utils::capture.output(model_prior <- eval(cl)) + utils::capture.output({model_prior <- eval(cl)}) )) model_prior diff --git a/R/utils.R b/R/utils.R index 4ba6c590c..de396a834 100644 --- a/R/utils.R +++ b/R/utils.R @@ -46,7 +46,7 @@ direction <- Value[tolower(direction[1])] if (is.na(direction)) { - stop("Unrecognized 'direction' argument.", call. = FALSE) + insight::format_error("Unrecognized 'direction' argument.") } direction } @@ -61,24 +61,34 @@ return(temp) } if (isTRUE(is_stan_mv)) { + # for models with multiple responses, we create a separate response column temp$Response <- gsub("(b\\[)*(.*)\\|(.*)", "\\2", temp$Parameter) + # from the parameter names, we can now remove the name of the respone variables for (i in unique(temp$Response)) { temp$Parameter <- gsub(sprintf("%s|", i), "", temp$Parameter, fixed = TRUE) } merge_by <- c("Parameter", "Effects", "Component", "Response") remove_cols <- c("Group", "Cleaned_Parameter", "Function", ".roworder") } else if (isTRUE(is_brms_mv)) { + # for models with multiple responses, we create a separate response column temp$Response <- gsub("(.*)_(.*)_(.*)", "\\2", temp$Parameter) - # temp$Parameter <- gsub("(.*)_(.*)_(.*)", "\\1_\\3", temp$Parameter) merge_by <- c("Parameter", "Effects", "Component", "Response") remove_cols <- c("Group", "Cleaned_Parameter", "Function", ".roworder") } else { + # By default, we only merge by these three columns merge_by <- c("Parameter", "Effects", "Component") remove_cols <- c("Group", "Cleaned_Parameter", "Response", "Function", ".roworder") } + + # in "temp", we have the data frame from the related functions (like + # `point_estimate()`, `ci()` etc.). "cleaned_parameters" is a data frame + # only with original parameter names, model components and "cleaned" + # parameter names (retrieved from `insight::clean_parameters()`). + merge_by <- intersect(merge_by, colnames(temp)) temp$.roworder <- seq_len(nrow(temp)) out <- merge(x = temp, y = cleaned_parameters, by = merge_by, all.x = TRUE) + # hope this works for stanmvreg... if ((isTRUE(is_stan_mv) || isTRUE(is_brms_mv)) && all(is.na(out$Effects)) && all(is.na(out$Component))) { out$Effects <- cleaned_parameters$Effects[seq_len(nrow(out))] diff --git a/R/utils_bayesfactor.R b/R/utils_bayesfactor.R index fa49752d1..1ec903490 100644 --- a/R/utils_bayesfactor.R +++ b/R/utils_bayesfactor.R @@ -70,29 +70,25 @@ if (is.null(prior)) { prior <- posterior - warning( - "Prior not specified! ", - "Please provide the original model to get meaningful results." - ) + insight::format_warning("Prior not specified! Please provide the original model to get meaningful results.") } if (!inherits(prior, "emmGrid")) { # then is it a model on.exit( - stop( + insight::format_error(paste0( "Unable to reconstruct prior estimates.\n", "Perhaps the emmGrid object has been transformed or regrid()-ed?\n", "See function details.\n\n", "Instead, you can reestimate the emmGrid with a prior model, Try:\n", "\tprior_model <- unupdate(mode)\n", - "\tprior_emmgrid <- emmeans(prior_model, ...) # pass this as the 'prior' argument.", - call. = FALSE - ) + "\tprior_emmgrid <- emmeans(prior_model, ...) # pass this as the 'prior' argument." + )) ) if (inherits(prior, "brmsfit")) { - stop("Cannot rebuild prior emmGrid from a brmsfit model.", call. = FALSE) + insight::format_error("Cannot rebuild prior emmGrid from a brmsfit model.") } @@ -107,16 +103,14 @@ "See '?bayesfactor_parameters' for more information.\n" ) } - stop(prior, call. = FALSE) + insight::format_error(prior) } prior <- emmeans::ref_grid(prior) prior <- prior@post.beta if (!isTRUE(all.equal(colnames(prior), colnames(posterior@post.beta)))) { - stop("post.beta and prior.beta are non-conformable arguments.", - call. = FALSE - ) + insight::format_error("post.beta and prior.beta are non-conformable arguments.") } prior <- stats::update(posterior, post.beta = prior) on.exit() # undo general error message @@ -135,17 +129,14 @@ verbose = TRUE) { if (is.null(prior)) { prior <- posterior - warning( - "Prior not specified! ", - "Please provide the original model to get meaningful results." - ) + insight::format_warning("Prior not specified! Please provide the original model to get meaningful results.") } if (!inherits(prior, "emm_list")) { # prior is a model if (inherits(prior, "brmsfit")) { - stop("Cannot rebuild prior emm_list from a brmsfit model.", call. = FALSE) + insight::format_error("Cannot rebuild prior emm_list from a brmsfit model.") } prior <- try(unupdate(prior, verbose = verbose), silent = TRUE) @@ -158,7 +149,7 @@ "See '?bayesfactor_parameters' for more information.\n" ) } - stop(prior, call. = FALSE) + insight::format_error(prior) } } # prior is now a model, or emm_list diff --git a/R/utils_check_collinearity.R b/R/utils_check_collinearity.R index 69a50f14d..8fb8636d2 100644 --- a/R/utils_check_collinearity.R +++ b/R/utils_check_collinearity.R @@ -45,14 +45,26 @@ results <- results[results$corr > threshold & results$corr <= 0.9, ] if (nrow(results) > 0) { where <- paste0("between ", toString(paste0(results$where, " (r = ", round(results$corr, 2), ")")), "") - message("Possible multicollinearity ", where, ". This might lead to inappropriate results. See 'Details' in '?", method, "'.") + insight::format_alert(paste0( + "Possible multicollinearity ", + where, + ". This might lead to inappropriate results. See 'Details' in '?", + method, + "'." + )) } # Filter by second threshold results <- results[results$corr > 0.9, ] if (nrow(results) > 0) { where <- paste0("between ", toString(paste0(results$where, " (r = ", round(results$corr, 2), ")")), "") - warning("Probable multicollinearity ", where, ". This might lead to inappropriate results. See 'Details' in '?", method, "'.", call. = FALSE) + insight::format_warning(paste0( + "Probable multicollinearity ", + where, + ". This might lead to inappropriate results. See 'Details' in '?", + method, + "'." + )) } } } diff --git a/man/describe_posterior.Rd b/man/describe_posterior.Rd index e190a61c7..71681d95c 100644 --- a/man/describe_posterior.Rd +++ b/man/describe_posterior.Rd @@ -74,9 +74,14 @@ for other classes mostly resemble the arguments of the \code{.numeric} method.} \item{...}{Additional arguments to be passed to or from methods.} -\item{centrality}{The point-estimates (centrality indices) to compute. Character (vector) or list with one or more of these options: \code{"median"}, \code{"mean"}, \code{"MAP"} or \code{"all"}.} - -\item{dispersion}{Logical, if \code{TRUE}, computes indices of dispersion related to the estimate(s) (\code{SD} and \code{MAD} for \code{mean} and \code{median}, respectively).} +\item{centrality}{The point-estimates (centrality indices) to compute. Character +(vector) or list with one or more of these options: \code{"median"}, \code{"mean"}, \code{"MAP"} +(see \code{\link[=map_estimate]{map_estimate()}}), \code{"trimmed"} (which is just \code{mean(x, trim = threshold)}), +\code{"mode"} or \code{"all"}.} + +\item{dispersion}{Logical, if \code{TRUE}, computes indices of dispersion related +to the estimate(s) (\code{SD} and \code{MAD} for \code{mean} and \code{median}, respectively). +Dispersion is not available for \code{"MAP"} or \code{"mode"} centrality indices.} \item{ci}{Value or vector of probability of the CI (between 0 and 1) to be estimated. Default to \code{.95} (\verb{95\%}).} diff --git a/man/map_estimate.Rd b/man/map_estimate.Rd index 6d5842933..e63ad348b 100644 --- a/man/map_estimate.Rd +++ b/man/map_estimate.Rd @@ -65,12 +65,19 @@ for the output.} A numeric value if \code{x} is a vector. If \code{x} is a model-object, returns a data frame with following columns: \itemize{ -\item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a vector, this column is missing. -\item \code{MAP_Estimate} The MAP estimate for the posterior or each model parameter. +\item \code{Parameter}: The model parameter(s), if \code{x} is a model-object. If \code{x} is a +vector, this column is missing. +\item \code{MAP_Estimate}: The MAP estimate for the posterior or each model parameter. } } \description{ -Find the \strong{Highest Maximum A Posteriori probability estimate (MAP)} of a posterior, i.e., the value associated with the highest probability density (the "peak" of the posterior distribution). In other words, it is an estimation of the \emph{mode} for continuous parameters. Note that this function relies on \link{estimate_density}, which by default uses a different smoothing bandwidth (\code{"SJ"}) compared to the legacy default implemented the base R \link{density} function (\code{"nrd0"}). +Find the \strong{Highest Maximum A Posteriori probability estimate (MAP)} of a +posterior, i.e., the value associated with the highest probability density +(the "peak" of the posterior distribution). In other words, it is an estimation +of the \emph{mode} for continuous parameters. Note that this function relies on +\code{\link[=estimate_density]{estimate_density()}}, which by default uses a different smoothing bandwidth +(\code{"SJ"}) compared to the legacy default implemented the base R \code{\link[=density]{density()}} +function (\code{"nrd0"}). } \examples{ \dontrun{ diff --git a/man/mediation.Rd b/man/mediation.Rd index 5efd7e7c3..fd28f5a6b 100644 --- a/man/mediation.Rd +++ b/man/mediation.Rd @@ -57,7 +57,10 @@ second response variable (for the treatment model, with the mediator as additional predictor), \code{y}, is not transformed. Then we could use \code{response} like this: \code{mediation(model, response = c(m = "m_center", y = "y"))}.} -\item{centrality}{The point-estimates (centrality indices) to compute. Character (vector) or list with one or more of these options: \code{"median"}, \code{"mean"}, \code{"MAP"} or \code{"all"}.} +\item{centrality}{The point-estimates (centrality indices) to compute. Character +(vector) or list with one or more of these options: \code{"median"}, \code{"mean"}, \code{"MAP"} +(see \code{\link[=map_estimate]{map_estimate()}}), \code{"trimmed"} (which is just \code{mean(x, trim = threshold)}), +\code{"mode"} or \code{"all"}.} \item{ci}{Value or vector of probability of the CI (between 0 and 1) to be estimated. Default to \code{.95} (\verb{95\%}).} diff --git a/man/point_estimate.Rd b/man/point_estimate.Rd index 4a0727a5b..d4baa14c1 100644 --- a/man/point_estimate.Rd +++ b/man/point_estimate.Rd @@ -44,11 +44,18 @@ resemble the arguments of the \code{.numeric} or \code{.data.frame}methods.} \item{...}{Additional arguments to be passed to or from methods.} -\item{centrality}{The point-estimates (centrality indices) to compute. Character (vector) or list with one or more of these options: \code{"median"}, \code{"mean"}, \code{"MAP"} or \code{"all"}.} +\item{centrality}{The point-estimates (centrality indices) to compute. Character +(vector) or list with one or more of these options: \code{"median"}, \code{"mean"}, \code{"MAP"} +(see \code{\link[=map_estimate]{map_estimate()}}), \code{"trimmed"} (which is just \code{mean(x, trim = threshold)}), +\code{"mode"} or \code{"all"}.} -\item{dispersion}{Logical, if \code{TRUE}, computes indices of dispersion related to the estimate(s) (\code{SD} and \code{MAD} for \code{mean} and \code{median}, respectively).} +\item{dispersion}{Logical, if \code{TRUE}, computes indices of dispersion related +to the estimate(s) (\code{SD} and \code{MAD} for \code{mean} and \code{median}, respectively). +Dispersion is not available for \code{"MAP"} or \code{"mode"} centrality indices.} -\item{threshold}{For \code{centrality = "trimmed"} (i.e. trimmed mean), indicates the fraction (0 to 0.5) of observations to be trimmed from each end of the vector before the mean is computed.} +\item{threshold}{For \code{centrality = "trimmed"} (i.e. trimmed mean), indicates +the fraction (0 to 0.5) of observations to be trimmed from each end of the +vector before the mean is computed.} \item{effects}{Should results for fixed effects, random effects or both be returned? Only applies to mixed models. May be abbreviated.} @@ -110,5 +117,7 @@ point_estimate(bf, centrality = c("median", "MAP")) } \references{ -Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. (2019). \emph{Indices of Effect Existence and Significance in the Bayesian Framework}. Frontiers in Psychology 2019;10:2767. \doi{10.3389/fpsyg.2019.02767} +Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. +(2019). \emph{Indices of Effect Existence and Significance in the Bayesian Framework}. +Frontiers in Psychology 2019;10:2767. \doi{10.3389/fpsyg.2019.02767} } diff --git a/man/spi.Rd b/man/spi.Rd index d5f2adbc5..085be3a30 100644 --- a/man/spi.Rd +++ b/man/spi.Rd @@ -78,8 +78,7 @@ uncertainty of (posterior) distributions. The SPI is said to be more stable than the HDI, because, the \emph{"HDI can be noisy (that is, have a high Monte Carlo error)"} (Liu et al. 2015). Furthermore, the HDI is sensitive to additional assumptions, in particular assumptions related to the different estimation methods, which -can make the HDI less accurate or reliable (see also discussion -\href{https://twitter.com/betanalpha/status/1479107186030624771}{here}). +can make the HDI less accurate or reliable. } \note{ The code to compute the SPI was adapted from the \strong{SPIn} package, diff --git a/man/unupdate.Rd b/man/unupdate.Rd index 42f037555..2d38361c8 100644 --- a/man/unupdate.Rd +++ b/man/unupdate.Rd @@ -25,7 +25,8 @@ unupdate(model, verbose = TRUE, ...) \item{...}{Not used} -\item{newdata}{List of \code{data.frames} to update the model with new data. Required even if the original data should be used.} +\item{newdata}{List of \code{data.frames} to update the model with new data. +Required even if the original data should be used.} } \value{ A model un-fitted to the data, representing the prior model. diff --git a/vignettes/web_only/indicesEstimationComparison.Rmd b/vignettes/web_only/indicesEstimationComparison.Rmd index b1be59b35..4e8cad52c 100644 --- a/vignettes/web_only/indicesEstimationComparison.Rmd +++ b/vignettes/web_only/indicesEstimationComparison.Rmd @@ -198,7 +198,8 @@ dat <- reshape_longer( out <- glm(true_effect ~ outcome_type / estimate / value, data = dat, family = "binomial") out <- parameters(out, ci_method = "wald") out <- data_select(out, c("Parameter", "Coefficient", "p")) -out <- data_filter(out, grep("^outcome_type(.*):value$", x = out$Parameter)) +rows <- grep("^outcome_type(.*):value$", x = out$Parameter) +out <- data_filter(out, rows) out <- out[order(out$Coefficient, decreasing = TRUE), ] knitr::kable(out, digits = 2) ```