R/wilcox_effsize.R

In rstatix: Pipe-Friendly Framework for Basic Statistical Tests

#' @include utilities.R utilities_two_sample_test.R
NULL

#'Wilcoxon Effect Size
#'@description Compute Wilcoxon effect size (\code{r}) for: \itemize{ \item
#'  one-sample test (Wilcoxon one-sample signed-rank test); \item paired
#'  two-samples test (Wilcoxon two-sample paired signed-rank test) and \item
#'  independent two-samples test ( Mann-Whitney, two-sample rank-sum test). }
#'
#'  It can also returns confidence intervals by bootstap.
#'
#'  The effect size \code{r} is calculated as \code{Z} statistic divided by
#'  square root of the sample size (N) (\eqn{Z/\sqrt{N}}). The \code{Z} value is
#'  extracted from either \code{coin::wilcoxsign_test()} (case of one- or
#'  paired-samples test) or \code{coin::wilcox_test()} (case of independent
#'  two-samples test).
#'
#'  Here, \code{N} is the number of independent observations contributing to the
#'  test: the total sample size for the independent two-samples test, and the
#'  \strong{number of pairs} (equivalently, the number of difference scores) for
#'  the one-sample and paired tests. This is because the paired test reduces to a
#'  one-sample signed-rank test on the pairwise differences, so each pair counts
#'  once. This convention matches the default of
#'  \code{rcompanion::wilcoxonPairedR()} (its \code{cases = TRUE} setting).
#'
#'  Some references instead define \code{N} as the total number of observations,
#'  i.e. twice the number of pairs (Field, 2012; Tomczak & Tomczak, 2014), which
#'  yields a smaller \code{r}. If you need that convention for a paired test,
#'  divide the reported \code{r} (or the \code{Z}) by \eqn{\sqrt 2}; it is also
#'  available via \code{rcompanion::wilcoxonPairedR(..., cases = FALSE)}.
#'
#'  The \code{r} value varies from 0 to close to 1. The interpretation values
#'  for r commonly in published litterature and on the internet are: \code{0.10
#'  - < 0.3} (small effect), \code{0.30 - < 0.5} (moderate effect) and \code{>=
#'  0.5} (large effect).
#'
#'  See the Datanovia tutorial
#'  \href{https://www.datanovia.com/learn/biostatistics/two-groups/wilcoxon-test-in-r}{Wilcoxon Test in R}
#'  for a worked walkthrough.
#'
#'@inheritParams wilcox_test
#'@param ci If TRUE, returns confidence intervals by bootstrap. May be slow.
#'@param conf.level The level for the confidence interval.
#'@param ci.type The type of confidence interval to use. Can be any of "norm",
#'  "basic", "perc", or "bca". Passed to \code{boot::boot.ci}.
#'@param nboot The number of replications to use for bootstrap.
#'@param boot.parallel The type of parallel operation to be used when computing
#'  the bootstrap confidence interval. Allowed values are \code{"no"} (default),
#'  \code{"multicore"} and \code{"snow"}. Passed to \code{\link[boot]{boot}()}.
#'  Defaults to \code{getOption("boot.parallel", "no")}, so it can also be set
#'  globally with \code{options(boot.parallel = "multicore")}. Only used when
#'  \code{ci = TRUE}.
#'@param boot.ncpus Integer. The number of processes to be used in the parallel
#'  bootstrap. Defaults to \code{getOption("boot.ncpus", 1L)}. Note that
#'  \code{boot.parallel} has no effect unless \code{boot.ncpus > 1}. Only used
#'  when \code{ci = TRUE}.
#'@param method the effect-size metric. Either \code{"r"} (default) for the
#'  rank correlation \code{r = Z / sqrt(N)}, or \code{"rank_biserial"} for the
#'  rank-biserial correlation --- Cliff's delta for an independent-samples test
#'  (equal to \code{\link{cliff_delta}()}) or the matched-pairs rank-biserial for
#'  a paired test, both equal to \code{effectsize::rank_biserial()}. The
#'  independent-samples case is labelled with the Romano et al. magnitude
#'  thresholds (those thresholds define Cliff's delta); the paired case carries
#'  no \code{magnitude} (\code{NA}), because no threshold set is calibrated for
#'  the matched-pairs rank-biserial. The confidence interval (\code{ci = TRUE})
#'  is a percentile bootstrap.
#'@param detailed logical value. Default is FALSE. If TRUE, and
#'  \code{method = "r"}, the output additionally includes the \code{Z}
#'  \code{statistic} (extracted from the \code{coin} package and used to compute
#'  \code{r = Z/sqrt(N)}), the p-value (\code{p}) and the test
#'  \code{method}/\code{alternative}, so the effect size and the underlying Z are
#'  reported together in one data frame. The rank-biserial metric
#'  (\code{method = "rank_biserial"}) has no underlying \code{Z}, so those extra
#'  columns are not meaningful for it.
#'@param ... Additional arguments passed to the functions
#'  \code{coin::wilcoxsign_test()} (case of one- or paired-samples test) or
#'  \code{coin::wilcox_test()} (case of independent two-samples test).
#'@return return a data frame with some of the following columns: \itemize{
#'  \item \code{.y.}: the y variable used in the test. \item
#'  \code{group1,group2}: the compared groups in the pairwise tests. \item
#'  \code{n,n1,n2}: Sample counts. \item \code{effsize}: estimate of the effect
#'  size (\code{r} value). \item \code{magnitude}: magnitude of effect size.
#'  \item \code{conf.low,conf.high}: lower and upper bound of the effect size
#'  confidence interval. \item \code{statistic}: the \code{Z} statistic and
#'  \code{p}: the p-value (only when \code{detailed = TRUE}).}
#'@references Maciej Tomczak and Ewa Tomczak. The need to report effect size
#'  estimates revisited. An overview of some recommended measures of effect
#'  size. Trends in Sport Sciences. 2014; 1(21):19-25.
#' @examples
#' if(require("coin")){
#'
#' # One-sample Wilcoxon test effect size
#' ToothGrowth %>% wilcox_effsize(len ~ 1, mu = 0)
#'
#' # Independent two-samples wilcoxon effect size
#' ToothGrowth %>% wilcox_effsize(len ~ supp)
#'
#'
#' # Paired-samples wilcoxon effect size
#' ToothGrowth %>% wilcox_effsize(len ~ supp, paired = TRUE)
#'
#' # Pairwise comparisons
#' ToothGrowth %>% wilcox_effsize(len ~ dose)
#'
#' # Grouped data
#' ToothGrowth %>%
#'   group_by(supp) %>%
#'   wilcox_effsize(len ~ dose)
#'
#' }
#' @seealso The Datanovia tutorial: \href{https://www.datanovia.com/learn/biostatistics/two-groups/wilcoxon-test-in-r}{Wilcoxon Test in R}.
#'@export
wilcox_effsize <- function(data, formula, comparisons = NULL, ref.group = NULL,
                                paired = FALSE, alternative = "two.sided",
                                mu = 0, ci = FALSE, conf.level = 0.95, ci.type = "perc",
                                nboot = 1000, detailed = FALSE, ...,
                                boot.parallel = getOption("boot.parallel", "no"),
                                boot.ncpus = getOption("boot.ncpus", 1L),
                                method = c("r", "rank_biserial")){
  method <- match.arg(method)
  env <- as.list(environment()) %>% remove_item("method")
  # See cohens_d(): the bootstrap-execution arguments are not part of the
  # statistical call, so they are excluded from the stashed args. `method` is
  # kept out of the stashed args unless it is the non-default value, so the
  # default (r = Z/sqrt(N)) leaves attr(x, "args") unchanged.
  args <- env %>%
    remove_item(c("boot.parallel", "boot.ncpus")) %>%
    .add_item(method = "wilcox_effsize")
  if(method != "r") args <- args %>% .add_item(effsize.method = method)
  # method = "r": Z/sqrt(N) via coin (unchanged). method = "rank_biserial": the
  # rank-biserial correlation -- Cliff's delta for an independent test (equals
  # cliff_delta()) or the matched-pairs rank-biserial for a paired test -- with
  # the Romano magnitude thresholds those metrics use.
  stat.method <- if(method == "r") "coin.wilcox.test"
                 else if(isTRUE(paired)) "rank.biserial" else "cliff.delta"
  # An independent-samples rank_biserial IS Cliff's delta, so its Romano
  # thresholds apply. The matched-pairs rank-biserial (paired) has no calibrated
  # threshold set, so -- like wilcox_test(effect.size = TRUE, paired = TRUE) --
  # no magnitude is assigned (the column stays, filled with NA).
  magnitude.fun <- if(method == "r") get_wilcox_effsize_magnitude
                   else if(isTRUE(paired)) no_effsize_magnitude
                   else get_cliff_delta_magnitude
  params <- c(env, list(...)) %>%
    remove_null_items() %>%
    add_item(method = stat.method, detailed = detailed)

  outcome <- get_formula_left_hand_side(formula)
  group <- get_formula_right_hand_side(formula)
  number.of.groups <- guess_number_of_groups(data, group)
  if(number.of.groups > 2 & !is.null(ref.group)){
    if(ref.group %in% c("all", ".all.")){
      params$data <- create_data_with_all_ref_group(data, outcome, group)
      params$ref.group <- "all"
    }
  }
  test.func <- two_sample_test
  if(number.of.groups > 2) test.func <- pairwise_two_sample_test
  res <- do.call(test.func, params) %>%
    select(all_of(c(".y.", "group1", "group2", "estimate")), everything()) %>%
    rename(effsize = "estimate") %>%
    mutate(magnitude = magnitude.fun(.data$effsize)) %>%
    set_attrs(args = args) %>%
    add_class(c("rstatix_test", "wilcox_effsize"))
  if(identical(stat.method, "rank.biserial")) warn_undefined_rank_biserial(res)
  warn_undefined_boot_ci(res, ci)
  res
}



