Latent-scale effect sizes for ordinal models, and effect sizes on the rows of their own contrasts (emmeans backend)

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,7 +1,35 @@
 # modelbased (devel)
 
+## Changes
+
+* `estimate_contrasts()` with `backend = "emmeans"` and
+  `effectsize = "emmeans"` now returns effect sizes 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 its error distribution:
+  `pi / sqrt(3)` for the logit link and `1` for the probit link. Other values
+  of `predict`, other links and `clm` models with a `scale` or `nominal`
+  formula give an error. `effectsize = "marginal"` and `effectsize = "boot"`
+  still give an error for ordinal models.
+
 ## Bug fixes
 
+* `estimate_contrasts()` with `effectsize = "emmeans"` now puts each effect
+  size on the row of its own contrast. Before, the contrasts were sometimes sorted
+  in another order than in *emmeans*, for example with factor levels that
+  were not in alphabetical order or with `comparison = "poly"`, and then some
+  effect sizes were on the wrong rows. With a `by` variable, the effect-size
+  columns held the values of the `by` variable and of other *emmeans* columns
+  instead of the effect sizes.
+
+* `estimate_contrasts()` with `effectsize = "emmeans"` now computes the
+  confidence limits of the effect sizes at the level of `ci`. Before, they
+  were always 95% limits.
+
+* `estimate_contrasts()` with an `effectsize` and `comparison = "joint"`, or
+  `comparison = "omnibus"` with `backend = "marginaleffects"`, now gives a
+  clear error. Before, the emmeans effect sizes gave an internal error, and
+  `effectsize = "marginal"` returned `NA` effect sizes.
+
 * 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,11 @@
 #' 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 need a Gaussian model without random
+#' effects. Only `"emmeans"` also works for `polr` and `clm` models with
+#' `predict = "latent"` and a logit or probit link, but not for `clm` models
+#' with a `scale` or `nominal` formula. Effect sizes are not possible for joint
+#' or omnibus tests.
 #' @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`
@@ -242,12 +246,26 @@
 #' 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
+#' `method = "identity"`. For Gaussian models, it uses the residual SD as
+#' `sigma` and the residual degrees of freedom as `edf` (from
+#' `insight::get_sigma()` and `insight::get_df()`). 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.
+#' _d_ will be overestimated. For `polr` and `clm` models with
+#' `predict = "latent"`, sigma is the SD of the error distribution of the
+#' latent variable: `pi / sqrt(3)` (the standard logistic distribution) for
+#' the logit link, and `1` (the standard normal distribution) for the probit
+#' link. With a `rescale` argument for *emmeans*, sigma is multiplied by the
+#' absolute value of `rescale[2]`, so the size of the effect sizes does not
+#' change. A negative `rescale[2]` reverses the sign of the contrasts and of
+#' the effect sizes. For these models, sigma is
+#' taken as known, with `edf = Inf`, so the effect sizes are the contrasts
+#' divided by sigma, and their confidence limits are the unadjusted confidence
+#' limits of the contrasts divided by sigma. For all models, the confidence
+#' limits of the effect sizes are not adjusted for multiple comparisons, also
+#' when `p_adjust` is not `"none"`.
 #'
 #' `effectsize = "marginal"` uses the following formula to compute effect
 #' size: `d_adj <- difference * (1- R2)/ sigma`. This standardizes
@@ -561,9 +579,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,16 @@
     ))
   }
 
+  # For ordinal models, `effectsize = "emmeans"` standardizes the contrasts of
+  # the latent variable by the SD of its error distribution
+  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 with `effectsize = "emmeans"`
+  # are the exception, because the latent variable has a known error SD.
+  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 +58,37 @@
   switch(
     effectsize,
     emmeans = {
+      if (latent_ordinal) {
+        sigma <- .ordinal_latent_sigma(model, attributes(estimated)$predict)
+        # `rescale` (passed to emmeans) scales the latent variable, and so its
+        # 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")
+      # the contrasts can be sorted in another order than the emmeans grid
+      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 +145,59 @@
 }
 
 
+# SD of the error distribution of the latent variable of ordinal models, which
+# is the standard logistic (SD pi / sqrt(3)) for the logit link, and the
+# standard normal (SD 1) for the probit link
+.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."
+    ))
+  )
+}
+
+
+# the `rescale` argument for emmeans, which emmeans also accepts partly
+# matched, e.g. `resc = c(0, 2)`
+.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

@@ -256,6 +256,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) {
@@ -276,6 +277,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)
@@ -294,6 +298,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")
@@ -356,5 +362,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,65 @@
+skip_on_cran()
+skip_if_not_installed("emmeans")
+
+test_that("estimate_contrasts, effectsize emmeans, effect sizes follow their contrasts", {
+  # the levels are not in alphabetical order, so the contrast rows are sorted
+  # in another order than the rows of emmeans
+  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)
+  # e.g. versicolor - setosa: 0.930 / sigma(g) = 0.930 / 0.5148 = 1.807
+  expect_equal(out$partial_d, out$Difference / sigma(g), 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)
+  # e.g. cyl 6 - cyl 4 at am 0: -3.775 / sigma(m) = -3.775 / 3.032 = -1.245
+  expect_type(out$partial_d, "double")
+  expect_equal(out$partial_d, out$Difference / sigma(m), tolerance = 1e-6)
+})
+
+
+test_that("estimate_contrasts, effectsize emmeans, polynomial contrasts and ci level", {
+  # with 6 levels, the polynomial contrasts ("linear" to "degree 5") are
+  # sorted in another order 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)
+
+  g <- lm(Sepal.Length ~ Species, data = iris)
+
+  # the effect-size limits are at the level of `ci`
+  e <- emmeans::eff_size(
+    emmeans::emmeans(g, "Species"),
+    sigma = sigma(g),
+    edf = df.residual(g)
+  )
+  e80 <- as.data.frame(confint(e, level = 0.8, adjust = "none"))
+  out80 <- estimate_contrasts(g, "Species", effectsize = "emmeans", backend = "emmeans", ci = 0.8)
+  # emmeans contrasts the levels in the other direction, so the rows are
+  # matched by the absolute effect size, and the half-widths are compared
+  idx <- match(round(abs(out80$partial_d), 6), round(abs(e80$effect.size), 6))
+  expect_false(anyNA(idx))
+  expect_equal(
+    (out80$es_CI_high - out80$es_CI_low) / 2,
+    ((e80$upper.CL - e80$lower.CL) / 2)[idx],
+    tolerance = 1e-6
+  )
+  out95 <- estimate_contrasts(g, "Species", effectsize = "emmeans", backend = "emmeans")
+  expect_true(all(out80$es_CI_high - out80$es_CI_low < out95$es_CI_high - out95$es_CI_low))
+})