Filter the response categories of ordinal models named in by (emmeans backend)

DESCRIPTION

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

NEWS.md

@@ -59,6 +59,18 @@
 
 ## Bug fixes
 
+* `estimate_means()`, `estimate_contrasts()` and `estimate_slopes()` now
+  filter the response categories of ordinal models, for example from
+  `MASS::polr()`, `ordinal::clm()` or `ordinal::clmm()`. The filter is a `by`
+  term that names the response with values, for example
+  `by = c("Infl", "Sat = 'Low'")`. The output then has only the rows of these
+  response categories. The filter works with `backend = "emmeans"` and
+  `predict = "prob"` (the default). Before, means and contrasts of factors
+  ignored the filter, and slopes and contrasts of slopes stopped with an
+  unclear error. If `predict` has another value, these functions now stop with
+  a clear error. If a value is not a response category, they also stop with an
+  error.
+
 * For non-linear models, the footer of predictions or contrasts on the
   link-scale no longer says "FALSE-transformed" or "TRUE-transformed". This
   footer appears, for example, with `estimate_link()`, or with

R/estimate_means.R

@@ -19,7 +19,13 @@
 #' [`insight::get_datagrid()`] to learn more about how to create data grids for
 #' predictors of interest. For `estimate_means()`, `by` defaults to `"auto"`,
 #' which automatically selects the first focal predictor found in the model. If
-#' `by = NULL`, the grand mean of the response is predicted.
+#' `by = NULL`, the grand mean of the response is predicted. For frequentist
+#' ordinal models (for example from `MASS::polr()`, `ordinal::clm()` or
+#' `ordinal::clmm()`), `by` can also name the response with values, for example
+#' `by = c("x", "y = c('low', 'high')")`. `estimate_means()`,
+#' `estimate_contrasts()` and `estimate_slopes()` then keep only the rows of
+#' these response categories. This filter works with `backend = "emmeans"` and
+#' `predict = "prob"` (the default).
 #' @param predict Is passed to the `type` argument in `emmeans::emmeans()` (when
 #' `backend = "emmeans"`) or in `marginaleffects::avg_predictions()` (when
 #' `backend = "marginaleffects"`). Valid options for `predict` are:

R/get_emcontrasts.R

@@ -70,6 +70,7 @@ get_emcontrasts <- function(
     } else {
       fun_args$mode <- ordinal$mode
     }
+    .check_ordinal_response_filter(my_args$response_filter, ordinal$mode)
   }
 
   # add dots
@@ -160,6 +161,7 @@ get_emcontrasts <- function(
   attr(out, "transform") <- TRUE
   attr(out, "keep_iterations") <- keep_iterations
   attr(out, "ordinal_mode") <- ordinal$mode
+  attr(out, "response_filter") <- my_args$response_filter
 
   out
 }
@@ -252,7 +254,12 @@ get_emcontrasts <- function(
     }
     # the response categories or thresholds are renamed first, because the
     # response can have the name of an output column (e.g. "response")
-    out <- .clean_pseudo_factor_emmeans_ordinal(out, model, ordinal_mode)
+    out <- .clean_pseudo_factor_emmeans_ordinal(
+      out,
+      model,
+      ordinal_mode,
+      attributes(estimated)$response_filter
+    )
     out <- .clean_names_frequentist(out, predict, m_info)
   }
   out$null <- NULL # introduced in emmeans 1.6.1 (#115)

R/get_emmeans.R

@@ -92,6 +92,7 @@ get_emmeans <- function(
     } else {
       fun_args$type <- predict
     }
+    .check_ordinal_response_filter(my_args$response_filter, ordinal_mode)
   }
 
   # add dots
@@ -128,6 +129,7 @@ get_emmeans <- function(
   attr(estimated, "transform") <- TRUE
   attr(estimated, "keep_iterations") <- keep_iterations
   attr(estimated, "ordinal_mode") <- ordinal_mode
+  attr(estimated, "response_filter") <- my_args$response_filter
 
   estimated
 }
@@ -318,13 +320,24 @@ get_emmeans <- function(
   if (!is.null(estimate_name) && estimate_name %in% colnames(means)) {
     colnames(means)[colnames(means) == estimate_name] <- .emmeans_ordinal_estimate_names[[ordinal_mode]]
   }
-  .clean_pseudo_factor_emmeans_ordinal(means, model, ordinal_mode)
+  .clean_pseudo_factor_emmeans_ordinal(
+    means,
+    model,
+    ordinal_mode,
+    attributes(x)$response_filter
+  )
 }
 
 
 # for ordinal models, renames the pseudo-factor column (response categories
-# or thresholds) to `Response` or `Threshold`
-.clean_pseudo_factor_emmeans_ordinal <- function(means, model, ordinal_mode) {
+# or thresholds) to `Response` or `Threshold`. With `response_filter`, keeps
+# only the rows of the chosen response categories
+.clean_pseudo_factor_emmeans_ordinal <- function(
+  means,
+  model,
+  ordinal_mode,
+  response_filter = NULL
+) {
   if (is.null(ordinal_mode)) {
     return(means)
   }
@@ -360,6 +373,9 @@ get_emmeans <- function(
         )
       }
     }
+    if (ordinal_mode == "prob") {
+      means <- .filter_ordinal_response_rows(means, pseudo_factor, response_filter)
+    }
     colnames(means)[colnames(means) == pseudo_factor] <- new_name
   }
   means
@@ -412,8 +428,115 @@ get_emmeans <- function(
   }
 }
 
