Add effect size (Cohen's d) to estimate_contrasts

DESCRIPTION

@@ -22,7 +22,12 @@ Authors@R:
              family = "Patil",
              role = "aut",
              email = "[email protected]",
-             comment = c(ORCID = "0000-0003-1995-6531")))
+             comment = c(ORCID = "0000-0003-1995-6531")),
+      person(given = "Rémi",
+             family = "Thériault",
+             role = c("aut"),
+             email = "[email protected]",
+             comment = c(ORCID = "0000-0003-4315-6788")))
 Maintainer: Dominique Makowski <[email protected]>
 Description: Implements a general interface for model-based estimations
     for a wide variety of models, used in the computation of
@@ -67,6 +72,7 @@ Suggests:
     rmarkdown,
     rstanarm,
     rtdists,
+    bootES,
     see (>= 0.8.4),
     testthat (>= 3.2.1)
 VignetteBuilder:

NEWS.md

@@ -23,6 +23,8 @@
 
 - Fixed issues related to updates of other *easystats* packages.
 
+- `estimate_contrasts`: now supports optional standardized effect sizes, one of "none" (default), "emmeans", or "bootES" (#227, @rempsyc).
+
 # modelbased 0.8.6
 
 ## Breaking Changes

R/estimate_contrasts.R

@@ -10,9 +10,49 @@
 #' comparisons. Can be one of `"holm"` (default), `"tukey"`, `"hochberg"`,
 #' `"hommel"`, `"bonferroni"`, `"BH"`, `"BY"`, `"fdr"` or `"none"`. See the
 #' p-value adjustment section in the `emmeans::test` documentation.
+#' @param effectsize Desired measure of standardized effect size, one of "none"
+#' (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`
+#' effect.type` argument of [bootES::bootES] for details.
+#' @param bootstraps The number of bootstrap resamples to perform.
 #'
 #' @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 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). 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.
+#'
+#' `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). 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 <- 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.
+#'
+#' `effectsize = "bootES"` uses bootstrapping (defaults to a low value of
+#' 200) through [bootES::bootES]. Adjust for contrasts, but not for covariates.
+#'
 #' @examplesIf all(insight::check_if_installed(c("lme4", "emmeans", "rstanarm"), quietly = TRUE))
 #' \dontrun{
 #' # Basic usage
@@ -78,6 +118,9 @@ estimate_contrasts <- function(model,
                                ci = 0.95,
                                p_adjust = "holm",
                                method = "pairwise",
+                               effectsize = "none",
+                               bootstraps = 200,
+                               bootES_type = "cohens.d",
                                backend = "emmeans",
                                ...) {
   if (backend == "emmeans") {
@@ -90,7 +133,7 @@ estimate_contrasts <- function(model,
       adjust = p_adjust,
       ...
     )
-    out <- .format_emmeans_contrasts(model, estimated, ci, transform, p_adjust, ...)
+    out <- .format_emmeans_contrasts(model, estimated, ci, transform, p_adjust, effectsize, ...)
     info <- attributes(estimated)
   } else {
     # Marginalmeans ------------------------------------------------------------
@@ -138,7 +181,7 @@ estimate_contrasts <- function(model,
 # Table formatting emmeans ----------------------------------------------------
 
 
-.format_emmeans_contrasts <- function(model, estimated, ci, transform, p_adjust, ...) {
+.format_emmeans_contrasts <- function(model, estimated, ci, transform, p_adjust, effectsize, ...) {
   # Summarize and clean
   if (insight::model_info(model)$is_bayesian) {
     out <- cbind(estimated@grid, bayestestR::describe_posterior(estimated, ci = ci, verbose = FALSE, ...))
@@ -171,6 +214,63 @@ estimate_contrasts <- function(model,
   cbind(level_cols, out)
 }
 
+  # Add standardized effect size
+  if (!effectsize %in% c("none", "emmeans", "marginal", "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("partial_d", "es_CI_low", "es_CI_high")
+    contrasts <- cbind(contrasts, eff)
+
+  } else if (effectsize == "marginal") {
+    # 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") {
+      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)
+      group <- names([email protected]$xlev)
+      contrast <- estimated@misc$con.coef
+
+      contrast <- lapply(seq_len(nrow(contrast)), function(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(dat),
+          R = bootstraps,
+          data.col = resp,
+          group.col = group,
+          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"))
+
+      contrasts <- cbind(contrasts, eff)
+  }
 
 
 # Table formatting marginal effects -------------------------------------------