From 31f1596798275c460eb0e3743fa8e63fd08d0efc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?R=C3=A9mi=20Th=C3=A9riault?=
 <13123390+rempsyc@users.noreply.github.com>
Date: Sun, 7 May 2023 20:03:44 -0400
Subject: [PATCH 1/8] add effect size (cohen's d) to estimate_contrasts

---
 DESCRIPTION                              |  2 +-
 R/estimate_contrasts.R                   | 25 ++++++++++++++++++++++++
 tests/testthat/test-estimate_contrasts.R | 22 ++++++++++-----------
 3 files changed, 37 insertions(+), 12 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index fed248f4..c1acc50b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -68,7 +68,7 @@ VignetteBuilder:
     knitr
 Encoding: UTF-8
 Language: en-US
-RoxygenNote: 7.2.3.9000
+RoxygenNote: 7.2.3
 Config/testthat/edition: 3
 Config/testthat/parallel: true
 Roxygen: list(markdown = TRUE)
diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index f5aed15b..fe3aa59b 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -135,6 +135,31 @@ estimate_contrasts <- function(model,
   contrasts$contrast <- NULL
   contrasts <- cbind(level_cols, contrasts)
 
+  # Add effect size (Cohen's d)
+  if (is.null(contrast) && is.null(fixed) && is.null(at)) {
+    dat <- insight::get_data(model)
+    resp <- insight::find_response(model)
+
+    if (is.numeric(dat[[resp]])) {
+      dat <- datawizard::data_select(dat, c(resp, info$contrast))
+
+      list.dat <- lapply(seq_len(nrow(contrasts)), function(i) {
+        log.vec <- which(dat[[info$contrast]] == unlist(info$misc$orig.grid)[[i]])
+        dat.temp <- datawizard::data_filter(dat, log.vec)
+      })
+      list.dat <- stats::setNames(list.dat, unique(unlist(level_cols)))
+
+      eff <- lapply(seq_len(nrow(contrasts)), function(i) {
+        effectsize::cohens_d(x = list.dat[[contrasts$Level1[i]]][[resp]],
+                             y = list.dat[[contrasts$Level2[i]]][[resp]],
+                             verbose = FALSE)
+      })
+
+      eff <- do.call(rbind, eff)
+      names(eff)[-1] <- paste0("Cohens_d_", names(eff)[-1])
+      contrasts <- cbind(contrasts, eff)
+    }
+  }
 
   # Table formatting
   attr(contrasts, "table_title") <- c("Marginal Contrasts Analysis", "blue")
diff --git a/tests/testthat/test-estimate_contrasts.R b/tests/testthat/test-estimate_contrasts.R
index 9a2db2b0..dd505878 100644
--- a/tests/testthat/test-estimate_contrasts.R
+++ b/tests/testthat/test-estimate_contrasts.R
@@ -8,7 +8,7 @@ test_that("estimate_contrasts - Frequentist", {
   model <- lm(Sepal.Width ~ Species, data = dat)
 
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 9))
+  expect_equal(dim(estim), c(3, 13))
 
   estim <- suppressMessages(estimate_contrasts(model, at = "Species=c('versicolor', 'virginica')"))
   expect_equal(dim(estim), c(1, 9))
@@ -21,16 +21,16 @@ test_that("estimate_contrasts - Frequentist", {
   model <- lm(Sepal.Width ~ Species * fac, data = dat)
 
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 9))
+  expect_equal(dim(estim), c(3, 13))
   estim <- suppressMessages(estimate_contrasts(model, levels = "Species"))
-  expect_equal(dim(estim), c(3, 9))
+  expect_equal(dim(estim), c(3, 13))
   estim <- suppressMessages(estimate_contrasts(model, fixed = "fac"))
   expect_equal(dim(estim), c(3, 10))
 
   # One factor and one continuous
   model <- lm(Sepal.Width ~ Species * Petal.Width, data = iris)
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 9))
+  expect_equal(dim(estim), c(3, 13))
   estim <- suppressMessages(estimate_contrasts(model, fixed = "Petal.Width"))
   expect_equal(dim(estim), c(3, 10))
   estim <- suppressMessages(estimate_contrasts(model, at = "Petal.Width", length = 4))
@@ -84,7 +84,7 @@ test_that("estimate_contrasts - Frequentist", {
 
   model <- lme4::lmer(Sepal.Width ~ Species + (1 | Petal.Length_factor), data = data)
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 9))
+  expect_equal(dim(estim), c(3, 13))
 
 
   # GLM - binomial
@@ -107,7 +107,7 @@ test_that("estimate_contrasts - Frequentist", {
   model <- glm(counts ~ treatment, data = dat, family = poisson())
 
   estim <- suppressMessages(estimate_contrasts(model, transform = "response"))
-  expect_equal(dim(estim), c(3, 9))
+  expect_equal(dim(estim), c(3, 13))
 })
 
 
