New repeated_measures_d() (#618)

This PR introduces the new `repeated_measures_d()` function (#372).
easystats · Nov 7, 2023 · b0caa91 · b0caa91
1 parent 448e39b
commit b0caa91
Show file tree

Hide file tree

Showing 40 changed files with 3,484 additions and 306 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -211,11 +211,13 @@ export(rb_to_common_language)
 export(rb_to_p_superiority)
 export(rb_to_vda)
 export(rb_to_wmw_odds)
+export(repeated_measures_d)
 export(riskratio)
 export(riskratio_to_arr)
 export(riskratio_to_logoddsratio)
 export(riskratio_to_nnt)
 export(riskratio_to_oddsratio)
+export(rm_d)
 export(rules)
 export(sd_pooled)
 export(standardise)

diff --git a/NEWS.md b/NEWS.md
@@ -2,6 +2,11 @@
 
 - This release changes the licensing model of `{effectsize}` to an MIT license.
 
+## New features
+
+- `repeated_measures_d()` to compute standardized mean differences (SMD) for repeated measures data.
+  - Also supported in `effectsize(<t.test(paired = TRUE)>)`
+
 # effectsize 0.8.7
 
 ## New features

diff --git a/R/cohens_d.R b/R/cohens_d.R
@@ -25,7 +25,7 @@
 #'   variance). Else the mean SD from both groups is used instead.
 #' @param paired If `TRUE`, the values of `x` and `y` are considered as paired.
 #'   This produces an effect size that is equivalent to the one-sample effect
-#'   size on `x - y`.
+#'   size on `x - y`. See also [repeated_measures_d()] for more options.
 #' @param ... Arguments passed to or from other methods. When `x` is a formula,
 #'   these can be `subset` and `na.action`.
 #' @inheritParams chisq_to_phi
@@ -48,7 +48,7 @@
 #'   `Glass_delta`) and their CIs (`CI_low` and `CI_high`).
 #'
 #' @family standardized differences
-#' @seealso [sd_pooled()], [t_to_d()], [r_to_d()]
+#' @seealso [rm_d()], [sd_pooled()], [t_to_d()], [r_to_d()]
 #'
 #' @examples
 #' \donttest{
@@ -91,10 +91,12 @@
 #'
 #' # same as:
 #' # cohens_d(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2], paired = TRUE)
+#' # cohens_d(sleep$extra[sleep$group == 1] - sleep$extra[sleep$group == 2])
+#' # rm_d(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2], method = "z", adjust = FALSE)
 #'
 #' # More options:
-#' cohens_d(Pair(extra[group == 1], extra[group == 2]) ~ 1, data = sleep, mu = -1)
-#' hedges_g(Pair(extra[group == 1], extra[group == 2]) ~ 1, data = sleep)
+#' cohens_d(Pair(extra[group == 1], extra[group == 2]) ~ 1, data = sleep, mu = -1, verbose = FALSE)
+#' hedges_g(Pair(extra[group == 1], extra[group == 2]) ~ 1, data = sleep, verbose = FALSE)
 #'
 #'
 #' # Interpretation -----------------------
@@ -208,6 +210,12 @@ glass_delta <- function(x, y = NULL, data = NULL,
   y <- out[["y"]]
   paired <- out[["paired"]]
 
+  if (verbose && paired && !is.null(y)) {
+    insight::format_alert(
+      "For paired samples, 'repeated_measures_d()' provides more options."
+    )
+  }
+
   if (is.null(y)) {
     if (type == "delta") {
       insight::format_error("For Glass' Delta, please provide data from two samples.")
@@ -287,7 +295,7 @@ glass_delta <- function(x, y = NULL, data = NULL,
 
 
   if (type == "g") {
-    J <- exp(lgamma(df / 2) - log(sqrt(df / 2)) - lgamma((df - 1) / 2)) # exact method
+    J <- .J(df)
 
     out[, colnames(out) %in% c("Hedges_g", "CI_low", "CI_high")] <-
       out[, colnames(out) %in% c("Hedges_g", "CI_low", "CI_high")] * J
@@ -300,3 +308,8 @@ glass_delta <- function(x, y = NULL, data = NULL,
   )
   return(out)
 }
+
+#' @keywords internal
+.J <- function(df) {
+  exp(lgamma(df / 2) - log(sqrt(df / 2)) - lgamma((df - 1) / 2)) # exact method
+}
diff --git a/R/common_language.R b/R/common_language.R
@@ -7,6 +7,9 @@
 #' any reported [`stats::t.test()`] or [`stats::wilcox.test()`].
 #'
 #' @inheritParams cohens_d
+#' @param paired If `TRUE`, the values of `x` and `y` are considered as paired.
+#'   This produces an effect size that is equivalent to the one-sample effect
+#'   size on `x - y`.
 #' @param parametric Use parametric estimation (see [cohens_d()]) or
 #'   non-parametric estimation (see [rank_biserial()]). See details.
 #' @param iterations The number of bootstrap replicates for computing confidence
@@ -126,6 +129,10 @@ p_superiority <- function(x, y = NULL, data = NULL,
   paired <- data[["paired"]]
 
   if (parametric) {
+    if (paired) {
+      x <- x - y
+      y <- NULL
+    }
     d <- cohens_d(
       x = x,
       y = y,

diff --git a/R/datasets.R b/R/datasets.R
@@ -28,6 +28,36 @@
 #'
 NULL
 
+#' Jeff Rouder's Example Dataset for Repeated Measures
+#'
+#' A dataset "with 25 people each observing 50 trials in 2 conditions",
+#' published as `effectSizePuzzler.txt` by Jeff Rouder on March 24, 2016
+#' (_http://jeffrouder.blogspot.com/2016/03/the-effect-size-puzzler.html_).
+#' \cr\cr
+#' The data is used in examples and tests of [rm_d()].
+#'
+#' @docType data
+#'
+#' @name rouder2016
+#'
+#' @keywords data
+#'
+#' @format A data frame with 2500 rows and 3 variables:
+#' \describe{
+#'   \item{id}{participant: 1...25}
+#'   \item{cond}{condition: 1,2}
+#'   \item{rt}{response time in seconds}
+#' }
+#'
+#' ```{r}
+#' data("rouder2016")
+#' head(rouder2016, n = 5)
+#' ```
+#' @family effect size datasets
+#'
+NULL
+
+
 #' Results from 2 Screening Tests
 #'
 #' A sample (simulated) dataset, used in tests and some examples.

diff --git a/R/effectsize.R b/R/effectsize.R
@@ -14,6 +14,7 @@
 #'
 #' - For an object of class `htest`, data is extracted via [insight::get_data()], and passed to the relevant function according to:
 #'   - A **t-test** depending on `type`: `"cohens_d"` (default), `"hedges_g"`, or one of `"p_superiority"`, `"u1"`, `"u2"`, `"u3"`, `"overlap"`.
+#'     - For a **Paired t-test**: depending on `type`: `"rm"`, `"av"`, `"b"`, `"d"`, `"z"`.
 #'   - A **Chi-squared tests of independence** or **Fisher's Exact Test**, depending on `type`: `"cramers_v"` (default), `"tschuprows_t"`, `"phi"`, `"cohens_w"`, `"pearsons_c"`, `"cohens_h"`, `"oddsratio"`, `"riskratio"`, `"arr"`, or `"nnt"`.
 #'   - A **Chi-squared tests of goodness-of-fit**, depending on `type`: `"fei"` (default) `"cohens_w"`, `"pearsons_c"`
 #'   - A **One-way ANOVA test**, depending on `type`: `"eta"` (default), `"omega"` or `"epsilon"` -squared, `"f"`, or `"f2"`.
@@ -50,6 +51,13 @@
 #' Tt <- t.test(1:10, y = c(7:20), alternative = "less")
 #' effectsize(Tt)
 #'
+#' Tt <- t.test(
+#'   x = c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30),
+#'   y = c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29),
+#'   paired = TRUE
+#' )
+#' effectsize(Tt, type = "rm_b")
+#'
 #' Aov <- oneway.test(extra ~ group, data = sleep, var.equal = TRUE)
 #' effectsize(Aov)
 #' effectsize(Aov, type = "omega")

diff --git a/R/effectsize.htest.R b/R/effectsize.htest.R
@@ -44,11 +44,11 @@ effectsize.htest <- function(model, type = NULL, verbose = TRUE, ...) {
   dots$alternative <- model$alternative
   dots$ci <- attr(model$conf.int, "conf.level")
   dots$mu <- model$null.value
-  dots$paired <- !grepl("Two", model$method, fixed = TRUE)
+  dots$paired <- grepl("Paired", model$method, fixed = TRUE)
   dots$verbose <- verbose
 
   if (!type %in% c("d", "g")) {
-    .fail_if_approx(approx, "cles")
+    .fail_if_approx(approx, if (startsWith(type, "rm")) "rm_d" else "cles")
   }
 
   if (approx) {
@@ -81,6 +81,11 @@ effectsize.htest <- function(model, type = NULL, verbose = TRUE, ...) {
         d = cohens_d,
         g = hedges_g
       )
+    } else if (dots$paired && startsWith(type, "rm")) {
+      args[c("x", "y")] <- split(args$x, args$y)
+      dots$paired <- args$pooled_sd <- NULL
+      args$method <- gsub("^rm\\_", "", type)
+      f <- rm_d
     } else {
       if (!dots$paired && !args$pooled_sd) {
         insight::format_error("Common language effect size only applicable to Cohen's d with pooled SD.")

diff --git a/R/mahalanobis_D.R b/R/mahalanobis_D.R
@@ -85,7 +85,6 @@ mahalanobis_d <- function(x, y = NULL, data = NULL,
                           pooled_cov = TRUE, mu = 0,
                           ci = 0.95, alternative = "greater",
                           verbose = TRUE, ...) {
-  # TODO add one sample case DV1 + DV2 ~ 1
   # TODO add paired samples case DV1 + DV2 ~ 1 | ID
   alternative <- .match.alt(alternative, FALSE)
   data <- .get_data_multivariate(x, y, data, verbose = verbose, ...)

diff --git a/R/print.effectsize_table.R b/R/print.effectsize_table.R
@@ -161,7 +161,8 @@ print.effectsize_difference <- function(x, digits = 2, append_CLES = NULL, ...)
 
 #' @export
 format.effectsize_difference <- function(x, digits = 2, ...) {
-  caption <- subtitle <- footer <- NULL
+  caption <- subtitle <- NULL
+  footer <- attr(x, "table_footer")
 
   ## Add footer
   mu <- attr(x, "mu")

diff --git a/R/rank_diff.R b/R/rank_diff.R
@@ -6,6 +6,7 @@
 #' to the [Common Language Effect Sizes][cohens_u3]. Pair with any reported
 #' [`stats::wilcox.test()`].
 #'
+#' @inheritParams p_superiority
 #' @inheritParams cohens_d
 #' @param x,y A numeric or ordered vector, or a character name of one in `data`.
 #'   Any missing values (`NA`s) are dropped from the resulting vector. `x` can