# Wilcoxon test using coin R package; returns effect size
coin.wilcox.test <- function(x, y = NULL, mu = 0, paired = FALSE, alternative = c("two.sided", "less", "greater"),
                     ci = FALSE, conf.level = 0.95,  ci.type = "perc", nboot = 1000, ...,
                     boot.parallel = getOption("boot.parallel", "no"),
                     boot.ncpus = getOption("boot.ncpus", 1L)){
  required_package("coin")

  alternative <- match.arg(alternative)
  check_two_samples_test_args(
    x = x, y = y, mu = mu, paired = paired,
    conf.level = conf.level
    )

  if (!is.null(y)) {
    DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
    if (paired) {
      # Transform paired test into one-sample test problem
      OK <- complete.cases(x, y)
      x <- x[OK] - y[OK]
      y <- NULL
      METHOD <- "Paired Wilcoxon test (coin)"
    }
    else {
      x <- x[is.finite(x)]
      y <- y[is.finite(y)]
      METHOD <- "Independent Wilcoxon test (coin)"
    }
  }
  else {
    DNAME <- deparse(substitute(x))
    METHOD <- "One-sample Wilcoxon test (coin)"
    x <- x[is.finite(x)]
  }

  if(is.null(y)){
    y <- rep(mu, length(x))
    test.type <- "symmetry"
  }
  else{
    group <- rep(c("grp1", "grp2"), times = c(length(x), length(y))) %>%
      factor()
    x <- c(x, y)
    y <- group
    test.type <- "independence"
  }
  data <- data.frame(x, y)
  results <- coin_wilcox_test(
    data, x ~ y, type = test.type,
    alternative = alternative,  ...
    )
  # Confidence interval of the effect size r
  if (ci == TRUE) {
    stat.func <- function(data, subset) {
      coin_wilcox_test(
        data, formula = x ~ y, subset = subset,
        type = test.type, alternative = alternative, ...
        )$r
    }
    CI <- get_boot_ci(
      data, stat.func, conf.level = conf.level,
      type = ci.type, nboot = nboot, parallel = boot.parallel, ncpus = boot.ncpus
      )
    results <- results %>% mutate(conf.low = CI[1], conf.high = CI[2])
  }
  # Note: no 'parameter' is set. The Wilcoxon test has no degrees of freedom; the
  # sample size is already reported as n1/n2, and tidying a 'parameter' here would
  # surface a spurious 'df' column equal to N in the detailed output (#122).
  RVAL <- list(statistic = results$z, p.value = results$p,
               null.value = mu, alternative = alternative, method = METHOD,
               data.name = DNAME, estimate = results$r)
  if (ci) {
    attr(CI, "conf.level") <- conf.level
    RVAL <- c(RVAL, list(conf.int = CI))
  }
  names(RVAL$statistic) <- "Z"
  names(RVAL$estimate) <- "Effect size (r)"
  class(RVAL) <- "htest"
  RVAL
}