@@ -145,7 +145,7 @@ test_that("estimate_contrasts - Bayesian", {
     )
   )
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 7))
+  expect_equal(dim(estim), c(3, 11))
   estim <- suppressMessages(estimate_contrasts(model, fixed = "Petal.Width"))
   expect_equal(dim(estim), c(3, 8))
   estim <- suppressMessages(estimate_contrasts(model, at = "Petal.Width", length = 4))
@@ -161,14 +161,14 @@ test_that("estimate_contrasts - Bayesian", {
   ))
 
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 7))
+  expect_equal(dim(estim), c(3, 11))
   estim <- suppressMessages(estimate_contrasts(model, transform = "response"))
-  expect_equal(dim(estim), c(3, 7))
+  expect_equal(dim(estim), c(3, 11))
 
   estim <- suppressWarnings(suppressMessages(estimate_contrasts(model, test = "bf")))
-  expect_equal(dim(estim), c(3, 6))
+  expect_equal(dim(estim), c(3, 10))
   estim <- suppressWarnings(suppressMessages(estimate_contrasts(model, transform = "response", test = "bf")))
-  expect_equal(dim(estim), c(3, 6))
+  expect_equal(dim(estim), c(3, 10))
 })
 
 

From d439676627e8d6476ce8799b5a86e55b0672b186 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?R=C3=A9mi=20Th=C3=A9riault?=
 <13123390+rempsyc@users.noreply.github.com>
Date: Sat, 1 Jul 2023 16:25:35 -0400
Subject: [PATCH 2/8] `estimate_contrasts`: now supports optional standardized
 effect sizes, one of "none" (default), "emmeans", or "bootES"

---
 DESCRIPTION                              |  8 ++-
 NEWS.md                                  |  2 +
 R/estimate_contrasts.R                   | 66 +++++++++++++++++-------
 man/estimate_contrasts.Rd                | 11 ++++
 man/modelbased-package.Rd                |  1 +
 tests/testthat/test-estimate_contrasts.R | 22 ++++----
 6 files changed, 79 insertions(+), 31 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index c1acc50b..5d328649 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -22,7 +22,12 @@ Authors@R:
              family = "Patil",
              role = "aut",
              email = "patilindrajeet.science@gmail.com",