+# For ordinal models, a `by` term can name the response, to keep only some
+# response categories (e.g. `"Sat='Low'"`). emmeans ignores `at` for the
+# response categories with `mode = "prob"`, so we remove the term from `by`,
+# and filter the formatted rows on the response categories. Returns `args`
+# with the term removed from `by`, and with `response_filter`: the name of
+# the response, and the chosen categories (`NULL` when the term has no values)
+.extract_ordinal_response_filter <- function(model, args, data) {
+  if (!is.character(args$by) || !length(args$by)) {
+    return(args)
+  }
+  response <- insight::find_response(model)
+  if (length(response) != 1) {
+    return(args)
+  }
+  is_response <- .grep_cleaned_by_vars(args$by) == response
+  if (
+    !any(is_response) ||
+      !response %in% colnames(data) ||
+      !.is_emmeans_ordinal(model)
+  ) {
+    return(args)
+  }
+  # terms with values, e.g. "Sat='Low'" or "rating=c(1, 3)". Numeric values
+  # are matched with the names of the response categories as character. The
+  # categories of all terms are kept, and a term without values adds none
+  value_terms <- args$by[is_response & grepl("=", args$by, fixed = TRUE)]
+  values <- NULL
+  if (length(value_terms)) {
+    values <- unique(unlist(lapply(value_terms, function(term) {
+      grid <- tryCatch(
+        suppressWarnings(insight::get_datagrid(data[response], by = term)),
+        error = function(e) NULL
+      )
+      term_values <- if (is.null(grid)) NA else as.character(grid[[response]])
+      if (!length(term_values) || anyNA(term_values)) {
+        insight::format_error(paste0(
+          "Could not read the response categories in `",
+          term,
+          "`. Please name the categories in quotes, e.g. `",
+          response,
+          " = c('a', 'b')`."
+        ))
+      }
+      term_values
+    })))
+    # the categories of the model, which can differ from the factor levels in
+    # the data (e.g. `clm()` drops unused levels)
+    response_levels <- model$y.levels
+    if (is.null(response_levels)) {
+      response_levels <- model$lev
+    }
+    if (is.null(response_levels)) {
+      response_levels <- levels(droplevels(as.factor(data[[response]])))
+    }
+    unknown <- setdiff(values, response_levels)
+    if (length(unknown)) {
+      insight::format_error(paste0(
+        "The response `",
+        response,
+        "` has no ",
+        if (length(unknown) == 1) "category " else "categories ",
+        datawizard::text_concatenate(unknown, enclose = "`"),
+        ". Its categories are ",
+        datawizard::text_concatenate(response_levels, enclose = "`"),
+        "."
+      ))
+    }
+  }
+  args$by <- args$by[!is_response]
+  if (!length(args$by)) {
+    args$by <- NULL
+  }
+  args$response_filter <- list(response = response, values = values)
+  args
+}
+
+
+# the response categories are only in the output for `predict = "prob"`, so
+# we cannot use the response in `by` for the other modes. `ordinal_mode` is
+# the mode the user asked for, not the mode that we use internally for slopes
+.check_ordinal_response_filter <- function(response_filter, ordinal_mode) {
+  if (!is.null(response_filter) && !identical(ordinal_mode, "prob")) {
+    insight::format_error(paste0(
+      "The response `",
+      response_filter$response,
+      "` can only be used in `by` with `predict = \"prob\"` (the default), which returns the probability of each response category."
+    ))
+  }
+}
+
+
+# keeps the rows of the response categories that a `by` filter chose (see
+# `.extract_ordinal_response_filter()`). `column` holds the response categories
+.filter_ordinal_response_rows <- function(x, column, response_filter) {
+  values <- response_filter$values
+  if (is.null(values) || !column %in% colnames(x)) {
+    return(x)
+  }
+  x <- x[as.character(x[[column]]) %in% values, , drop = FALSE]
+  rownames(x) <- NULL
+  x
+}
+
+
 #' @keywords internal
 .process_emmeans_arguments <- function(model, args, data, ...) {
+  # ordinal models: a `by` term on the response filters the response categories
+  args <- .extract_ordinal_response_filter(model, args, data)
+
   # Create the data_matrix
   # ---------------------------
   # data <- insight::get_data(model, verbose = FALSE)

R/get_emtrends.R

@@ -64,6 +64,7 @@ get_emtrends <- function(
     } else {
       fun_args <- .emtrends_ordinal_arguments(fun_args, ordinal)
     }
+    .check_ordinal_response_filter(my_args$response_filter, ordinal_mode)
   }
 
   # Run emtrends
@@ -89,6 +90,7 @@ get_emtrends <- function(
   attr(estimated, "transform") <- TRUE
   attr(estimated, "keep_iterations") <- keep_iterations
   attr(estimated, "ordinal_mode") <- ordinal_mode
+  attr(estimated, "response_filter") <- my_args$response_filter
 
   estimated
 }
@@ -249,7 +251,12 @@ get_emtrends <- function(
   # them right before the slopes
   ordinal_mode <- attributes(estimated)$ordinal_mode
   if (!is.null(ordinal_mode)) {
-    trends <- .clean_pseudo_factor_emmeans_ordinal(trends, model, ordinal_mode)
+    trends <- .clean_pseudo_factor_emmeans_ordinal(
+      trends,
+      model,
+      ordinal_mode,
+      attributes(estimated)$response_filter
+    )
     pseudo_column <- intersect(c("Response", "Threshold"), colnames(trends))
     if (length(pseudo_column)) {
       trends <- datawizard::data_relocate(trends, pseudo_column, before = "Slope")