Contrasts of ordinal slopes when the by filter keeps one level

NEWS.md

@@ -2,6 +2,14 @@
 
 ## Bug fixes
 
+* If the filter in `by` keeps one level, `estimate_contrasts()` with
+  `backend = "marginaleffects"` now returns contrasts of slopes for ordinal
+  models, such as `MASS::polr()` or `ordinal::clm()`. An example is
+  `estimate_contrasts(m, "x", by = "Infl = 'Low'")`. The output compares the
+  slopes of the response categories at that level. Before, calls with
+  `estimate = "typical"` stopped with the error "No contrasts to show", and
+  calls with `estimate = "average"` ignored the filter.
+
 * 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/format.R

@@ -345,23 +345,43 @@ format.marginaleffects_contrasts <- function(
     # If we have, for example, `contrast = c("vs", "am"), by = "gear='5'"`, the
     # `by` column is the one with one unique value only, we thus have to update
     # `by`, and also `contrast` (the latter not(!) for numerics)...
+    all_by <- by
     by <- by[lengths(lapply(dgrid[by], unique)) > 1]
 
     # for slopes, we have our "levels" in the by-variable, because `by` indicates
     # the levels for which slopes are contrasted - thus, we replace contrast with by
+    one_level_slopes <- FALSE
     if (inherits(x, "estimate_slopes")) {
       contrast <- by
+      # models with response categories (e.g., ordinal models) have one slope
+      # per category. When the filter in `by` keeps one level, e.g.
+      # `by = "Infl='Low'"`, the slopes of the categories at that level are
+      # still contrasted. The labels from marginaleffects then hold only the
+      # response categories, e.g. "(Medium) - (Low)". Custom hypotheses, such
+      # as "b1=b2", and group formulas, such as `~pairwise | Infl`, are not
+      # labelled this way, so these still stop with an error below
+      one_level_slopes <- !length(contrast) &&
+        length(all_by) > 0 &&
+        all(all_by %in% colnames(dgrid)) &&
+        all(lengths(lapply(dgrid[all_by], unique)) == 1) &&
+        !.is_custom_comparison(comparison) &&
+        !(inherits(comparison, "formula") &&
+          grepl("|", insight::safe_deparse(comparison), fixed = TRUE)) &&
+        .is_response_level_model(model, x)
+      if (one_level_slopes) {
+        contrast <- all_by
+      }
       by <- NULL
+    } else {
+      # for contrasts, we also filter variables with one unique value, but we
+      # keep numeric variables. When these are hold constant in the data grid,
+      # they are set to their mean value - meaning, they only have one unique
+      # value in the data grid, anyway. so we need to keep them.
+      keep_contrasts <- lengths(lapply(dgrid[contrast], unique)) > 1 |
+        vapply(dgrid[contrast], is.numeric, logical(1)) # nolint
+      contrast <- contrast[keep_contrasts]
     }
 
-    # for contrasts, we also filter variables with one unique value, but we
-    # keep numeric variables. When these are hold constant in the data grid,
-    # they are set to their mean value - meaning, they only have one unique
-    # value in the data grid, anyway. so we need to keep them.
-    keep_contrasts <- lengths(lapply(dgrid[contrast], unique)) > 1 |
-      vapply(dgrid[contrast], is.numeric, logical(1)) # nolint
-    contrast <- contrast[keep_contrasts]
-
     # set to NULL, if all by-values have been removed here
     if (!length(by)) {
       by <- NULL
@@ -373,8 +393,12 @@ format.marginaleffects_contrasts <- function(
       insight::format_error("No contrasts to show. Please adjust `contrast`.")
     }
 
-    # for more than one term, we have comma-separated levels.
-    if (length(contrast) > 1) {
+    # for slopes with a one-level filter in `by`, the two columns hold the
+    # contrasted response categories. For more than one term, we have
+    # comma-separated levels.
+    if (one_level_slopes) {
+      new_colnames <- c("Response1", "Response2")
+    } else if (length(contrast) > 1) {
       # levels may contain the separator char. to be 100% certain we extract
       # levels correctly, we now replace levels with a special "token", and later
       # replace those tokens with the original levels again
@@ -487,8 +511,11 @@ format.marginaleffects_contrasts <- function(
       # all contrast terms that belong to one "contrast group", separated
       # by comma, and each the two "new contrast groups" go into separate
       # columns named "Level1" and "Level2". For one contrast term, we just
-      # need to rename the columns
-      if (length(contrast) == 1) {
+      # need to rename the columns. For slopes with a one-level filter in `by`,
+      # "Level1" and "Level2" get that level, next to the response categories
+      if (one_level_slopes) {
+        params <- .one_level_slope_labels(params, model, dgrid[all_by])
+      } else if (length(contrast) == 1) {
         # we define this object to avoid scoping issues in data_rename(),
         # see https://github.com/easystats/modelbased/issues/401
         rename_columns <- paste0(contrast, 1:2)
@@ -522,13 +549,17 @@ format.marginaleffects_contrasts <- function(
       # level and the contrasted factor level into the same string (e.g.,
       # "setosa 3"). We split those into separate "Response1" and "Response2"
       # columns, so contrasted factor levels and response levels are not
-      # confounded
-      params <- .split_contrast_response_levels(params, model, x)
+      # confounded. Slopes with a one-level filter in `by` already have these
+      # columns, because their labels hold only the response categories
+      if (!one_level_slopes) {
+        params <- .split_contrast_response_levels(params, model, x)
+      }
 
       # we need to update these variables, because these are the new column
       # names for contrasts and focal terms. We use `intersect()` here because
       # "Response1"/"Response2" are only present if the contrasted factor and
-      # response category levels were successfully split above
+      # response category levels were successfully split above, or for slopes
+      # with a one-level filter in `by`
       contrast <- intersect(
         c("Level1", "Response1", "Level2", "Response2"),
         colnames(params)
@@ -594,6 +625,61 @@ format.marginaleffects_contrasts <- function(
     identical(comparison, "ratio")
 }
 
+# TRUE for models where predictions are made separately for each level of the
+# response (categorical, ordinal/cumulative, multinomial or multivariate models)
+
+.is_response_level_model <- function(model, x = NULL) {
+  if (is.null(model)) {
+    return(FALSE)
+  }
+  # prefer the model info already stored as attribute (consistent with how
+  # other format-methods retrieve model information), and only compute it
+  # ourselves as a fallback
+  model_info <- attr(x, "model_info", exact = TRUE)
+  if (is.null(model_info)) {
+    model_info <- .safe(insight::model_info(model, verbose = FALSE))
+  }
+  isTRUE(model_info$is_multinomial) ||
+    isTRUE(model_info$is_categorical) ||
+    isTRUE(model_info$is_ordinal) ||
+    isTRUE(model_info$is_cumulative) ||
+    isTRUE(.safe(insight::is_multivariate(model)))
+}
+
+# for contrasts of slopes of models with response categories, when the filter
+# in `by` keeps one level: `params` has the response categories in "Response1"
+# and "Response2", and "Level1" and "Level2" get the one level of `by`, from
+# the data grid. Several `by` variables are combined with ", ", as in
+# contrasts of more than one term
+
+.one_level_slope_labels <- function(params, model, by_grid) {
+  by_level <- paste(
+    vapply(by_grid, function(i) as.character(unique(i)), character(1)),
+    collapse = ", "
+  )
+  # keep the order of the response levels, as in
+  # `.split_contrast_response_levels()`. Labels that are no response level,
+  # e.g. because marginaleffects removes parentheses from the labels, follow
+  # in the order in which they appear, so that no label becomes `NA`
+  labels <- unique(c(params$Response1, params$Response2))
+  response <- .safe(insight::get_response(model, verbose = FALSE))
+  if (is.factor(response)) {
+    known_levels <- levels(response)
+  } else if (is.null(dim(response))) {
+    known_levels <- unique(as.character(response))
+  } else {
+    known_levels <- colnames(response)
+  }
+  response_levels <- c(intersect(known_levels, labels), setdiff(labels, known_levels))
+  data.frame(
+    Level1 = factor(by_level),
+    Response1 = factor(params$Response1, levels = response_levels),
+    Level2 = factor(by_level),
+    Response2 = factor(params$Response2, levels = response_levels),
+    stringsAsFactors = FALSE
+  )
+}
+
 # since we combine levels from different factors, we have to make
 # sure levels are unique across different terms. If not, paste
 # variable names to levels. We first find the intersection of all
@@ -613,23 +699,10 @@ format.marginaleffects_contrasts <- function(
     return(params)
   }
 
-  # prefer the model info already stored as attribute (consistent with how
-  # other format-methods retrieve model information), and only compute it
-  # ourselves as a fallback
-  model_info <- attr(x, "model_info", exact = TRUE)
-  if (is.null(model_info)) {
-    model_info <- .safe(insight::model_info(model, verbose = FALSE))
-  }
-  is_multivariate <- isTRUE(.safe(insight::is_multivariate(model)))
-  is_response_model <- isTRUE(model_info$is_multinomial) ||
-    isTRUE(model_info$is_categorical) ||
-    isTRUE(model_info$is_ordinal) ||
-    isTRUE(model_info$is_cumulative) ||
-    is_multivariate
-
-  if (!is_response_model) {
+  if (!.is_response_level_model(model, x)) {
     return(params)
   }
+  is_multivariate <- isTRUE(.safe(insight::is_multivariate(model)))
 
   # response levels, used to detect and split off the response-category
   # label from the contrasted factor level. We keep the original order of

tests/testthat/test-ordinal.R

@@ -1338,6 +1338,76 @@ test_that("estimate_means, backend emmeans, ordinal, several by terms on the res
 })
 
 
+test_that("estimate_contrasts, ordinal slopes, by filter that keeps one level", {
+  skip_if_not_installed("ordinal")
+  d <- MASS::housing
+  set.seed(1)
+  d$x <- rnorm(nrow(d))
+  models <- list(
+    polr = MASS::polr(Sat ~ x * Infl + Type, data = d, Hess = TRUE),
+    clm = ordinal::clm(Sat ~ x * Infl + Type, data = d)
+  )
+  keys <- c("Level1", "Response1", "Level2", "Response2")
+
+  for (m in models) {
+    for (estimate in c("typical", "average")) {
+      out_all <- estimate_contrasts(m, "x", by = "Infl", estimate = estimate)
+
+      # one level: the contrasts between the response categories at that level
+      out <- estimate_contrasts(m, "x", by = "Infl='Low'", estimate = estimate)
+      expected <- out_all[out_all$Level1 == "Low" & out_all$Level2 == "Low", ]
+      expect_identical(nrow(out), 3L)
+      expect_identical(
+        lapply(as.data.frame(out)[keys], as.character),
+        lapply(as.data.frame(expected)[keys], as.character)
+      )
+      expect_equal(out$Difference, expected$Difference, tolerance = 1e-8)
+      expect_equal(out$SE, expected$SE, tolerance = 1e-8)
+
+      # two levels: the matching rows of the unfiltered output
+      out <- estimate_contrasts(m, "x", by = "Infl=c('Low','High')", estimate = estimate)
+      expect_identical(nrow(out), 15L)
+      out_keys <- do.call(paste, lapply(as.data.frame(out)[keys], as.character))
+      all_keys <- do.call(paste, lapply(as.data.frame(out_all)[keys], as.character))
+      matched <- out_all[match(out_keys, all_keys), ]
+      expect_false(anyNA(match(out_keys, all_keys)))
+      expect_equal(out$Difference, matched$Difference, tolerance = 1e-8)
+      expect_equal(out$SE, matched$SE, tolerance = 1e-8)
+    }
+  }
+
+  # response levels with parentheses keep a label
+  d2 <- d
+  levels(d2$Sat) <- c("Low (bad)", "Mid", "High (good)")
+  m <- MASS::polr(Sat ~ x * Infl + Type, data = d2, Hess = TRUE)
+  out <- estimate_contrasts(m, "x", by = "Infl='Low'")
+  expect_identical(nrow(out), 3L)
+  expect_false(anyNA(out$Response1))
+  expect_false(anyNA(out$Response2))
+  expect_length(unique(c(as.character(out$Response1), as.character(out$Response2))), 3L)
+
+  # custom hypotheses and group formulas are not labelled for one level
+  expect_error(
+    estimate_contrasts(models$polr, "x", by = "Infl='Low'", comparison = "b1=b2"),
+    "No contrasts to show",
+    fixed = TRUE
+  )
+  expect_error(
+    estimate_contrasts(models$polr, "x", by = "Infl='Low'", comparison = ~ pairwise | Infl),
+    "No contrasts to show",
+    fixed = TRUE
+  )
+
+  # models without response categories still need two slopes
+  mt <- transform(mtcars, g = factor(am))
+  expect_error(
+    estimate_contrasts(lm(mpg ~ wt * g, data = mt), "wt", by = "g='0'"),
+    "The filter in `by` leaves fewer than two slopes to compare",
+    fixed = TRUE
+  )
+})
+
+
 # Wald statistic of probability contrasts of a polr model, by the delta method.
 # For `polr(Sat ~ Infl + Type)`, P(Sat <= k) = plogis(zeta_k - eta), with eta
 # the linear predictor. The category probabilities are differences of these