# Perform wilcoxon test using coin package
coin_wilcox_test <- function(data, formula, subset = NULL, type = c("independence", "symmetry"),  ...){
  type <- match.arg(type)
  coin_wilcox_test_func <- switch (
    type,
    independence = coin::wilcox_test,
    symmetry = coin::wilcoxsign_test
  )
  if(!is.null(subset)) data <- data[subset, ]
  res.wilcox <-suppressWarnings(coin_wilcox_test_func(formula, data = data,...))
  n <- nrow(data)
  z <- as.vector(coin::statistic(res.wilcox, type = "standardized"))
  p <- coin::pvalue(res.wilcox)
  r <- abs(z)/sqrt(n) # Effect size
  tibble(n = n, z = z, r = r, p = p)
}

get_wilcox_effsize_magnitude <- function(d){
  magnitude.levels = c(0.3, 0.5, Inf)
  magnitude = c("small","moderate","large")
  d.index <- findInterval(abs(d), magnitude.levels)+1
  magnitude <- factor(magnitude[d.index], levels = magnitude, ordered = TRUE)
  magnitude
}

# No calibrated magnitude thresholds exist for the matched-pairs rank-biserial,
# so return an all-NA ordered factor rather than mis-applying the independent-
# sample Romano thresholds. The levels mirror get_cliff_delta_magnitude() (the
# independent rank_biserial case) so the magnitude column keeps one factor type
# whether the rank_biserial result is independent or paired.
no_effsize_magnitude <- function(d){
  factor(rep(NA_character_, length(d)),
         levels = c("negligible", "small", "medium", "large"), ordered = TRUE)
}