Effect sizes for ordinal models (emmeans backend), and effect-size fixes

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.18.0.3
+Version: 0.18.0.4
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",

NEWS.md

@@ -1,5 +1,12 @@
 # modelbased (devel)
 
+## Changes
+
+* `effectsize = "emmeans"` in `estimate_contrasts()` now also works for `polr`
+  and `clm` models with `predict = "latent"` and a logit or probit link. The
+  contrasts of the latent variable are divided by the SD of the latent error,
+  i.e. `pi / sqrt(3)` (logit) or `1` (probit).
+
 ## Bug fixes
 
 * If the filter in `by` keeps one level, `estimate_contrasts()` with
@@ -10,6 +17,16 @@
   `estimate = "typical"` stopped with the error "No contrasts to show", and
   calls with `estimate = "average"` ignored the filter.
 
+* Fixed issue where `effectsize = "emmeans"` in `estimate_contrasts()` returned
+  effect sizes in the wrong rows (e.g., when factor levels were not sorted
+  alphabetically) or in the wrong columns (when `by` was specified).
+
+* The confidence intervals of effect sizes from `effectsize = "emmeans"` now
+  use the level of `ci`. Before, they were always 95% intervals.
+
+* `estimate_contrasts()` now gives an informative error when effect sizes are
+  requested for joint or omnibus tests.
+
 * Filters on the response categories of frequentist ordinal models with
   `backend = "emmeans"`, where `Sat` is the response:
   * `estimate_means()`, `estimate_contrasts()` and `estimate_slopes()` now

R/estimate_contrasts.R

@@ -103,7 +103,9 @@
 #' specified post-process tests. See 'Examples'.
 #' @param effectsize Desired measure of standardized effect size, one of
 #' `"emmeans"`, `"marginal"`, or `"boot"`. Default is `NULL`, i.e. no effect
-#' size will be computed.
+#' size will be computed. Effect sizes are only available for Gaussian models
+#' without random effects, and (only for `"emmeans"`) for ordinal models on the
+#' latent scale (see 'Details').
 #' @param es_type Specifies the type of effect-size measure to estimate when
 #' using `effectsize = "boot"`. One of `"unstandardized"`, `"cohens.d"`,
 #' `"hedges.g"`, `"cohens.d.sigma"`, `"r"`, or `"akp.robust.d"`. See `effect.type`
@@ -249,6 +251,13 @@
 #' covariates should not be removed from the SD used to standardize. Otherwise,
 #' _d_ will be overestimated.
 #'
+#' For `polr` and `clm` models with `predict = "latent"` and a logit or probit
+#' link, `effectsize = "emmeans"` uses the SD of the latent error as sigma, i.e.
+#' `pi / sqrt(3)` (logit) or `1` (probit), with `edf = Inf`. A `rescale`
+#' argument passed to *emmeans* also rescales sigma. `clm` models with a `scale`
+#' or `nominal` formula are not supported. The confidence intervals of the
+#' effect sizes are not adjusted for multiple comparisons.
+#'
 #' `effectsize = "marginal"` uses the following formula to compute effect
 #' size: `d_adj <- difference * (1- R2)/ sigma`. This standardizes
 #' using the response SD with only the between-groups variance on the focal
@@ -561,9 +570,12 @@ estimate_contrasts.default <- function(
       effectsize = effectsize,
       bootstraps = bootstraps,
       bootES_type = es_type,
-      backend = backend
+      backend = backend,
+      ci = ci,
+      rescale = .rescale_argument(...)
     )
   }
+  attr(out, "emmeans_rows") <- NULL
 
   # sanity check - did method return standard errors?
   .check_standard_errors(

R/estimate_contrasts_effectsize.R

@@ -5,11 +5,17 @@
   effectsize,
   bootstraps,
   bootES_type,
-  backend
+  backend,
+  ci = 0.95,
+  rescale = NULL
 ) {
   # Add standardized effect size
   insight::validate_argument(effectsize, c("none", "emmeans", "marginal", "boot"))
 
+  if (effectsize != "none" && isTRUE(attributes(estimated)$joint_test)) {
+    insight::format_error("Effect sizes are not possible for joint or omnibus tests.")
+  }
+
   if (effectsize == "emmeans" && backend != "emmeans") {
     insight::format_error(
       "`effectsize = \"emmeans\"` only possible with `backend = \"emmeans\"`"
@@ -33,11 +39,15 @@
     ))
   }
 
