diff --git a/NEWS.md b/NEWS.md index 87219d7cd..3bccaba79 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,12 @@ # bayestestR 0.13.2 +## Breaking Changes + +* `pd_to_p()` now returns 1 and a warning for pds smaller than 0.5. +* `map_estimate()`, `p_direction()`, `p_map()`, and `p_significance()` now + return a data-frame when the input is a numeric vector. (making the output + consistently a data frame for all inputs.) + ## Changes * Retrieving models from the environment was improved. diff --git a/R/convert_pd_to_p.R b/R/convert_pd_to_p.R index 50e538a81..e63e52f80 100644 --- a/R/convert_pd_to_p.R +++ b/R/convert_pd_to_p.R @@ -6,17 +6,52 @@ #' @param p A p-value. #' @param direction What type of p-value is requested or provided. Can be #' `"two-sided"` (default, two tailed) or `"one-sided"` (one tailed). +#' @param verbose Toggle off warnings. #' @param ... Arguments passed to or from other methods. #' +#' @details +#' Conversion is done using the following equation (see Makowski et al., 2019): +#' \cr\cr +#' When `direction = "two-sided"` - +#' \cr\cr +#' \deqn{p = 2 \times (1 - p_d)}{p = 2 * (1 - pd)} +#' When `direction = "one-sided"` - +#' \cr\cr +#' \deqn{p = 1 - p_d}{p = 1 - pd} +#' \cr\cr +#' Note that this conversion is only valid when the lowest possible values of pd +#' is 0.5 - i.e., when the posterior represents continuous parameter space (see +#' [p_direction]). If any pd < 0.5 are detected, they are converted to a p of 1, +#' and a warning is given. +#' +#' @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} +#' #' @examples #' pd_to_p(pd = 0.95) #' pd_to_p(pd = 0.95, direction = "one-sided") +#' #' @export -pd_to_p <- function(pd, direction = "two-sided", ...) { - p <- 1 - pmax(pd, 1 - pd) +pd_to_p <- function(pd, direction = "two-sided", verbose = TRUE, ...) { + p <- 1 - pd if (.get_direction(direction) == 0) { p <- 2 * p } + + less_than_0.5 <- pd < 0.5 + if (any(less_than_0.5)) { + if (verbose) { + insight::format_warning(paste( + "pd values smaller than 0.5 detected.", + "pd-to-p conversion assumes a continious parameter space;", + "see help('p_direction') for more info." + )) + } + p[less_than_0.5] <- 1 + } + p } diff --git a/R/map_estimate.R b/R/map_estimate.R index ac06708ec..c0ee9345c 100644 --- a/R/map_estimate.R +++ b/R/map_estimate.R @@ -49,21 +49,20 @@ map_estimate <- function(x, precision = 2^10, method = "kernel", ...) { #' @rdname map_estimate #' @export map_estimate.numeric <- function(x, precision = 2^10, method = "kernel", ...) { - d <- estimate_density(x, precision = precision, method = method, ...) - - hdp_x <- d$x[which.max(d$y)] - hdp_y <- max(d$y) - - out <- hdp_x - attr(out, "MAP_density") <- hdp_y - + out <- map_estimate(data.frame(x = x), + precision, method = method, ...) + out[[1]] <- NULL attr(out, "data") <- x - attr(out, "centrality") <- "map" - class(out) <- unique(c("map_estimate", "see_point_estimate", class(out))) - out } +.map_estimate <- function(x, precision = 2^10, method = "kernel", ...) { + d <- estimate_density(x, precision = precision, method = method, ...) + + out <- d$x[which.max(d$y)] + attr(out, "MAP_density") <- max(d$y) + out +} # other models ----------------------- @@ -98,7 +97,7 @@ map_estimate.mcmc.list <- map_estimate.bayesQR #' @keywords internal .map_estimate_models <- function(x, precision, method, ...) { - l <- sapply(x, map_estimate, precision = precision, method = method, simplify = FALSE, ...) + l <- sapply(x, .map_estimate, precision = precision, method = method, simplify = FALSE, ...) out <- data.frame( Parameter = colnames(x), diff --git a/R/p_direction.R b/R/p_direction.R index b743ce9d2..288e773c7 100644 --- a/R/p_direction.R +++ b/R/p_direction.R @@ -1,81 +1,103 @@ #' Probability of Direction (pd) #' -#' Compute the **Probability of Direction** (***pd***, also known -#' as the Maximum Probability of Effect - *MPE*). It varies between `50%` -#' and `100%` (*i.e.*, `0.5` and `1`) and can be interpreted as -#' the probability (expressed in percentage) that a parameter (described by its -#' posterior distribution) is strictly positive or negative (whichever is the -#' most probable). It is mathematically defined as the proportion of the -#' posterior distribution that is of the median's sign. Although differently -#' expressed, this index is fairly similar (*i.e.*, is strongly correlated) -#' to the frequentist **p-value**. -#' \cr\cr -#' Note that in some (rare) cases, especially when used with model averaged -#' posteriors (see [weighted_posteriors()] or -#' `brms::posterior_average`), `pd` can be smaller than `0.5`, -#' reflecting high credibility of `0`. +#' Compute the **Probability of Direction** (***pd***, also known as the Maximum +#' Probability of Effect - *MPE*). This can be interpreted as the probability +#' that a parameter (described by its posterior distribution) is strictly +#' positive or negative (whichever is the most probable). Although differently +#' expressed, this index is fairly similar (*i.e.*, is strongly correlated) to +#' the frequentist **p-value** (see details). #' -#' @param x Vector representing a posterior distribution. Can also be a Bayesian model (`stanreg`, `brmsfit` or `BayesFactor`). -#' @param method Can be `"direct"` or one of methods of [density estimation][estimate_density], such as `"kernel"`, `"logspline"` or `"KernSmooth"`. If `"direct"` (default), the computation is based on the raw ratio of samples superior and inferior to 0. Else, the result is based on the [Area under the Curve (AUC)][auc] of the estimated [density][estimate_density] function. -#' @param null The value considered as a "null" effect. Traditionally 0, but could also be 1 in the case of ratios. +#' @param x A vector representing a posterior distribution, a data frame of +#' posterior draws (samples be parameter). Can also be a Bayesian model. +#' @param method Can be `"direct"` or one of methods of [`estimate_density()`], +#' such as `"kernel"`, `"logspline"` or `"KernSmooth"`. See details. +#' @param null The value considered as a "null" effect. Traditionally 0, but +#' could also be 1 in the case of ratios of change (OR, IRR, ...). #' @inheritParams hdi #' #' @details -#' \subsection{What is the *pd*?}{ -#' The Probability of Direction (pd) is an index of effect existence, ranging -#' from `50%` to `100%`, representing the certainty with which an effect goes in -#' a particular direction (*i.e.*, is positive or negative). Beyond its -#' simplicity of interpretation, understanding and computation, this index also -#' presents other interesting properties: -#' \itemize{ -#' \item It is independent from the model: It is solely based on the posterior +#' ## What is the *pd*? +#' The Probability of Direction (pd) is an index of effect existence, representing +#' the certainty with which an effect goes in a particular direction (i.e., is +#' positive or negative / has a sign), typically ranging from 0.5 to 1 (but see +#' next section for cases where it can range between 0 and 1). Beyond +#' its simplicity of interpretation, understanding and computation, this index +#' also presents other interesting properties: +#' - Like other posterior-based indices, *pd* is solely based on the posterior #' distributions and does not require any additional information from the data -#' or the model. -#' \item It is robust to the scale of both the response variable and the predictors. -#' \item It is strongly correlated with the frequentist p-value, and can thus +#' or the model (e.g., such as priors, as in the case of Bayes factors). +#' - It is robust to the scale of both the response variable and the predictors. +#' - It is strongly correlated with the frequentist p-value, and can thus #' be used to draw parallels and give some reference to readers non-familiar -#' with Bayesian statistics. -#' } -#' } -#' \subsection{Relationship with the p-value}{ -#' In most cases, it seems that the *pd* has a direct correspondence with the frequentist one-sided *p*-value through the formula \ifelse{html}{\out{pone sided = 1 - p(d)/100}}{\eqn{p_{one sided}=1-\frac{p_{d}}{100}}} and to the two-sided p-value (the most commonly reported one) through the formula \ifelse{html}{\out{ptwo sided = 2 * (1 - p(d)/100)}}{\eqn{p_{two sided}=2*(1-\frac{p_{d}}{100})}}. Thus, a two-sided p-value of respectively `.1`, `.05`, `.01` and `.001` would correspond approximately to a *pd* of `95%`, `97.5%`, `99.5%` and `99.95%`. See also [pd_to_p()]. -#' } -#' \subsection{Methods of computation}{ -#' The most simple and direct way to compute the *pd* is to 1) look at the -#' median's sign, 2) select the portion of the posterior of the same sign and -#' 3) compute the percentage that this portion represents. This "simple" method -#' is the most straightforward, but its precision is directly tied to the -#' number of posterior draws. The second approach relies on [density -#' estimation][estimate_density]. It starts by estimating the density function -#' (for which many methods are available), and then computing the [area under -#' the curve][area_under_curve] (AUC) of the density curve on the other side of -#' 0. -#' } -#' \subsection{Strengths and Limitations}{ -#' **Strengths:** Straightforward computation and interpretation. Objective -#' property of the posterior distribution. 1:1 correspondence with the -#' frequentist p-value. -#' \cr \cr -#' **Limitations:** Limited information favoring the null hypothesis. -#' } +#' with Bayesian statistics (Makowski et al., 2019). #' -#' @return -#' Values between 0.5 and 1 corresponding to the probability of direction (pd). +#' ## Relationship with the p-value +#' In most cases, it seems that the *pd* has a direct correspondence with the +#' frequentist one-sided *p*-value through the formula (for two-sided *p*): +#' \deqn{p = 2 \times (1 - p_d)}{p = 2 * (1 - pd)} +#' Thus, a two-sided p-value of respectively `.1`, `.05`, `.01` and `.001` would +#' correspond approximately to a *pd* of `95%`, `97.5%`, `99.5%` and `99.95%`. +#' See [pd_to_p()] for details. +#' +#' ## Possible Range of Values +#' The largest value *pd* can take is 1 - the posterior is strictly directional. +#' However, the smallest value *pd* can take depends on the parameter space +#' represented by the posterior. +#' \cr\cr +#' **For a continuous parameter space**, exact values of 0 (or any point null +#' value) are not possible, and so 100% of the posterior has _some_ sign, some +#' positive, some negative. Therefore, the smallest the *pd* can be is 0.5 - +#' with an equal posterior mass of positive and negative values. Values close to +#' 0.5 _cannot_ be used to support the null hypothesis (that the parameter does +#' _not_ have a direction) is a similar why to how large p-values cannot be used +#' to support the null hypothesis (see [`pd_tp_p()`]; Makowski et al., 2019). #' \cr\cr -#' Note that in some (rare) cases, especially when used with model averaged -#' posteriors (see [weighted_posteriors()] or -#' `brms::posterior_average`), `pd` can be smaller than `0.5`, -#' reflecting high credibility of `0`. To detect such cases, the -#' `method = "direct"` must be used. +#' **For a discrete parameter space or a parameter space that is a mixture +#' between discrete and continuous spaces**, exact values of 0 (or any point +#' null value) _are_ possible! Therefore, the smallest the *pd* can be is 0 - +#' with 100% of the posterior mass on 0. Thus values close to 0 can be used to +#' support the null hypothesis (see van den Bergh et al., 2021). +#' \cr\cr +#' Examples of posteriors representing discrete parameter space: +#' - When a parameter can only take discrete values. +#' - When a mixture prior/posterior is used (such as the spike-and-slab prior; +#' see van den Bergh et al., 2021). +#' - When conducting Bayesian model averaging (e.g., [weighted_posteriors()] or +#' `brms::posterior_average`). +#' +#' ## Methods of computation +#' The *pd* is defined as: +#' \deqn{p_d = max({Pr(\hat{\theta} < \theta_{null}), Pr(\hat{\theta} > \theta_{null})})}{pd = max(mean(x < null), mean(x > null))} +#' \cr\cr +#' The most simple and direct way to compute the *pd* is to compute the +#' proportion of positive (or larger than `null`) posterior samples, the +#' proportion of negative (or smaller than `null`) posterior samples, and take +#' the larger of the two. This "simple" method is the most straightforward, but +#' its precision is directly tied to the number of posterior draws. +#' \cr\cr +#' The second approach relies on [`density estimation()`]: It starts by +#' estimating the continuous-smooth density function (for which many methods are +#' available), and then computing the [area under the curve][area_under_curve] +#' (AUC) of the density curve on either side of `null` and taking the maximum +#' between them. Note the this approach assumes a continuous density function, +#' and so **when the posterior represents a (partially) discrete parameter +#' space, only the direct method _must_ be used** (see above). +#' +#' @return +#' Values between 0.5 and 1 *or* between 0 and 1 (see above) corresponding to +#' the probability of direction (pd). #' #' @seealso [pd_to_p()] to convert between Probability of Direction (pd) and p-value. #' #' @note There is also a [`plot()`-method](https://easystats.github.io/see/articles/bayestestR.html) implemented in the \href{https://easystats.github.io/see/}{\pkg{see}-package}. #' #' @references -#' Makowski D, Ben-Shachar MS, Chen SHA, Lüdecke D (2019) Indices of Effect -#' Existence and Significance in the Bayesian Framework. Frontiers in Psychology -#' 2019;10:2767. \doi{10.3389/fpsyg.2019.02767} +#' - Makowski, D., Ben-Shachar, M. S., Chen, S. A., & Lüdecke, D. (2019). +#' Indices of effect existence and significance in the Bayesian framework. +#' Frontiers in psychology, 10, 2767. \doi{10.3389/fpsyg.2019.02767} +#' - van den Bergh, D., Haaf, J. M., Ly, A., Rouder, J. N., & Wagenmakers, E. J. +#' (2021). A cautionary note on estimating effect size. Advances in Methods +#' and Practices in Psychological Science, 4(1). \doi{10.1177/2515245921992035} #' #' @examples #' library(bayestestR) @@ -144,52 +166,10 @@ p_direction.default <- function(x, ...) { #' @rdname p_direction #' @export p_direction.numeric <- function(x, method = "direct", null = 0, ...) { - if (method == "direct") { - pdir <- max( - c( - length(x[x > null]) / length(x), # pd positive - length(x[x < null]) / length(x) # pd negative - ) - ) - } else { - dens <- estimate_density(x, method = method, precision = 2^10, extend = TRUE, ...) - if (length(x[x > null]) > length(x[x < null])) { - dens <- dens[dens$x > null, ] - } else { - dens <- dens[dens$x < null, ] - } - pdir <- area_under_curve(dens$x, dens$y, method = "spline") - if (pdir >= 1) pdir <- 1 # Enforce bounds - } - - attr(pdir, "method") <- method - attr(pdir, "data") <- x - - class(pdir) <- unique(c("p_direction", "see_p_direction", class(pdir))) - - pdir -} - - - - - -#' @export -p_direction.parameters_model <- function(x, ...) { - out <- data.frame( - "Parameter" = x$Parameter, - "pd" = p_to_pd(p = x[["p"]]), - row.names = NULL, - stringsAsFactors = FALSE - ) - - if (!is.null(x$Component)) { - out$Component <- x$Component - } - - attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x)) - class(out) <- unique(c("p_direction", "see_p_direction", class(out))) - + obj_name <- insight::safe_deparse_symbol(substitute(x)) + out <- p_direction(data.frame(x = x), method = method, null = null, ...) + out[[1]] <- NULL + attr(out, "object_name") <- obj_name out } @@ -201,9 +181,9 @@ p_direction.data.frame <- function(x, method = "direct", null = 0, ...) { x <- .select_nums(x) if (ncol(x) == 1) { - pd <- p_direction(x[[1]], method = method, null = null, ...) + pd <- .p_direction(x[[1]], method = method, null = null, ...) } else { - pd <- sapply(x, p_direction, method = method, null = null, simplify = TRUE, ...) + pd <- sapply(x, .p_direction, method = method, null = null, simplify = TRUE, ...) } out <- data.frame( @@ -465,6 +445,49 @@ p_direction.get_predicted <- function(x, ...) { out } +#' @export +p_direction.parameters_model <- function(x, ...) { + out <- data.frame( + "Parameter" = x$Parameter, + "pd" = p_to_pd(p = x[["p"]]), + row.names = NULL, + stringsAsFactors = FALSE + ) + + if (!is.null(x$Component)) { + out$Component <- x$Component + } + + attr(out, "object_name") <- insight::safe_deparse_symbol(substitute(x)) + class(out) <- unique(c("p_direction", "see_p_direction", class(out))) + + out +} + +#' @keywords internal +.p_direction <- function(x, method = "direct", null = 0, ...) { + if (method == "direct") { + pdir <- max( + length(x[x > null]), # pd positive + length(x[x < null]) # pd negative + ) / length(x) + } else { + dens <- estimate_density(x, method = method, precision = 2^10, extend = TRUE, ...) + if (length(x[x > null]) > length(x[x < null])) { + dens <- dens[dens$x > null, ] + } else { + dens <- dens[dens$x < null, ] + } + pdir <- area_under_curve(dens$x, dens$y, method = "spline") + if (pdir >= 1) { + # Enforce bounds + pdir <- 1 + } + } + + pdir +} + # Methods ----------------------------------------------------------------- diff --git a/R/p_map.R b/R/p_map.R index 0e6bc7485..d16f1dc11 100644 --- a/R/p_map.R +++ b/R/p_map.R @@ -82,17 +82,9 @@ p_pointnull <- p_map #' @export p_map.numeric <- function(x, null = 0, precision = 2^10, method = "kernel", ...) { - # Density at MAP - map <- attributes(map_estimate(x, precision = precision, method = method, ...))$MAP_density - - # Density at 0 - d_0 <- density_at(x, null, precision = precision, method = method, ...) - if (is.na(d_0)) d_0 <- 0 - - # Odds - p <- d_0 / map - class(p) <- c("p_map", class(p)) - p + out <- p_map(data.frame(x = x), null = null, precision = precision, method = method, ...) + out[[1]] <- NULL + out } @@ -102,9 +94,9 @@ p_map.data.frame <- function(x, null = 0, precision = 2^10, method = "kernel", . x <- .select_nums(x) if (ncol(x) == 1) { - p_MAP <- p_map(x[, 1], null = null, precision = precision, method = method, ...) + p_MAP <- .p_map(x[, 1], null = null, precision = precision, method = method, ...) } else { - p_MAP <- sapply(x, p_map, null = null, precision = precision, method = method, simplify = TRUE, ...) + p_MAP <- sapply(x, .p_map, null = null, precision = precision, method = method, simplify = TRUE, ...) } out <- data.frame( @@ -364,7 +356,19 @@ p_map.bayesQR <- function(x, null = 0, precision = 2^10, method = "kernel", ...) out } +#' @keywords internal +.p_map <- function(x, null = 0, precision = 2^10, method = "kernel", ...) { + # Density at MAP + map <- attributes(map_estimate(x, precision = precision, method = method, ...))$MAP_density + # Density at 0 + d_0 <- density_at(x, null, precision = precision, method = method, ...) + if (is.na(d_0)) d_0 <- 0 + + # Odds + p <- d_0 / map + p +} #' @rdname as.numeric.p_direction diff --git a/R/p_significance.R b/R/p_significance.R index d4ca4fa71..5432c3ab6 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -60,20 +60,10 @@ p_significance.default <- function(x, ...) { #' @export p_significance.numeric <- function(x, threshold = "default", ...) { threshold <- .select_threshold_ps(threshold = threshold) - - psig <- max( - c( - length(x[x > abs(threshold)]) / length(x), # ps positive - length(x[x < -abs(threshold)]) / length(x) # ps negative - ) - ) - - attr(psig, "threshold") <- threshold - attr(psig, "data") <- x - - class(psig) <- unique(c("p_significance", "see_p_significance", class(psig))) - - psig + out <- p_significance(data.frame(x = x), threshold = threshold) + out[[1]] <- NULL + attr(out, "data") <- x + out } @@ -84,9 +74,9 @@ p_significance.data.frame <- function(x, threshold = "default", ...) { x <- .select_nums(x) if (ncol(x) == 1) { - ps <- p_significance(x[, 1], threshold = threshold, ...) + ps <- .p_significance(x[, 1], threshold = threshold, ...) } else { - ps <- sapply(x, p_significance, threshold = threshold, simplify = TRUE, ...) + ps <- sapply(x, .p_significance, threshold = threshold, simplify = TRUE, ...) } out <- data.frame( @@ -260,7 +250,16 @@ p_significance.brmsfit <- function(x, out } +.p_significance <- function(x, threshold, ...) { + psig <- max( + c( + length(x[x > abs(threshold)]) / length(x), # ps positive + length(x[x < -abs(threshold)]) / length(x) # ps negative + ) + ) + psig +} # methods --------------------------- diff --git a/man/p_direction.Rd b/man/p_direction.Rd index cf90a3531..fec5356bc 100644 --- a/man/p_direction.Rd +++ b/man/p_direction.Rd @@ -48,13 +48,16 @@ pd(x, ...) \method{p_direction}{BFBayesFactor}(x, method = "direct", null = 0, ...) } \arguments{ -\item{x}{Vector representing a posterior distribution. Can also be a Bayesian model (\code{stanreg}, \code{brmsfit} or \code{BayesFactor}).} +\item{x}{A vector representing a posterior distribution, a data frame of +posterior draws (samples be parameter). Can also be a Bayesian model.} \item{...}{Currently not used.} -\item{method}{Can be \code{"direct"} or one of methods of \link[=estimate_density]{density estimation}, such as \code{"kernel"}, \code{"logspline"} or \code{"KernSmooth"}. If \code{"direct"} (default), the computation is based on the raw ratio of samples superior and inferior to 0. Else, the result is based on the \link[=auc]{Area under the Curve (AUC)} of the estimated \link[=estimate_density]{density} function.} +\item{method}{Can be \code{"direct"} or one of methods of \code{\link[=estimate_density]{estimate_density()}}, +such as \code{"kernel"}, \code{"logspline"} or \code{"KernSmooth"}. See details.} -\item{null}{The value considered as a "null" effect. Traditionally 0, but could also be 1 in the case of ratios.} +\item{null}{The value considered as a "null" effect. Traditionally 0, but +could also be 1 in the case of ratios of change (OR, IRR, ...).} \item{effects}{Should results for fixed effects, random effects or both be returned? Only applies to mixed models. May be abbreviated.} @@ -70,65 +73,95 @@ filtered by default, so only parameters that typically appear in the for the output.} } \value{ -Values between 0.5 and 1 corresponding to the probability of direction (pd). -\cr\cr -Note that in some (rare) cases, especially when used with model averaged -posteriors (see \code{\link[=weighted_posteriors]{weighted_posteriors()}} or -\code{brms::posterior_average}), \code{pd} can be smaller than \code{0.5}, -reflecting high credibility of \code{0}. To detect such cases, the -\code{method = "direct"} must be used. +Values between 0.5 and 1 \emph{or} between 0 and 1 (see above) corresponding to +the probability of direction (pd). } \description{ -Compute the \strong{Probability of Direction} (\emph{\strong{pd}}, also known -as the Maximum Probability of Effect - \emph{MPE}). It varies between \verb{50\%} -and \verb{100\%} (\emph{i.e.}, \code{0.5} and \code{1}) and can be interpreted as -the probability (expressed in percentage) that a parameter (described by its -posterior distribution) is strictly positive or negative (whichever is the -most probable). It is mathematically defined as the proportion of the -posterior distribution that is of the median's sign. Although differently -expressed, this index is fairly similar (\emph{i.e.}, is strongly correlated) -to the frequentist \strong{p-value}. -\cr\cr -Note that in some (rare) cases, especially when used with model averaged -posteriors (see \code{\link[=weighted_posteriors]{weighted_posteriors()}} or -\code{brms::posterior_average}), \code{pd} can be smaller than \code{0.5}, -reflecting high credibility of \code{0}. +Compute the \strong{Probability of Direction} (\emph{\strong{pd}}, also known as the Maximum +Probability of Effect - \emph{MPE}). This can be interpreted as the probability +that a parameter (described by its posterior distribution) is strictly +positive or negative (whichever is the most probable). Although differently +expressed, this index is fairly similar (\emph{i.e.}, is strongly correlated) to +the frequentist \strong{p-value} (see details). } \details{ \subsection{What is the \emph{pd}?}{ -The Probability of Direction (pd) is an index of effect existence, ranging -from \verb{50\%} to \verb{100\%}, representing the certainty with which an effect goes in -a particular direction (\emph{i.e.}, is positive or negative). Beyond its -simplicity of interpretation, understanding and computation, this index also -presents other interesting properties: + +The Probability of Direction (pd) is an index of effect existence, representing +the certainty with which an effect goes in a particular direction (i.e., is +positive or negative / has a sign), typically ranging from 0.5 to 1 (but see +next section for cases where it can range between 0 and 1). Beyond +its simplicity of interpretation, understanding and computation, this index +also presents other interesting properties: \itemize{ -\item It is independent from the model: It is solely based on the posterior +\item Like other posterior-based indices, \emph{pd} is solely based on the posterior distributions and does not require any additional information from the data -or the model. +or the model (e.g., such as priors, as in the case of Bayes factors). \item It is robust to the scale of both the response variable and the predictors. \item It is strongly correlated with the frequentist p-value, and can thus be used to draw parallels and give some reference to readers non-familiar -with Bayesian statistics. +with Bayesian statistics (Makowski et al., 2019). } } + \subsection{Relationship with the p-value}{ -In most cases, it seems that the \emph{pd} has a direct correspondence with the frequentist one-sided \emph{p}-value through the formula \ifelse{html}{\out{pone sided = 1 - p(d)/100}}{\eqn{p_{one sided}=1-\frac{p_{d}}{100}}} and to the two-sided p-value (the most commonly reported one) through the formula \ifelse{html}{\out{ptwo sided = 2 * (1 - p(d)/100)}}{\eqn{p_{two sided}=2*(1-\frac{p_{d}}{100})}}. Thus, a two-sided p-value of respectively \code{.1}, \code{.05}, \code{.01} and \code{.001} would correspond approximately to a \emph{pd} of \verb{95\%}, \verb{97.5\%}, \verb{99.5\%} and \verb{99.95\%}. See also \code{\link[=pd_to_p]{pd_to_p()}}. + +In most cases, it seems that the \emph{pd} has a direct correspondence with the +frequentist one-sided \emph{p}-value through the formula (for two-sided \emph{p}): +\deqn{p = 2 \times (1 - p_d)}{p = 2 * (1 - pd)} +Thus, a two-sided p-value of respectively \code{.1}, \code{.05}, \code{.01} and \code{.001} would +correspond approximately to a \emph{pd} of \verb{95\%}, \verb{97.5\%}, \verb{99.5\%} and \verb{99.95\%}. +See \code{\link[=pd_to_p]{pd_to_p()}} for details. } + +\subsection{Possible Range of Values}{ + +The largest value \emph{pd} can take is 1 - the posterior is strictly directional. +However, the smallest value \emph{pd} can take depends on the parameter space +represented by the posterior. +\cr\cr +\strong{For a continuous parameter space}, exact values of 0 (or any point null +value) are not possible, and so 100\% of the posterior has \emph{some} sign, some +positive, some negative. Therefore, the smallest the \emph{pd} can be is 0.5 - +with an equal posterior mass of positive and negative values. Values close to +0.5 \emph{cannot} be used to support the null hypothesis (that the parameter does +\emph{not} have a direction) is a similar why to how large p-values cannot be used +to support the null hypothesis (see \code{\link[=pd_tp_p]{pd_tp_p()}}; Makowski et al., 2019). +\cr\cr +\strong{For a discrete parameter space or a parameter space that is a mixture +between discrete and continuous spaces}, exact values of 0 (or any point +null value) \emph{are} possible! Therefore, the smallest the \emph{pd} can be is 0 - +with 100\% of the posterior mass on 0. Thus values close to 0 can be used to +support the null hypothesis (see van den Bergh et al., 2021). +\cr\cr +Examples of posteriors representing discrete parameter space: +\itemize{ +\item When a parameter can only take discrete values. +\item When a mixture prior/posterior is used (such as the spike-and-slab prior; +see van den Bergh et al., 2021). +\item When conducting Bayesian model averaging (e.g., \code{\link[=weighted_posteriors]{weighted_posteriors()}} or +\code{brms::posterior_average}). +} +} + \subsection{Methods of computation}{ -The most simple and direct way to compute the \emph{pd} is to 1) look at the -median's sign, 2) select the portion of the posterior of the same sign and -3) compute the percentage that this portion represents. This "simple" method -is the most straightforward, but its precision is directly tied to the -number of posterior draws. The second approach relies on \link[=estimate_density]{density estimation}. It starts by estimating the density function -(for which many methods are available), and then computing the \link[=area_under_curve]{area under the curve} (AUC) of the density curve on the other side of -0. -} -\subsection{Strengths and Limitations}{ -\strong{Strengths:} Straightforward computation and interpretation. Objective -property of the posterior distribution. 1:1 correspondence with the -frequentist p-value. -\cr \cr -\strong{Limitations:} Limited information favoring the null hypothesis. + +The \emph{pd} is defined as: +\deqn{p_d = max({Pr(\hat{\theta} < \theta_{null}), Pr(\hat{\theta} > \theta_{null})})}{pd = max(mean(x < null), mean(x > null))} +\cr\cr +The most simple and direct way to compute the \emph{pd} is to compute the +proportion of positive (or larger than \code{null}) posterior samples, the +proportion of negative (or smaller than \code{null}) posterior samples, and take +the larger of the two. This "simple" method is the most straightforward, but +its precision is directly tied to the number of posterior draws. +\cr\cr +The second approach relies on \code{\link[=density estimation]{density estimation()}}: It starts by +estimating the continuous-smooth density function (for which many methods are +available), and then computing the \link[=area_under_curve]{area under the curve} +(AUC) of the density curve on either side of \code{null} and taking the maximum +between them. Note the this approach assumes a continuous density function, +and so \strong{when the posterior represents a (partially) discrete parameter +space, only the direct method \emph{must} be used} (see above). } } \note{ @@ -184,9 +217,14 @@ if (require("BayesFactor")) { } } \references{ -Makowski D, Ben-Shachar MS, Chen SHA, Lüdecke D (2019) Indices of Effect -Existence and Significance in the Bayesian Framework. Frontiers in Psychology -2019;10:2767. \doi{10.3389/fpsyg.2019.02767} +\itemize{ +\item Makowski, D., Ben-Shachar, M. S., Chen, S. A., & Lüdecke, D. (2019). +Indices of effect existence and significance in the Bayesian framework. +Frontiers in psychology, 10, 2767. \doi{10.3389/fpsyg.2019.02767} +\item van den Bergh, D., Haaf, J. M., Ly, A., Rouder, J. N., & Wagenmakers, E. J. +(2021). A cautionary note on estimating effect size. Advances in Methods +and Practices in Psychological Science, 4(1). \doi{10.1177/2515245921992035} +} } \seealso{ \code{\link[=pd_to_p]{pd_to_p()}} to convert between Probability of Direction (pd) and p-value. diff --git a/man/p_map.Rd b/man/p_map.Rd index 9e9e625fe..94d9c21f2 100644 --- a/man/p_map.Rd +++ b/man/p_map.Rd @@ -41,7 +41,8 @@ of models (see, for example, \code{methods("hdi")}) and not all of those are documented in the 'Usage' section, because methods for other classes mostly resemble the arguments of the \code{.numeric} or \code{.data.frame}methods.} -\item{null}{The value considered as a "null" effect. Traditionally 0, but could also be 1 in the case of ratios.} +\item{null}{The value considered as a "null" effect. Traditionally 0, but +could also be 1 in the case of ratios of change (OR, IRR, ...).} \item{precision}{Number of points of density data. See the \code{n} parameter in \code{density}.} diff --git a/man/pd_to_p.Rd b/man/pd_to_p.Rd index 19c229157..d33e1f961 100644 --- a/man/pd_to_p.Rd +++ b/man/pd_to_p.Rd @@ -7,13 +7,13 @@ \alias{convert_pd_to_p} \title{Convert between Probability of Direction (pd) and p-value.} \usage{ -pd_to_p(pd, direction = "two-sided", ...) +pd_to_p(pd, direction = "two-sided", verbose = TRUE, ...) p_to_pd(p, direction = "two-sided", ...) convert_p_to_pd(p, direction = "two-sided", ...) -convert_pd_to_p(pd, direction = "two-sided", ...) +convert_pd_to_p(pd, direction = "two-sided", verbose = TRUE, ...) } \arguments{ \item{pd}{A Probability of Direction (pd) value (between 0 and 1).} @@ -21,6 +21,8 @@ convert_pd_to_p(pd, direction = "two-sided", ...) \item{direction}{What type of p-value is requested or provided. Can be \code{"two-sided"} (default, two tailed) or \code{"one-sided"} (one tailed).} +\item{verbose}{Toggle off warnings.} + \item{...}{Arguments passed to or from other methods.} \item{p}{A p-value.} @@ -28,7 +30,28 @@ convert_pd_to_p(pd, direction = "two-sided", ...) \description{ Enables a conversion between Probability of Direction (pd) and p-value. } +\details{ +Conversion is done using the following equation (see Makowski et al., 2019): +\cr\cr +When \code{direction = "two-sided"} - +\cr\cr +\deqn{p = 2 \times (1 - p_d)}{p = 2 * (1 - pd)} +When \code{direction = "one-sided"} - +\cr\cr +\deqn{p = 1 - p_d}{p = 1 - pd} +\cr\cr +Note that this conversion is only valid when the lowest possible values of pd +is 0.5 - i.e., when the posterior represents continuous parameter space (see +\link{p_direction}). If any pd < 0.5 are detected, they are converted to a p of 1, +and a warning is given. +} \examples{ pd_to_p(pd = 0.95) pd_to_p(pd = 0.95, direction = "one-sided") + +} +\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} } diff --git a/man/sexit.Rd b/man/sexit.Rd index 8b0f8029a..74d4cd5ce 100644 --- a/man/sexit.Rd +++ b/man/sexit.Rd @@ -7,7 +7,8 @@ sexit(x, significant = "default", large = "default", ci = 0.95, ...) } \arguments{ -\item{x}{Vector representing a posterior distribution. Can also be a Bayesian model (\code{stanreg}, \code{brmsfit} or \code{BayesFactor}).} +\item{x}{A vector representing a posterior distribution, a data frame of +posterior draws (samples be parameter). Can also be a Bayesian model.} \item{significant, large}{The threshold values to use for significant and large probabilities. If left to 'default', will be selected through diff --git a/tests/testthat/test-bayesfactor_models.R b/tests/testthat/test-bayesfactor_models.R index 24fe273ae..27294b874 100644 --- a/tests/testthat/test-bayesfactor_models.R +++ b/tests/testthat/test-bayesfactor_models.R @@ -118,8 +118,6 @@ test_that("bayesfactor_models BRMS", { skip_on_cran() skip_on_ci() - skip_if_not_or_load_if_installed("lme4") - skip_if_not_or_load_if_installed("rstanarm") skip_if_not_or_load_if_installed("bridgesampling") skip_if_not_or_load_if_installed("brms") diff --git a/tests/testthat/test-format.R b/tests/testthat/test-format.R index c8e0377a5..08297ed2a 100644 --- a/tests/testthat/test-format.R +++ b/tests/testthat/test-format.R @@ -23,12 +23,12 @@ test_that("p_significance", { ) expect_equal( format(p_direction(x)), - data.frame(x = "0.51", stringsAsFactors = FALSE), + data.frame(x = "51.00%", stringsAsFactors = FALSE), ignore_attr = TRUE ) expect_equal( format(p_map(x)), - data.frame(x = "0.97", stringsAsFactors = FALSE), + data.frame(x = "0.973", stringsAsFactors = FALSE), ignore_attr = TRUE ) expect_equal( diff --git a/tests/testthat/test-map_estimate.R b/tests/testthat/test-map_estimate.R index f66431818..f99b6c764 100644 --- a/tests/testthat/test-map_estimate.R +++ b/tests/testthat/test-map_estimate.R @@ -1,13 +1,12 @@ # numeric ---------------------- test_that("map_estimate", { - skip_if_offline() - skip_if_not_or_load_if_installed("rstanarm") - - expect_equal( - as.numeric(map_estimate(distribution_normal(1000))), - 0, - tolerance = 0.01 - ) + x <- distribution_normal(1000, 1) + MAP <- map_estimate(x) + expect_equal(as.numeric(MAP), 0.997, tolerance = 0.001) + expect_s3_class(MAP, "map_estimate") + expect_s3_class(MAP, "data.frame") + expect_equal(dim(MAP), c(1, 1)) + expect_output(print(MAP), "MAP Estimate: 1.00") }) # stanreg ---------------------- diff --git a/tests/testthat/test-p_direction.R b/tests/testthat/test-p_direction.R index f0c182d71..6d38d0d7c 100644 --- a/tests/testthat/test-p_direction.R +++ b/tests/testthat/test-p_direction.R @@ -1,15 +1,19 @@ test_that("p_direction", { - skip_if_not_or_load_if_installed("rstanarm") - skip_if_not_or_load_if_installed("brms") - set.seed(333) x <- distribution_normal(10000, 1, 1) pd <- p_direction(x) expect_equal(as.numeric(pd), 0.842, tolerance = 0.1) expect_equal(as.numeric(p_direction(x, method = "kernel")), 0.842, tolerance = 0.1) - expect_equal(nrow(p_direction(data.frame(replicate(4, rnorm(100))))), 4) expect_s3_class(pd, "p_direction") - expect_equal(tail(capture.output(print(pd)), 1), "Probability of Direction: 0.84") + expect_s3_class(pd, "data.frame") + expect_equal(dim(pd), c(1L, 1L)) + expect_output(print(pd), regexp = "Probability of Direction: 84.13%", fixed = TRUE) + + df <- data.frame(replicate(4, rnorm(100))) + pd <- p_direction(df) + expect_s3_class(pd, "p_direction") + expect_s3_class(pd, "data.frame") + expect_equal(dim(pd), c(4, 2)) }) diff --git a/tests/testthat/test-p_map.R b/tests/testthat/test-p_map.R index f47f4703f..88514c446 100644 --- a/tests/testthat/test-p_map.R +++ b/tests/testthat/test-p_map.R @@ -1,6 +1,11 @@ test_that("p_map", { - skip_if_offline() - skip_if_not_or_load_if_installed("rstanarm") + x <- distribution_normal(1000, 0.4) + pmap <- p_map(x) + expect_equal(as.numeric(pmap), 0.9285376, tolerance = 0.001) + expect_s3_class(pmap, "p_map") + expect_s3_class(pmap, "data.frame") + expect_equal(dim(pmap), c(1, 1)) + expect_output(print(pmap), "MAP-based p-value: 0.929") expect_equal(as.numeric(p_map(distribution_normal(1000))), 1, tolerance = 0.1) expect_equal(as.numeric(p_map(distribution_normal(1000, 1, 1))), 0.62, tolerance = 0.1) @@ -42,6 +47,6 @@ test_that("p_map", { test_that("p_map | null", { x <- distribution_normal(4000, mean = 1) - expect_equal(p_map(x), 0.6194317, ignore_attr = TRUE, tolerance = 0.01) - expect_equal(p_map(x, null = 1), 1, ignore_attr = TRUE, tolerance = 0.01) + expect_equal(as.numeric(p_map(x)), 0.6194317, ignore_attr = TRUE, tolerance = 0.01) + expect_equal(as.numeric(p_map(x, null = 1)), 1, ignore_attr = TRUE, tolerance = 0.01) }) diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R index 1423a67e9..8463a06db 100644 --- a/tests/testthat/test-p_significance.R +++ b/tests/testthat/test-p_significance.R @@ -1,15 +1,17 @@ test_that("p_significance", { - skip_if_offline() - skip_if_not_or_load_if_installed("rstanarm") - # numeric set.seed(333) x <- distribution_normal(10000, 1, 1) ps <- p_significance(x) expect_equal(as.numeric(ps), 0.816, tolerance = 0.1) - expect_equal(nrow(p_significance(data.frame(replicate(4, rnorm(100))))), 4) expect_s3_class(ps, "p_significance") - expect_equal(tail(capture.output(print(ps)), 1), "Practical Significance (threshold: 0.10): 0.82") + expect_s3_class(ps, "data.frame") + expect_equal(dim(ps), c(1, 1)) + expect_output(print(ps), "Practical Significance \\(threshold: 0.10\\): 0.82") + + x <- data.frame(replicate(4, rnorm(100))) + pd <- p_significance(x) + expect_equal(dim(pd), c(4, 2)) }) test_that("stanreg", { diff --git a/tests/testthat/test-pd_to_p.R b/tests/testthat/test-pd_to_p.R new file mode 100644 index 000000000..3fc513b10 --- /dev/null +++ b/tests/testthat/test-pd_to_p.R @@ -0,0 +1,8 @@ +test_that("pd_to_p", { + pds <- c(0.7, 0.95, 0.99, 0.5) + expect_equal(pd_to_p(pds), c(0.6, 0.1, 0.02, 1)) + expect_equal(pd_to_p(pds, direction = 1), c(0.3, 0.05, 0.01, 0.5)) + + expect_warning(p <- pd_to_p(0.3), "0.5") + expect_equal(p, 1) +})