Add effect size (Cohen's d) to estimate_contrasts

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",
@@ -22,7 +22,12 @@ Authors@R:
              family = "Patil",
              role = "aut",
              email = "[email protected]",
-             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 = "[email protected]",
+             comment = c(ORCID = "0000-0003-4315-6788", Twitter = "@rempsyc")))
 Maintainer: Dominique Makowski <[email protected]>
 Description: Implements a general interface for model-based estimations
     for a wide variety of models (see list of supported models using the
@@ -63,12 +68,13 @@ Suggests:
     rstanarm,
     rtdists,
     see (>= 0.7.4),
+    bootES,
     testthat
 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)

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

R/estimate_contrasts.R

@@ -11,9 +11,32 @@
 #'   "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
 #'
+#' @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)
@@ -83,6 +106,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,6 +161,49 @@ estimate_contrasts <- function(model,
   contrasts$contrast <- NULL
   contrasts <- cbind(level_cols, contrasts)
 
+  # 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)
+
+    } else if (effectsize == "bootES") {
+      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) {
+        contrast[x, ]
+        })
+
+      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"))
+
+      contrasts <- cbind(contrasts, eff)
+  }
 
   # Table formatting
   attr(contrasts, "table_title") <- c("Marginal Contrasts Analysis", "blue")