+  # ordinal models are standardized by the SD of the latent error
+  latent_ordinal <- effectsize == "emmeans" && inherits(model, c("polr", "clm"))
+
   # Check if the model's response variable follows a Gaussian (normal)
   # distribution. Effect size calculations implemented in this function are
   # designed for Gaussian models and may not be appropriate or meaningful for
-  # other types of distributions.
-  if (insight::get_family(model)$family != "gaussian") {
+  # other types of distributions. Ordinal models on the latent scale are an
+  # exception for `effectsize = "emmeans"`.
+  if (!latent_ordinal && insight::get_family(model)$family != "gaussian") {
     insight::format_error(paste0(
       "We strongly recommend not using the `effectsize` ",
       "argument with non-Gaussian models."
@@ -47,15 +57,36 @@
   switch(
     effectsize,
     emmeans = {
+      if (latent_ordinal) {
+        sigma <- .ordinal_latent_sigma(model, attributes(estimated)$predict)
+        # `rescale` also rescales the latent error SD
+        if (!is.null(rescale)) {
+          sigma <- sigma * abs(rescale[2])
+        }
+        edf <- Inf
+      } else {
+        sigma <- insight::get_sigma(model, verbose = FALSE)
+        edf <- insight::get_df(model)
+      }
       eff <- emmeans::eff_size(
         estimated,
-        sigma = insight::get_sigma(model, verbose = FALSE),
-        edf = insight::get_df(model),
+        sigma = sigma,
+        edf = edf,
         method = "identity"
       )
-      eff <- as.data.frame(eff)
-      eff <- eff[c(2, 5:6)]
+      eff <- as.data.frame(stats::confint(eff, level = ci, adjust = "none"))
+      # select columns by name, because `by` variables add columns
+      es_col <- intersect(c("effect.size", "estimate"), colnames(eff))[1]
+      ci_low <- intersect(c("lower.CL", "asymp.LCL"), colnames(eff))[1]
+      ci_high <- intersect(c("upper.CL", "asymp.UCL"), colnames(eff))[1]
+      eff <- eff[c(es_col, ci_low, ci_high)]
       names(eff) <- c("partial_d", "es_CI_low", "es_CI_high")
+      # match effect sizes to the contrasts, which may be sorted differently
+      emmeans_rows <- attributes(contrasts_results)$emmeans_rows
+      if (!is.null(emmeans_rows)) {
+        eff <- eff[emmeans_rows, , drop = FALSE]
+        rownames(eff) <- NULL
+      }
       contrasts_results <- cbind(contrasts_results, eff)
     },
     marginal = {
@@ -112,6 +143,56 @@
 }
 
 
+# SD of the latent error of ordinal models: pi / sqrt(3) for logit, 1 for probit
+.ordinal_latent_sigma <- function(model, predict) {
+  if (!identical(predict, "latent")) {
+    insight::format_error(paste0(
+      "Effect sizes for ordinal models are only possible on the latent scale. ",
+      "Please use `predict = \"latent\"`."
+    ))
+  }
+  if (inherits(model, "clm") && !is.null(model$formulas$scale)) {
+    insight::format_error(
+      "Effect sizes are not possible for ordinal models with a `scale` formula."
+    )
+  }
+  # nominal effects make the thresholds differ between groups
+  if (inherits(model, "clm") && !is.null(model$formulas$nominal)) {
+    insight::format_error(
+      "Effect sizes are not possible for ordinal models with a `nominal` formula."
+    )
+  }
+  link <- if (inherits(model, "polr")) model$method else model$link
+  switch(
+    link,
+    logistic = ,
+    logit = pi / sqrt(3),
+    probit = 1,
+    insight::format_error(paste0(
+      "Effect sizes for ordinal models are only possible for the `logit` and ",
+      "`probit` links, but the model uses the `",
+      link,
+      "` link."
+    ))
+  )
+}
+
+
+# `rescale` argument for emmeans from `...`, also when partly matched
+.rescale_argument <- function(...) {
+  dots <- list(...)
+  dot_names <- names(dots)
+  if (is.null(dot_names)) {
+    return(NULL)
+  }
+  i <- which(nzchar(dot_names) & startsWith("rescale", dot_names))
+  if (length(i) == 0) {
+    return(NULL)
+  }
+  dots[[i[1]]]
+}
+
+
 .get_group_variable <- function(estimated, backend) {
   if (backend == "emmeans") {
     names(estimated@model.info$xlev)

R/get_emcontrasts.R

@@ -267,6 +267,7 @@ get_emcontrasts <- function(
 .format_emmeans_contrasts <- function(model, estimated, ci, p_adjust, ...) {
   predict <- attributes(estimated)$predict
   m_info <- insight::model_info(model, response = 1)
+  emmeans_rows <- NULL
 
   # Summarize and clean
   if (m_info$is_bayesian) {
@@ -287,6 +288,9 @@ get_emcontrasts <- function(
       out <- as.data.frame(estimated)
     } else {
       estimated_data <- as.data.frame(estimated)
+      # `merge()` sorts the rows, so we keep the row of each contrast in the
+      # emmeans grid, to match other emmeans results (e.g. effect sizes)
+      estimated_data$.emmeans_row <- seq_len(nrow(estimated_data))
       out <- as.data.frame(merge(
         estimated_data,
         stats::confint(estimated, level = ci, adjust = p_adjust)
@@ -305,6 +309,8 @@ get_emcontrasts <- function(
         out <- out[match(row_keys(estimated_data), row_keys(out)), , drop = FALSE]
         rownames(out) <- NULL
       }
+      emmeans_rows <- out$.emmeans_row
+      out$.emmeans_row <- NULL
     }
     # the response categories or thresholds are renamed first, because the
     # response can have the name of an output column (e.g. "response")
@@ -367,5 +373,9 @@ get_emcontrasts <- function(
   }
 
   # add posterior draws?
-  .add_posterior_draws_emmeans(attributes(estimated), out)
+  out <- .add_posterior_draws_emmeans(attributes(estimated), out)
+  if (!is.null(emmeans_rows)) {
+    attr(out, "emmeans_rows") <- emmeans_rows
+  }
+  out
 }

tests/testthat/test-estimate_contrasts_effectsize_order.R

@@ -0,0 +1,29 @@
+skip_on_cran()
+skip_if_not_installed("emmeans")
+
+test_that("estimate_contrasts, effectsize emmeans, effect sizes follow their contrasts", {
+  # levels not in alphabetical order
+  d <- iris
+  d$Species <- factor(d$Species, levels = c("virginica", "setosa", "versicolor"))
+  g <- lm(Sepal.Length ~ Species, data = d)
+  out <- estimate_contrasts(g, "Species", effectsize = "emmeans", backend = "emmeans")
+  expect_identical(nrow(out), 3L)
+  expect_equal(out$partial_d, out$Difference / sigma(g), tolerance = 1e-6)
+
+  # with 6 levels, polynomial contrasts are sorted differently than in emmeans
+  w <- lm(weight ~ feed, data = chickwts)
+  out <- estimate_contrasts(w, "feed", comparison = "poly", effectsize = "emmeans", backend = "emmeans")
+  expect_identical(nrow(out), 5L)
+  expect_equal(out$partial_d, out$Difference / sigma(w), tolerance = 1e-6)
+})
+
+
+test_that("estimate_contrasts, effectsize emmeans, effect sizes with a by variable", {
+  d <- mtcars
+  d$cyl <- factor(d$cyl)
+  d$am <- factor(d$am)
+  m <- lm(mpg ~ cyl * am, data = d)
+  out <- estimate_contrasts(m, "cyl", by = "am", effectsize = "emmeans", backend = "emmeans")
+  expect_identical(nrow(out), 6L)
+  expect_equal(out$partial_d, out$Difference / sigma(m), tolerance = 1e-6)
+})

tests/testthat/test-ordinal.R

@@ -1986,3 +1986,91 @@ test_that("estimate_contrasts, backend emmeans, ordinal, slopes with a by filter
     c("Difference", "SE")
   )
 })
+
+
+test_that("estimate_contrasts, backend emmeans, ordinal, latent-scale effect sizes", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("ordinal")
+  data(housing, package = "MASS")
+  es_cols <- c("partial_d", "es_CI_low", "es_CI_high")
+  ci_cols <- c("Difference", "CI_low", "CI_high")
+
+  # logit: the latent error SD is pi / sqrt(3)
+  m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE)
+  out <- estimate_contrasts(m, "Infl", predict = "latent", effectsize = "emmeans", backend = "emmeans")
+  expect_identical(nrow(out), 3L)
+  expect_equal(out[es_cols], out[ci_cols] / (pi / sqrt(3)), tolerance = 1e-6, ignore_attr = TRUE)
+
+  # probit: the latent error SD is 1
+  mp <- ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing, link = "probit")
+  out <- estimate_contrasts(mp, "Infl", predict = "latent", effectsize = "emmeans", backend = "emmeans")
+  expect_equal(out[es_cols], out[ci_cols], tolerance = 1e-6, ignore_attr = TRUE)
+
+  # the effect-size limits are at the level of `ci`
+  out <- estimate_contrasts(
+    m,
+    "Infl",
+    predict = "latent",
+    effectsize = "emmeans",
+    backend = "emmeans",
+    ci = 0.8
+  )
+  expect_equal(out[es_cols], out[ci_cols] / (pi / sqrt(3)), tolerance = 1e-6, ignore_attr = TRUE)
+
+  # a partly matched, negative `rescale` also rescales the latent error SD
+  out_r <- estimate_contrasts(
+    m,
+    "Infl",
+    predict = "latent",
+    effectsize = "emmeans",
+    backend = "emmeans",
+    resc = c(1, -3)
+  )
+  expect_equal(out_r$Difference, -3 * out$Difference, tolerance = 1e-6)
+  expect_equal(out_r$partial_d, -out$partial_d, tolerance = 1e-6)
+})
+
+
+test_that("estimate_contrasts, backend emmeans, ordinal, effect size errors", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("ordinal")
+  data(housing, package = "MASS")
+  m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE)
+
+  expect_error(
+    estimate_contrasts(m, "Infl", effectsize = "emmeans", backend = "emmeans"),
+    "latent"
+  )
+  expect_error(
+    estimate_contrasts(m, "Infl", predict = "latent", effectsize = "marginal", backend = "emmeans"),
+    "non-Gaussian"
+  )
+  expect_error(
+    estimate_contrasts(
+      m,
+      "Infl",
+      by = "Type",
+      comparison = "joint",
+      predict = "latent",
+      effectsize = "emmeans",
+      backend = "emmeans"
+    ),
+    "omnibus"
+  )
+
+  mc <- ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing, link = "cloglog")
+  expect_error(
+    estimate_contrasts(mc, "Infl", predict = "latent", effectsize = "emmeans", backend = "emmeans"),
+    "cloglog"
+  )
+  ms <- ordinal::clm(Sat ~ Infl + Type, scale = ~Cont, weights = Freq, data = housing)
+  expect_error(
+    estimate_contrasts(ms, "Infl", predict = "latent", effectsize = "emmeans", backend = "emmeans"),
+    "scale"
+  )
+  mn <- ordinal::clm(Sat ~ Infl, nominal = ~Type, weights = Freq, data = housing)
+  expect_error(
+    estimate_contrasts(mn, "Type", predict = "latent", effectsize = "emmeans", backend = "emmeans"),
+    "nominal"
+  )
+})