-             comment = c(ORCID = "0000-0003-1995-6531", Twitter = "@patilindrajeets")))
+             comment = c(ORCID = "0000-0003-1995-6531", Twitter = "@patilindrajeets")),
+      person(given = "Rémi", 
+             family = "Thériault",
+             role = c("aut"), 
+             email = "remi.theriault@mail.mcgill.ca",
+             comment = c(ORCID = "0000-0003-4315-6788", Twitter = "@rempsyc")))
 Maintainer: Dominique Makowski <dom.makowski@gmail.com>
 Description: Implements a general interface for model-based estimations
     for a wide variety of models (see list of supported models using the
@@ -63,6 +68,7 @@ Suggests:
     rstanarm,
     rtdists,
     see (>= 0.7.4),
+    bootES,
     testthat
 VignetteBuilder:
     knitr
diff --git a/NEWS.md b/NEWS.md
index f18b73e7..7f94f96a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # modelbased (development version)
 
+- `estimate_contrasts`: now supports optional standardized effect sizes, one of "none" (default), "emmeans", or "bootES" (#227, @rempsyc).
+
 # modelbased 0.8.6
 
 ## Breaking Changes
diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index fe3aa59b..353fd7cc 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -11,6 +11,13 @@
 #'   "bonferroni", "BH", "BY", "fdr" or "none". See the p-value adjustment
 #'   section in the `emmeans::test` documentation.
 #' @param adjust Deprecated in favour of `p_adjust`.
+#' @param effectsize Desired measure of standardized effect size, one of "none"
+#' (default), "emmeans", "bootES".
+#' @param bootES_type Specifies the type of effect-size measure to
+#' estimate when using `effectsize = "bootES"`. One of `c("unstandardized",
+#' "cohens.d", "hedges.g", "cohens.d.sigma", "r", "akp.robust.d")`. See`
+#' effect.type` argument of [bootES::bootES] for details.
+#' @param bootstraps The number of bootstrap resamples to perform.
 #'
 #' @inherit estimate_slopes details
 #'
@@ -83,6 +90,9 @@ estimate_contrasts <- function(model,
                                p_adjust = "holm",
                                method = "pairwise",
                                adjust = NULL,
+                               effectsize = "none",
+                               bootstraps = 200,
+                               bootES_type = "cohens.d",
                                ...) {
   # Deprecation
   if (!is.null(adjust)) {
@@ -135,30 +145,48 @@ estimate_contrasts <- function(model,
   contrasts$contrast <- NULL
   contrasts <- cbind(level_cols, contrasts)
 
-  # Add effect size (Cohen's d)
-  if (is.null(contrast) && is.null(fixed) && is.null(at)) {
-    dat <- insight::get_data(model)
-    resp <- insight::find_response(model)
+  # Add standardized effect size
+  if (!effectsize %in% c("none", "emmeans", "bootES")) {
+    message("Unsupported effect size '", effectsize, "', returning none.")
+    }
+
+  if (effectsize == "emmeans") {
+    eff <- emmeans::eff_size(
+      estimated, sigma = stats::sigma(model),
+      edf = stats::df.residual(model), method = "identity")
+    eff <- as.data.frame(eff)
+    eff <- eff[c(2, 5:6)]
+    names(eff) <- c("effect_size", "es_CI_low", "es_CI_high")
+    contrasts <- cbind(contrasts, eff)
 
-    if (is.numeric(dat[[resp]])) {
-      dat <- datawizard::data_select(dat, c(resp, info$contrast))
+    } else if (effectsize == "bootES") {
+      insight::check_if_installed("bootES")
+      dat <- insight::get_data(model)
+      resp <- insight::find_response(model)
+      group <- names(estimated@model.info$xlev)
+      contrast <- estimated@misc$con.coef
 
-      list.dat <- lapply(seq_len(nrow(contrasts)), function(i) {
-        log.vec <- which(dat[[info$contrast]] == unlist(info$misc$orig.grid)[[i]])
-        dat.temp <- datawizard::data_filter(dat, log.vec)
-      })
-      list.dat <- stats::setNames(list.dat, unique(unlist(level_cols)))
+      contrast <- lapply(seq_len(nrow(contrast)), function(x) {
+        contrast[x, ]
+        })
 
-      eff <- lapply(seq_len(nrow(contrasts)), function(i) {
-        effectsize::cohens_d(x = list.dat[[contrasts$Level1[i]]][[resp]],
-                             y = list.dat[[contrasts$Level2[i]]][[resp]],
-                             verbose = FALSE)
-      })
+      es.lists <- lapply(contrast, function(x) {
+        y <- bootES::bootES(
+          data = stats::na.omit(insight::get_data(model)),
+          R = bootstraps,
+          data.col = insight::find_response(model),
+          group.col = insight::find_predictors(model)[[1]],
+          contrast = x,
+          effect.type = bootES_type
+          )
+        y <- as.data.frame(summary(y))})
+
+      eff <- do.call(rbind, es.lists)
+      eff <- eff[1:3]
+      names(eff) <- c(bootES_type, paste0(bootES_type, "_CI_low"),
+                      paste0(bootES_type, "es_CI_high"))
 
-      eff <- do.call(rbind, eff)
-      names(eff)[-1] <- paste0("Cohens_d_", names(eff)[-1])
       contrasts <- cbind(contrasts, eff)
-    }
   }
 
   # Table formatting
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index 250bf016..4abb4624 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -14,6 +14,9 @@ estimate_contrasts(
   p_adjust = "holm",
   method = "pairwise",
   adjust = NULL,
+  effectsize = "none",
+  bootstraps = 200,
+  bootES_type = "cohens.d",
   ...
 )
 }
@@ -53,6 +56,14 @@ section in the \code{emmeans::test} documentation.}
 
 \item{adjust}{Deprecated in favour of \code{p_adjust}.}
 
+\item{effectsize}{Desired measure of standardized effect size, one of "none"
+(default), "emmeans", "bootES".}
+
+\item{bootstraps}{The number of bootstrap resamples to perform.}
+
+\item{bootES_type}{Specifies the type of effect-size measure to
+estimate when using \code{effectsize = "bootES"}. One of \code{c("unstandardized", "cohens.d", "hedges.g", "cohens.d.sigma", "r", "akp.robust.d")}. See\code{ effect.type} argument of \link[bootES:bootES]{bootES::bootES} for details.}
+
 \item{...}{Other arguments passed for instance to \code{\link[insight:get_datagrid]{insight::get_datagrid()}}.}
 }
 \value{
diff --git a/man/modelbased-package.Rd b/man/modelbased-package.Rd
index 6ad16b23..a38d0a76 100644
--- a/man/modelbased-package.Rd
+++ b/man/modelbased-package.Rd
@@ -29,6 +29,7 @@ Authors:
   \item Daniel Lüdecke \email{d.luedecke@uke.de} (\href{https://orcid.org/0000-0002-8895-3206}{ORCID}) (@strengejacke)
   \item Mattan S. Ben-Shachar \email{matanshm@post.bgu.ac.il} (\href{https://orcid.org/0000-0002-4287-4801}{ORCID}) (@mattansb)
   \item Indrajeet Patil \email{patilindrajeet.science@gmail.com} (\href{https://orcid.org/0000-0003-1995-6531}{ORCID}) (@patilindrajeets)
+  \item Rémi Thériault \email{remi.theriault@mail.mcgill.ca} (\href{https://orcid.org/0000-0003-4315-6788}{ORCID}) (@rempsyc)
 }
 
 }
diff --git a/tests/testthat/test-estimate_contrasts.R b/tests/testthat/test-estimate_contrasts.R
index dd505878..9a2db2b0 100644
--- a/tests/testthat/test-estimate_contrasts.R
+++ b/tests/testthat/test-estimate_contrasts.R
@@ -8,7 +8,7 @@ test_that("estimate_contrasts - Frequentist", {
   model <- lm(Sepal.Width ~ Species, data = dat)
 
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 13))
+  expect_equal(dim(estim), c(3, 9))
 
   estim <- suppressMessages(estimate_contrasts(model, at = "Species=c('versicolor', 'virginica')"))
   expect_equal(dim(estim), c(1, 9))
@@ -21,16 +21,16 @@ test_that("estimate_contrasts - Frequentist", {
   model <- lm(Sepal.Width ~ Species * fac, data = dat)
 
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 13))
+  expect_equal(dim(estim), c(3, 9))
   estim <- suppressMessages(estimate_contrasts(model, levels = "Species"))
-  expect_equal(dim(estim), c(3, 13))
+  expect_equal(dim(estim), c(3, 9))
   estim <- suppressMessages(estimate_contrasts(model, fixed = "fac"))
   expect_equal(dim(estim), c(3, 10))
 
   # One factor and one continuous
   model <- lm(Sepal.Width ~ Species * Petal.Width, data = iris)
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 13))
+  expect_equal(dim(estim), c(3, 9))
   estim <- suppressMessages(estimate_contrasts(model, fixed = "Petal.Width"))
   expect_equal(dim(estim), c(3, 10))
   estim <- suppressMessages(estimate_contrasts(model, at = "Petal.Width", length = 4))
@@ -84,7 +84,7 @@ test_that("estimate_contrasts - Frequentist", {
 
   model <- lme4::lmer(Sepal.Width ~ Species + (1 | Petal.Length_factor), data = data)
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 13))
+  expect_equal(dim(estim), c(3, 9))
 
 
   # GLM - binomial
@@ -107,7 +107,7 @@ test_that("estimate_contrasts - Frequentist", {
   model <- glm(counts ~ treatment, data = dat, family = poisson())
 
   estim <- suppressMessages(estimate_contrasts(model, transform = "response"))
-  expect_equal(dim(estim), c(3, 13))
+  expect_equal(dim(estim), c(3, 9))
 })
 
 
@@ -145,7 +145,7 @@ test_that("estimate_contrasts - Bayesian", {
     )
   )
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 11))
+  expect_equal(dim(estim), c(3, 7))
   estim <- suppressMessages(estimate_contrasts(model, fixed = "Petal.Width"))
   expect_equal(dim(estim), c(3, 8))
   estim <- suppressMessages(estimate_contrasts(model, at = "Petal.Width", length = 4))
@@ -161,14 +161,14 @@ test_that("estimate_contrasts - Bayesian", {
   ))
 
   estim <- suppressMessages(estimate_contrasts(model))
-  expect_equal(dim(estim), c(3, 11))
+  expect_equal(dim(estim), c(3, 7))
   estim <- suppressMessages(estimate_contrasts(model, transform = "response"))
-  expect_equal(dim(estim), c(3, 11))
+  expect_equal(dim(estim), c(3, 7))
 
   estim <- suppressWarnings(suppressMessages(estimate_contrasts(model, test = "bf")))
-  expect_equal(dim(estim), c(3, 10))
+  expect_equal(dim(estim), c(3, 6))
   estim <- suppressWarnings(suppressMessages(estimate_contrasts(model, transform = "response", test = "bf")))
-  expect_equal(dim(estim), c(3, 10))
+  expect_equal(dim(estim), c(3, 6))
 })
 
 

From eadb7987f6bcb3dd0b6c96efabaf750ad0cd2733 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?R=C3=A9mi=20Th=C3=A9riault?=
 <13123390+rempsyc@users.noreply.github.com>
Date: Sat, 1 Jul 2023 21:25:43 -0400
Subject: [PATCH 3/8] update docs [skip ci]

---
 DESCRIPTION               |  2 +-
 R/estimate_contrasts.R    | 16 ++++++++++++++++
 man/estimate_contrasts.Rd | 18 ++++++++++++++++++
 3 files changed, 35 insertions(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 5d328649..340aef2d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.8.6.3
+Version: 0.8.6.4
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index 353fd7cc..2b9a7f3d 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -21,6 +21,22 @@
 #'
 #' @inherit estimate_slopes details
 #'
+#' @section Effect Size: By default, `estimate_contrasts` reports no standardized effect size
+#' on purpose. Should one request one, some things are to keep in mind. As the
+#' authors of `emmeans` write, "There is substantial disagreement among
+#' practitioners on what is the appropriate sigma to use in computing effect
+#' sizes; or, indeed, whether any effect-size measure is appropriate for some
+#' situations. The user is completely responsible for specifying appropriate
+#' parameters (or for failing to do so)."
+#'
+#' In particular, effect size methods `"emmeans"` and `"bootES"` do not correct
+#' for covariates in the model, so should probably only be used when there is
+#' just one categorical predictor (with however many levels). If there are
+#' multiple predictors or any covariates, it is important to re-compute sigma
+#' adding back in the response variance associated with the variables that
+#' aren't part of the contrast (or else the Cohen's *d* scale does not really
+#' make sense).
+#'
 #' @examplesIf require("emmeans", quietly = TRUE)
 #' # Basic usage
 #' model <- lm(Sepal.Width ~ Species, data = iris)
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index 4abb4624..b88f2318 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -130,6 +130,24 @@ these levels (using \code{\link[=estimate_contrasts]{estimate_contrasts()}}). Fi
 the effect of x averaged over all conditions, or instead within each
 condition (\code{using [estimate_slopes]}).
 }
+\section{Effect Size}{
+ By default, \code{estimate_contrasts} reports no standardized effect size
+on purpose. Should one request one, some things are to keep in mind. As the
+authors of \code{emmeans} write, "There is substantial disagreement among
+practitioners on what is the appropriate sigma to use in computing effect
+sizes; or, indeed, whether any effect-size measure is appropriate for some
+situations. The user is completely responsible for specifying appropriate
+parameters (or for failing to do so)."
+
+In particular, effect size methods \code{"emmeans"} and \code{"bootES"} do not correct
+for covariates in the model, so should probably only be used when there is
+just one categorical predictor (with however many levels). If there are
+multiple predictors or any covariates, it is important to re-compute sigma
+adding back in the response variance associated with the variables that
+aren't part of the contrast (or else the Cohen's \emph{d} scale does not really
+make sense).
+}
+
 \examples{
 \dontshow{if (require("emmeans", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
 # Basic usage

From 7d0c5bc07792ef96d8449d1fd2bb2cce112f9edf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?R=C3=A9mi=20Th=C3=A9riault?=
 <13123390+rempsyc@users.noreply.github.com>
Date: Sun, 2 Jul 2023 14:24:58 -0400
Subject: [PATCH 4/8] add method "marginal"

---
 R/estimate_contrasts.R    | 63 +++++++++++++++++++++++++++------------
 man/estimate_contrasts.Rd | 43 ++++++++++++++++----------
 2 files changed, 72 insertions(+), 34 deletions(-)

diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index 2b9a7f3d..e008c7a7 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -12,7 +12,7 @@
 #'   section in the `emmeans::test` documentation.
 #' @param adjust Deprecated in favour of `p_adjust`.
 #' @param effectsize Desired measure of standardized effect size, one of "none"
-#' (default), "emmeans", "bootES".
+#' (default), "emmeans", "marginal", or "bootES".
 #' @param bootES_type Specifies the type of effect-size measure to
 #' estimate when using `effectsize = "bootES"`. One of `c("unstandardized",
 #' "cohens.d", "hedges.g", "cohens.d.sigma", "r", "akp.robust.d")`. See`
@@ -21,21 +21,34 @@
 #'
 #' @inherit estimate_slopes details
 #'
-#' @section Effect Size: By default, `estimate_contrasts` reports no standardized effect size
-#' on purpose. Should one request one, some things are to keep in mind. As the
-#' authors of `emmeans` write, "There is substantial disagreement among
-#' practitioners on what is the appropriate sigma to use in computing effect
-#' sizes; or, indeed, whether any effect-size measure is appropriate for some
-#' situations. The user is completely responsible for specifying appropriate
-#' parameters (or for failing to do so)."
+#' @section Effect Size: By default, `estimate_contrasts` reports no
+#' standardized effect size on purpose. Should one request one, some things
+#' are to keep in mind. As the authors of `emmeans` write, "There is
+#' substantial disagreement among practitioners on what is the appropriate
+#' sigma to use in computing effect sizes; or, indeed, whether any effect-size
+#' measure is appropriate for some situations. The user is completely
+#' responsible for specifying appropriate parameters (or for failing to do
+#' so)."
 #'
-#' In particular, effect size methods `"emmeans"` and `"bootES"` do not correct
+#' In particular, effect size method `"bootES"` does not correct
 #' for covariates in the model, so should probably only be used when there is
-#' just one categorical predictor (with however many levels). If there are
-#' multiple predictors or any covariates, it is important to re-compute sigma
-#' adding back in the response variance associated with the variables that
-#' aren't part of the contrast (or else the Cohen's *d* scale does not really
-#' make sense).
+#' just one categorical predictor (with however many levels). Some believe that
+#' if there are multiple predictors or any covariates, it is important to
+#' re-compute sigma adding back in the response variance associated with the
+#' variables that aren't part of the contrast.
+#'
+#' Note also that standardizing the components of a composite variable will
+#' probably lead to nonsense results.
+#'
+#' `effectsize = "emmeans"` uses [emmeans::eff_size] with
+#' `sigma = stats::sigma(model)`, `edf = stats::df.residual(model)` and
+#' `method = "identity")`.
+#'
+#' `effectsize = "marginal"` uses the following formula to compute effect
+#' size: `d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)`.
+#'
+#' `effectsize = "bootES"` uses bootstrapping (defaults to a low value of
+#' 200) through [bootES::bootES]. Does not adjust for covariates.
 #'
 #' @examplesIf require("emmeans", quietly = TRUE)
 #' # Basic usage
@@ -162,7 +175,7 @@ estimate_contrasts <- function(model,
   contrasts <- cbind(level_cols, contrasts)
 
   # Add standardized effect size
-  if (!effectsize %in% c("none", "emmeans", "bootES")) {
+  if (!effectsize %in% c("none", "emmeans", "marginal", "bootES")) {
     message("Unsupported effect size '", effectsize, "', returning none.")
     }
 
@@ -175,7 +188,17 @@ estimate_contrasts <- function(model,
     names(eff) <- c("effect_size", "es_CI_low", "es_CI_high")
     contrasts <- cbind(contrasts, eff)
 
+  } else if (effectsize == "marginal") {
+    # d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)
+    R2_cov <- summary(model)$r.squared
+    d_adj <- contrasts$t * contrasts$SE / sigma(model) * sqrt(1 - R2_cov)
+    contrasts <- cbind(contrasts, marginal_d = d_adj)
+
     } else if (effectsize == "bootES") {
+      if (bootstraps < 500) {
+        message("Number of bootstraps probably too low. Consider increasing it.")
+      }
+
       insight::check_if_installed("bootES")
       dat <- insight::get_data(model)
       resp <- insight::find_response(model)
@@ -183,15 +206,17 @@ estimate_contrasts <- function(model,
       contrast <- estimated@misc$con.coef
 
       contrast <- lapply(seq_len(nrow(contrast)), function(x) {
-        contrast[x, ]
+        z <- contrast[x, ]
+        names(z) <- levels(as.factor(dat[[group]]))
+        z
         })
 
       es.lists <- lapply(contrast, function(x) {
         y <- bootES::bootES(
-          data = stats::na.omit(insight::get_data(model)),
+          data = stats::na.omit(dat),
           R = bootstraps,
-          data.col = insight::find_response(model),
-          group.col = insight::find_predictors(model)[[1]],
+          data.col = resp,
+          group.col = group,
           contrast = x,
           effect.type = bootES_type
           )
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index b88f2318..ac08d3e5 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -57,7 +57,7 @@ section in the \code{emmeans::test} documentation.}
 \item{adjust}{Deprecated in favour of \code{p_adjust}.}
 
 \item{effectsize}{Desired measure of standardized effect size, one of "none"
-(default), "emmeans", "bootES".}
+(default), "emmeans", "marginal", or "bootES".}
 
 \item{bootstraps}{The number of bootstrap resamples to perform.}
 
@@ -131,21 +131,34 @@ the effect of x averaged over all conditions, or instead within each
 condition (\code{using [estimate_slopes]}).
 }
 \section{Effect Size}{
- By default, \code{estimate_contrasts} reports no standardized effect size
-on purpose. Should one request one, some things are to keep in mind. As the
-authors of \code{emmeans} write, "There is substantial disagreement among
-practitioners on what is the appropriate sigma to use in computing effect
-sizes; or, indeed, whether any effect-size measure is appropriate for some
-situations. The user is completely responsible for specifying appropriate
-parameters (or for failing to do so)."
-
-In particular, effect size methods \code{"emmeans"} and \code{"bootES"} do not correct
+ By default, \code{estimate_contrasts} reports no
+standardized effect size on purpose. Should one request one, some things
+are to keep in mind. As the authors of \code{emmeans} write, "There is
+substantial disagreement among practitioners on what is the appropriate
+sigma to use in computing effect sizes; or, indeed, whether any effect-size
+measure is appropriate for some situations. The user is completely
+responsible for specifying appropriate parameters (or for failing to do
+so)."
+
+In particular, effect size method \code{"bootES"} does not correct
 for covariates in the model, so should probably only be used when there is
-just one categorical predictor (with however many levels). If there are
-multiple predictors or any covariates, it is important to re-compute sigma
-adding back in the response variance associated with the variables that
-aren't part of the contrast (or else the Cohen's \emph{d} scale does not really
-make sense).
+just one categorical predictor (with however many levels). Some believe that
+if there are multiple predictors or any covariates, it is important to
+re-compute sigma adding back in the response variance associated with the
+variables that aren't part of the contrast.
+
+Note also that standardizing the components of a composite variable will
+probably lead to nonsense results.
+
+\code{effectsize = "emmeans"} uses \link[emmeans:eff_size]{emmeans::eff_size} with
+\code{sigma = stats::sigma(model)}, \code{edf = stats::df.residual(model)} and
+\verb{method = "identity")}.
+
+\code{effectsize = "marginal"} uses the following formula to compute effect
+size: \code{d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)}.
+
+\code{effectsize = "bootES"} uses bootstrapping (defaults to a low value of
+200) through \link[bootES:bootES]{bootES::bootES}. Does not adjust for covariates.
 }
 
 \examples{

From 219935d241eff95fd2767c35e49f3bcfa39c79cf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?R=C3=A9mi=20Th=C3=A9riault?=
 <13123390+rempsyc@users.noreply.github.com>
Date: Sun, 2 Jul 2023 14:40:03 -0400
Subject: [PATCH 5/8] Buff details of methods [skip ci]

---
 R/estimate_contrasts.R    | 16 ++++++++++------
 man/estimate_contrasts.Rd | 16 ++++++++++------
 2 files changed, 20 insertions(+), 12 deletions(-)

diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index e008c7a7..dad2ec3f 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -37,18 +37,22 @@
 #' re-compute sigma adding back in the response variance associated with the
 #' variables that aren't part of the contrast.
 #'
-#' Note also that standardizing the components of a composite variable will
-#' probably lead to nonsense results.
-#'
 #' `effectsize = "emmeans"` uses [emmeans::eff_size] with
 #' `sigma = stats::sigma(model)`, `edf = stats::df.residual(model)` and
-#' `method = "identity")`.
+#' `method = "identity")`. This standardizes using the MSE (sigma). This works
+#' when the contrasts are the only predictors in the model, but not when there
+#' are covariates. The response variance accounted for by the covariates should
+#' not be removed from the SD used to standardize. Otherwise, _d_ will be
+#' overestimated.
 #'
 #' `effectsize = "marginal"` uses the following formula to compute effect
-#' size: `d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)`.
+#' size: `d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)`. This standardized
+#' using the response SD with only the between-groups variance on the focal
+#' factor/contrast removed. This allows for groups to be equated on their
+#' covariates, but creates an appropriate scale for standardizing the response.
 #'
 #' `effectsize = "bootES"` uses bootstrapping (defaults to a low value of
-#' 200) through [bootES::bootES]. Does not adjust for covariates.
+#' 200) through [bootES::bootES]. Adjust for contrasts, but not for covariates.
 #'
 #' @examplesIf require("emmeans", quietly = TRUE)
 #' # Basic usage
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index ac08d3e5..7d666b6d 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -147,18 +147,22 @@ if there are multiple predictors or any covariates, it is important to
 re-compute sigma adding back in the response variance associated with the
 variables that aren't part of the contrast.
 
-Note also that standardizing the components of a composite variable will
-probably lead to nonsense results.
-
 \code{effectsize = "emmeans"} uses \link[emmeans:eff_size]{emmeans::eff_size} with
 \code{sigma = stats::sigma(model)}, \code{edf = stats::df.residual(model)} and
-\verb{method = "identity")}.
+\verb{method = "identity")}. This standardizes using the MSE (sigma). This works
+when the contrasts are the only predictors in the model, but not when there
+are covariates. The response variance accounted for by the covariates should
+not be removed from the SD used to standardize. Otherwise, \emph{d} will be
+overestimated.
 
 \code{effectsize = "marginal"} uses the following formula to compute effect
-size: \code{d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)}.
+size: \code{d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)}. This standardized
+using the response SD with only the between-groups variance on the focal
+factor/contrast removed. This allows for groups to be equated on their
+covariates, but creates an appropriate scale for standardizing the response.
 
 \code{effectsize = "bootES"} uses bootstrapping (defaults to a low value of
-200) through \link[bootES:bootES]{bootES::bootES}. Does not adjust for covariates.
+200) through \link[bootES:bootES]{bootES::bootES}. Adjust for contrasts, but not for covariates.
 }
 
 \examples{

From 719f546ee1468147e61c6eea24895942adc770b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?R=C3=A9mi=20Th=C3=A9riault?=
 <13123390+rempsyc@users.noreply.github.com>
Date: Sat, 5 Aug 2023 14:50:23 -0400
Subject: [PATCH 6/8] update docs [skip ci]

---
 R/estimate_contrasts.R    | 10 +++++-----
 man/estimate_contrasts.Rd | 10 +++++-----
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index dad2ec3f..7ce53b98 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -39,11 +39,11 @@
 #'
 #' `effectsize = "emmeans"` uses [emmeans::eff_size] with
 #' `sigma = stats::sigma(model)`, `edf = stats::df.residual(model)` and
-#' `method = "identity")`. This standardizes using the MSE (sigma). This works
-#' when the contrasts are the only predictors in the model, but not when there
-#' are covariates. The response variance accounted for by the covariates should
-#' not be removed from the SD used to standardize. Otherwise, _d_ will be
-#' overestimated.
+#' `method = "identity")`. This standardizes using the MSE (sigma). Some believe
+#' this works when the contrasts are the only predictors in the model, but not
+#' when there are covariates. The response variance accounted for by the
+#' covariates should not be removed from the SD used to standardize. Otherwise,
+#' _d_ will be overestimated.
 #'
 #' `effectsize = "marginal"` uses the following formula to compute effect
 #' size: `d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)`. This standardized
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index 7d666b6d..52e8f112 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -149,11 +149,11 @@ variables that aren't part of the contrast.
 
 \code{effectsize = "emmeans"} uses \link[emmeans:eff_size]{emmeans::eff_size} with
 \code{sigma = stats::sigma(model)}, \code{edf = stats::df.residual(model)} and
-\verb{method = "identity")}. This standardizes using the MSE (sigma). This works
-when the contrasts are the only predictors in the model, but not when there
-are covariates. The response variance accounted for by the covariates should
-not be removed from the SD used to standardize. Otherwise, \emph{d} will be
-overestimated.
+\verb{method = "identity")}. This standardizes using the MSE (sigma). Some believe
+this works when the contrasts are the only predictors in the model, but not
+when there are covariates. The response variance accounted for by the
+covariates should not be removed from the SD used to standardize. Otherwise,
+\emph{d} will be overestimated.
 
 \code{effectsize = "marginal"} uses the following formula to compute effect
 size: \code{d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)}. This standardized

