Specify reference level when comparing two groups/levels (#666)

Closes #665

Author: bwiernik

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: effectsize
 Title: Indices of Effect Size
-Version: 1.0.0
+Version: 1.0.0.1
 Authors@R: 
     c(person(given = "Mattan S.",
              family = "Ben-Shachar",

NEWS.md

@@ -1,3 +1,7 @@
+# effectsize 1.0.x
+
+- `cohens_d()`, `p_superiority()`, `rank_biserial()` and their relatives gain a `reference=` argument to control which level of the group variable should be treated as the reference (thanks @profandyfield for the suggestion).
+
 # effectsize 1.0.0
 
 ***First stable release of `{effectsize}`!***

R/cohens_d.R

@@ -1,6 +1,6 @@
 #' Cohen's *d* and Other Standardized Differences
 #'
-#' Compute effect size indices for standardized differences: Cohen's *d*,
+#' Compute effect size indices for standardized mean differences: Cohen's *d*,
 #' Hedges' *g* and Glass’s *delta* (\eqn{\Delta}). (This function returns the
 #' **population** estimate.) Pair with any reported [`stats::t.test()`].
 #' \cr\cr
@@ -9,7 +9,7 @@
 #' correction for small-sample bias (using the exact method) to Cohen's *d*. For
 #' sample sizes > 20, the results for both statistics are roughly equivalent.
 #' Glass’s *delta* is appropriate when the standard deviations are significantly
-#' different between the populations, as it uses only the *second* group's
+#' different between the populations, as it uses only the reference group's
 #' standard deviation.
 #'
 #' @param x,y A numeric vector, or a character name of one in `data`.
@@ -29,6 +29,8 @@
 #' @param adjust Should the effect size be adjusted for small-sample bias using
 #'   Hedges' method? Note that `hedges_g()` is an alias for
 #'   `cohens_d(adjust = TRUE)`.
+#' @param reference (Optional) character value of the "group" used as the
+#'   reference. By default, the _second_ group is the reference group.
 #' @param ... Arguments passed to or from other methods. When `x` is a formula,
 #'   these can be `subset` and `na.action`.
 #' @inheritParams chisq_to_phi
@@ -136,6 +138,7 @@
 #' @export
 cohens_d <- function(x, y = NULL, data = NULL,
                      pooled_sd = TRUE, mu = 0, paired = FALSE,
+                     reference = NULL,
                      adjust = FALSE,
                      ci = 0.95, alternative = "two.sided",
                      verbose = TRUE, ...) {
@@ -147,6 +150,7 @@ cohens_d <- function(x, y = NULL, data = NULL,
     y = y, data = data,
     type = "d", adjust = adjust,
     pooled_sd = pooled_sd, mu = mu, paired = paired,
+    reference = reference,
     ci = ci, alternative = alternative,
     verbose = verbose,
     ...
@@ -157,6 +161,7 @@ cohens_d <- function(x, y = NULL, data = NULL,
 #' @export
 hedges_g <- function(x, y = NULL, data = NULL,
                      pooled_sd = TRUE, mu = 0, paired = FALSE,
+                     reference = NULL,
                      ci = 0.95, alternative = "two.sided",
                      verbose = TRUE, ...) {
   cl <- match.call()
@@ -169,13 +174,15 @@ hedges_g <- function(x, y = NULL, data = NULL,
 #' @export
 glass_delta <- function(x, y = NULL, data = NULL,
                         mu = 0, adjust = TRUE,
+                        reference = NULL,
                         ci = 0.95, alternative = "two.sided",
                         verbose = TRUE, ...) {
   .effect_size_difference(
     x,
     y = y, data = data,
     type = "delta",
     mu = mu, adjust = adjust,
+    reference = reference,
     ci = ci, alternative = alternative,
     verbose = verbose,
     pooled_sd = NULL, paired = FALSE,
@@ -189,10 +196,12 @@ glass_delta <- function(x, y = NULL, data = NULL,
 .effect_size_difference <- function(x, y = NULL, data = NULL,
                                     type = "d", adjust = FALSE,
                                     mu = 0, pooled_sd = TRUE, paired = FALSE,
+                                    reference = NULL,
                                     ci = 0.95, alternative = "two.sided",
                                     verbose = TRUE, ...) {
   if (type == "d" && adjust) type <- "g"
 
+  # TODO: Check if we can do anything with `reference` for these classes
   if (type != "delta") {
     if (.is_htest_of_type(x, "t-test")) {
       return(effectsize(x, type = type, verbose = verbose, data = data, ...))
@@ -203,7 +212,7 @@ glass_delta <- function(x, y = NULL, data = NULL,
 
 
   alternative <- .match.alt(alternative)
-  out <- .get_data_2_samples(x, y, data, paired = paired, verbose = verbose, ...)
+  out <- .get_data_2_samples(x, y, data, paired = paired, reference = reference, verbose = verbose, ...)
   x <- out[["x"]]
   y <- out[["y"]]
   paired <- out[["paired"]]
@@ -308,7 +317,7 @@ glass_delta <- function(x, y = NULL, data = NULL,
     paired, pooled_sd, mu, ci, ci_method, alternative, adjust,
     approximate = FALSE
   )
-  return(out)
+  out
 }
 
 #' @keywords internal

R/common_language.R

@@ -41,7 +41,7 @@
 #'
 #' Where \eqn{U_1}, \eqn{U_2}, and *Overlap* are agnostic to the direction of
 #' the difference between the groups, \eqn{U_3} and probability of superiority
-#' are not.
+#' are not (this can be controlled with the `reference` argument).
 #'
 #' The parametric version of these effects assumes normality of both populations
 #' and homoscedasticity. If those are not met, the non parametric versions
@@ -111,6 +111,7 @@
 #' @aliases cles
 p_superiority <- function(x, y = NULL, data = NULL,
                           mu = 0, paired = FALSE, parametric = TRUE,
+                          reference = NULL,
                           ci = 0.95, alternative = "two.sided",
                           verbose = TRUE, ...) {
   if (.is_htest_of_type(x, "(t-test|Wilcoxon)", "t-test or a Wilcoxon-test")) {
@@ -120,7 +121,7 @@ p_superiority <- function(x, y = NULL, data = NULL,
   }
 
   data <- .get_data_2_samples(x, y, data,
-    paired = paired,
+    paired = paired, reference = reference,
     allow_ordered = !parametric,
     verbose = verbose, ...
   )
@@ -244,6 +245,7 @@ cohens_u2 <- function(x, y = NULL, data = NULL,
 #' @rdname p_superiority
 cohens_u3 <- function(x, y = NULL, data = NULL,
                       mu = 0, parametric = TRUE,
+                      reference = NULL,
                       ci = 0.95, alternative = "two.sided", iterations = 200,
                       verbose = TRUE, ...) {
   if (.is_htest_of_type(x, "(t-test|Wilcoxon)", "t-test or a Wilcoxon-test")) {
@@ -254,7 +256,7 @@ cohens_u3 <- function(x, y = NULL, data = NULL,
 
 
   data <- .get_data_2_samples(x, y, data,
-    allow_ordered = !parametric,
+    allow_ordered = !parametric, reference = reference,
     verbose = verbose, ...
   )
   x <- data[["x"]]
@@ -455,13 +457,13 @@ wmw_odds <- function(x, y = NULL, data = NULL,
 
       out$CI <- ci
 
-      R <- boot::boot(
+      res <- boot::boot(
         data = d,
         statistic = est,
         R = iterations
       )
 
-      bCI <- boot::boot.ci(R, conf = ci, type = "perc")[["percent"]]
+      bCI <- boot::boot.ci(res, conf = ci, type = "perc")[["percent"]]
       bCI <- utils::tail(as.vector(bCI), 2)
       out$CI_low <- bCI[1]
       out$CI_high <- bCI[2]
@@ -478,5 +480,5 @@ wmw_odds <- function(x, y = NULL, data = NULL,
       approximate = TRUE,
       table_footer = "Non-parametric CLES"
     )
-    return(out)
+    out
   }

R/means_ratio.R

@@ -14,8 +14,8 @@
 #'
 #' @details
 #' The Means Ratio ranges from 0 to \eqn{\infty}, with values smaller than 1
-#' indicating that the second mean is larger than the first, values larger than
-#' 1 indicating that the second mean is smaller than the first, and values of 1
+#' indicating that the mean of the reference group is larger, values larger than
+#' 1 indicating that the mean of the reference group is smaller, and values of 1
 #' indicating that the means are equal.
 #'
 #' # Confidence (Compatibility) Intervals (CIs)
@@ -63,15 +63,16 @@
 #' @export
 means_ratio <- function(x, y = NULL, data = NULL,
                         paired = FALSE, adjust = TRUE, log = FALSE,
+                        reference = NULL,
                         ci = 0.95, alternative = "two.sided",
                         verbose = TRUE, ...) {
   alternative <- .match.alt(alternative)
 
   ## Prep data
   out <- .get_data_2_samples(
     x = x, y = y, data = data,
+    paired = paired, reference = reference,
     verbose = verbose,
-    paired = paired,
     ...
   )
   x <- out[["x"]]
@@ -104,14 +105,14 @@ means_ratio <- function(x, y = NULL, data = NULL,
 
     # Calc log RR
     log_val <- .logrom_calc(
-      paired = TRUE,
       m1 = m1,
       sd1 = sd1,
       m2 = m2,
       sd2 = sd2,
       n1 = n,
       r = r,
-      adjust = adjust
+      adjust = adjust,
+      paired = TRUE
     )
   } else {
     ## ------------------------ 2-sample case -------------------------
@@ -121,14 +122,14 @@ means_ratio <- function(x, y = NULL, data = NULL,
 
     # Calc log RR
     log_val <- .logrom_calc(
-      paired = FALSE,
       m1 = m1,
       sd1 = sd1,
       n1 = n1,
       m2 = m2,
       sd2 = sd2,
       n2 = n2,
-      adjust = adjust
+      adjust = adjust,
+      paired = FALSE
     )
   }
 
@@ -175,44 +176,44 @@ means_ratio <- function(x, y = NULL, data = NULL,
     mu = 0,
     approximate = TRUE
   )
-  return(out)
+  out
 }
 
 
 #' @keywords internal
-.logrom_calc <- function(paired = FALSE,
-                         m1,
+.logrom_calc <- function(m1,
                          sd1,
                          n1,
                          m2,
                          sd2,
                          n2 = n1,
                          r = NULL,
-                         adjust = TRUE) {
+                         adjust = TRUE,
+                         paired = FALSE) {
   if (isTRUE(paired)) {
-    yi <- log(m1 / m2)
-    vi <-
+    y_i <- log(m1 / m2)
+    v_i <-
       sd1^2 / (n1 * m1^2) +
       sd2^2 / (n1 * m2^2) -
       2 * r * sd1 * sd2 / (m1 * m2 * n1)
   } else {
-    yi <- log(m1 / m2)
+    y_i <- log(m1 / m2)
     ### large sample approximation to the sampling variance (does not assume homoscedasticity)
-    vi <- sd1^2 / (n1 * m1^2) + sd2^2 / (n2 * m2^2)
+    v_i <- sd1^2 / (n1 * m1^2) + sd2^2 / (n2 * m2^2)
   }
 
 
   if (isTRUE(adjust)) {
     J <- 0.5 * (sd1^2 / (n1 * m1^2) - sd2^2 / (n2 * m2^2))
-    yi <- yi + J
+    y_i <- y_i + J
 
     Jvar <- 0.5 * (sd1^4 / (n1^2 * m1^4) - sd2^4 / (n2^2 * m2^4))
-    vi <- vi + Jvar
+    v_i <- v_i + Jvar
   }
 
 
   list(
-    log_rom = yi,
-    var_rom = vi
+    log_rom = y_i,
+    var_rom = v_i
   )
 }

R/rank_diff.R

@@ -120,6 +120,7 @@
 #' @export
 rank_biserial <- function(x, y = NULL, data = NULL,
                           mu = 0, paired = FALSE,
+                          reference = NULL,
                           ci = 0.95, alternative = "two.sided",
                           verbose = TRUE, ...) {
   alternative <- .match.alt(alternative)
@@ -131,6 +132,7 @@ rank_biserial <- function(x, y = NULL, data = NULL,
   ## Prep data
   out <- .get_data_2_samples(x, y, data,
     paired = paired,
+    reference = reference,
     allow_ordered = TRUE,
     verbose = verbose, ...
   )
@@ -201,19 +203,21 @@ rank_biserial <- function(x, y = NULL, data = NULL,
   attr(out, "ci_method") <- ci_method
   attr(out, "approximate") <- FALSE
   attr(out, "alternative") <- alternative
-  return(out)
+  out
 }
 
 #' @export
 #' @rdname rank_biserial
 cliffs_delta <- function(x, y = NULL, data = NULL,
                          mu = 0,
+                         reference = NULL,
                          ci = 0.95, alternative = "two.sided",
                          verbose = TRUE, ...) {
   cl <- match.call()
   data <- .get_data_2_samples(x, y, data,
     verbose = verbose,
     allow_ordered = TRUE,
+    reference = reference,
     ...
   )
   x <- data$x
@@ -255,5 +259,5 @@ cliffs_delta <- function(x, y = NULL, data = NULL,
 
   u_ <- U1 / S
   f_ <- U2 / S
-  return(u_ - f_)
+  u_ - f_
 }

R/repeated_measures_d.R

@@ -157,6 +157,7 @@ repeated_measures_d <- function(x, y,
                                 data = NULL,
                                 mu = 0, method = c("rm", "av", "z", "b", "d", "r"),
                                 adjust = TRUE,
+                                reference = NULL,
                                 ci = 0.95, alternative = "two.sided",
                                 verbose = TRUE, ...) {
   method <- match.arg(method)
@@ -165,7 +166,11 @@ repeated_measures_d <- function(x, y,
   }
 
   alternative <- .match.alt(alternative)
-  data <- .get_data_paired(x, y, data = data, method = method, verbose = verbose, ...)
+  data <- .get_data_paired(x, y,
+    data = data, method = method,
+    reference = reference,
+    verbose = verbose, ...
+  )
 
   if (method %in% c("d", "r")) {
     values <- .replication_d(data, mu = mu, method = method)
@@ -217,7 +222,7 @@ repeated_measures_d <- function(x, y,
     mu, ci, ci_method, alternative,
     approximate = FALSE
   )
-  return(out)
+  out
 }
 
 #' @rdname repeated_measures_d