From 0c22e29fd7bdef9103b340f9d05451a27e28a529 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?R=C3=A9mi=20Th=C3=A9riault?=
 <13123390+rempsyc@users.noreply.github.com>
Date: Sat, 5 Aug 2023 16:47:11 -0400
Subject: [PATCH 7/8] rewrite as "difference * (1- R2)/ sigma", rename as
 "partial_effect_size" [skip ci]

---
 R/estimate_contrasts.R    | 12 +++++++-----
 man/estimate_contrasts.Rd |  2 +-
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index 7ce53b98..04ed00a3 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -46,7 +46,7 @@
 #' _d_ will be overestimated.
 #'
 #' `effectsize = "marginal"` uses the following formula to compute effect
-#' size: `d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)`. This standardized
+#' size: `d_adj <- difference * (1- R2)/ sigma`. This standardized
 #' using the response SD with only the between-groups variance on the focal
 #' factor/contrast removed. This allows for groups to be equated on their
 #' covariates, but creates an appropriate scale for standardizing the response.
@@ -189,13 +189,15 @@ estimate_contrasts <- function(model,
       edf = stats::df.residual(model), method = "identity")
     eff <- as.data.frame(eff)
     eff <- eff[c(2, 5:6)]
-    names(eff) <- c("effect_size", "es_CI_low", "es_CI_high")
+    names(eff) <- c("partial_effect_size", "es_CI_low", "es_CI_high")
     contrasts <- cbind(contrasts, eff)
 
   } else if (effectsize == "marginal") {
-    # d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)
-    R2_cov <- summary(model)$r.squared
-    d_adj <- contrasts$t * contrasts$SE / sigma(model) * sqrt(1 - R2_cov)
+    # Original: d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)
+    # d_adj <- contrasts$t * contrasts$SE / sigma(model) * sqrt(1 - R2)
+    # New: d_adj <- difference * (1- R2)/ sigma
+    R2 <- summary(model)$r.squared
+    d_adj <- contrasts$Difference * (1 - R2) / sigma(model)
     contrasts <- cbind(contrasts, marginal_d = d_adj)
 
     } else if (effectsize == "bootES") {
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index 52e8f112..870699b9 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -156,7 +156,7 @@ covariates should not be removed from the SD used to standardize. Otherwise,
 \emph{d} will be overestimated.
 
 \code{effectsize = "marginal"} uses the following formula to compute effect
-size: \code{d_adj <- t * se_b / sigma * sqrt(1 - R2_cov)}. This standardized
+size: \code{d_adj <- difference * (1- R2)/ sigma}. This standardized
 using the response SD with only the between-groups variance on the focal
 factor/contrast removed. This allows for groups to be equated on their
 covariates, but creates an appropriate scale for standardizing the response.

From a3aad0d521e7ee514df61f764bed5af42aa0ee47 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?R=C3=A9mi=20Th=C3=A9riault?=
 <13123390+rempsyc@users.noreply.github.com>
Date: Sun, 6 Aug 2023 14:13:54 -0400
Subject: [PATCH 8/8] rewrite output column as "partial_d" when using emmeans
 [skip ci]

---
 R/estimate_contrasts.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index 04ed00a3..fe144d42 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -189,7 +189,7 @@ estimate_contrasts <- function(model,
       edf = stats::df.residual(model), method = "identity")
     eff <- as.data.frame(eff)
     eff <- eff[c(2, 5:6)]
-    names(eff) <- c("partial_effect_size", "es_CI_low", "es_CI_high")
+    names(eff) <- c("partial_d", "es_CI_low", "es_CI_high")
     contrasts <- cbind(contrasts, eff)
 
   } else if (effectsize == "marginal") {