diff --git a/DESCRIPTION b/DESCRIPTION
index a9f38f15d..c2cd6668c 100644
--- a/DESCRIPTION
+++ b/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",
diff --git a/NEWS.md b/NEWS.md
index e46983722..a88fe66fa 100644
--- a/NEWS.md
+++ b/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
diff --git a/R/estimate_means.R b/R/estimate_means.R
index b5224de94..73010d172 100644
--- a/R/estimate_means.R
+++ b/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:
diff --git a/R/get_emcontrasts.R b/R/get_emcontrasts.R
index 899b90e2b..a0e301db2 100644
--- a/R/get_emcontrasts.R
+++ b/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)
diff --git a/R/get_emmeans.R b/R/get_emmeans.R
index bf8563017..12cb0f3ac 100644
--- a/R/get_emmeans.R
+++ b/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)
diff --git a/R/get_emtrends.R b/R/get_emtrends.R
index 2f10e41ec..3c568bab0 100644
--- a/R/get_emtrends.R
+++ b/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")
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index 444b80fb0..bf67bb755 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -92,7 +92,13 @@ case, \code{by} can also be list of named elements. See details in
 \code{\link[insight:get_datagrid]{insight::get_datagrid()}} to learn more about how to create data grids for
 predictors of interest. For \code{estimate_means()}, \code{by} defaults to \code{"auto"},
 which automatically selects the first focal predictor found in the model. If
-\code{by = NULL}, the grand mean of the response is predicted.}
+\code{by = NULL}, the grand mean of the response is predicted. For frequentist
+ordinal models (for example from \code{MASS::polr()}, \code{ordinal::clm()} or
+\code{ordinal::clmm()}), \code{by} can also name the response with values, for example
+\code{by = c("x", "y = c('low', 'high')")}. \code{estimate_means()},
+\code{estimate_contrasts()} and \code{estimate_slopes()} then keep only the rows of
+these response categories. This filter works with \code{backend = "emmeans"} and
+\code{predict = "prob"} (the default).}
 
 \item{predict}{Is passed to the \code{type} argument in \code{emmeans::emmeans()} (when
 \code{backend = "emmeans"}) or in \code{marginaleffects::avg_predictions()} (when
diff --git a/man/estimate_means.Rd b/man/estimate_means.Rd
index 61a36e3fc..86da04dcb 100644
--- a/man/estimate_means.Rd
+++ b/man/estimate_means.Rd
@@ -32,7 +32,13 @@ case, \code{by} can also be list of named elements. See details in
 \code{\link[insight:get_datagrid]{insight::get_datagrid()}} to learn more about how to create data grids for
 predictors of interest. For \code{estimate_means()}, \code{by} defaults to \code{"auto"},
 which automatically selects the first focal predictor found in the model. If
-\code{by = NULL}, the grand mean of the response is predicted.}
+\code{by = NULL}, the grand mean of the response is predicted. For frequentist
+ordinal models (for example from \code{MASS::polr()}, \code{ordinal::clm()} or
+\code{ordinal::clmm()}), \code{by} can also name the response with values, for example
+\code{by = c("x", "y = c('low', 'high')")}. \code{estimate_means()},
+\code{estimate_contrasts()} and \code{estimate_slopes()} then keep only the rows of
+these response categories. This filter works with \code{backend = "emmeans"} and
+\code{predict = "prob"} (the default).}
 
 \item{predict}{Is passed to the \code{type} argument in \code{emmeans::emmeans()} (when
 \code{backend = "emmeans"}) or in \code{marginaleffects::avg_predictions()} (when
diff --git a/man/estimate_slopes.Rd b/man/estimate_slopes.Rd
index 72966f0fa..c766719c7 100644
--- a/man/estimate_slopes.Rd
+++ b/man/estimate_slopes.Rd
@@ -46,7 +46,13 @@ case, \code{by} can also be list of named elements. See details in
 \code{\link[insight:get_datagrid]{insight::get_datagrid()}} to learn more about how to create data grids for
 predictors of interest. For \code{estimate_means()}, \code{by} defaults to \code{"auto"},
 which automatically selects the first focal predictor found in the model. If
-\code{by = NULL}, the grand mean of the response is predicted.}
+\code{by = NULL}, the grand mean of the response is predicted. For frequentist
+ordinal models (for example from \code{MASS::polr()}, \code{ordinal::clm()} or
+\code{ordinal::clmm()}), \code{by} can also name the response with values, for example
+\code{by = c("x", "y = c('low', 'high')")}. \code{estimate_means()},
+\code{estimate_contrasts()} and \code{estimate_slopes()} then keep only the rows of
+these response categories. This filter works with \code{backend = "emmeans"} and
+\code{predict = "prob"} (the default).}
 
 \item{predict}{Is passed to the \code{type} argument in \code{emmeans::emmeans()} (when
 \code{backend = "emmeans"}) or in \code{marginaleffects::avg_predictions()} (when
diff --git a/man/get_emmeans.Rd b/man/get_emmeans.Rd
index 2cf912c2e..fe43628cc 100644
--- a/man/get_emmeans.Rd
+++ b/man/get_emmeans.Rd
@@ -113,7 +113,13 @@ case, \code{by} can also be list of named elements. See details in
 \code{\link[insight:get_datagrid]{insight::get_datagrid()}} to learn more about how to create data grids for
 predictors of interest. For \code{estimate_means()}, \code{by} defaults to \code{"auto"},
 which automatically selects the first focal predictor found in the model. If
-\code{by = NULL}, the grand mean of the response is predicted.}
+\code{by = NULL}, the grand mean of the response is predicted. For frequentist
+ordinal models (for example from \code{MASS::polr()}, \code{ordinal::clm()} or
+\code{ordinal::clmm()}), \code{by} can also name the response with values, for example
+\code{by = c("x", "y = c('low', 'high')")}. \code{estimate_means()},
+\code{estimate_contrasts()} and \code{estimate_slopes()} then keep only the rows of
+these response categories. This filter works with \code{backend = "emmeans"} and
+\code{predict = "prob"} (the default).}
 
 \item{predict}{Is passed to the \code{type} argument in \code{emmeans::emmeans()} (when
 \code{backend = "emmeans"}) or in \code{marginaleffects::avg_predictions()} (when
diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R
index 3497382d6..71b5d20af 100644
--- a/tests/testthat/test-ordinal.R
+++ b/tests/testthat/test-ordinal.R
@@ -1084,3 +1084,255 @@ test_that("estimate_slopes, backend emmeans, ordinal glmmTMB", {
   expect_identical(paste(out$Infl, out$Response), paste(out_polr$Infl, out_polr$Response))
   expect_equal(out$Slope, out_polr$Slope, tolerance = 1e-4)
 })
+
+
+# `by` filters on the response categories, emmeans backend --------------------
+
+# the output of a call with a `by` filter on the response categories must be
+# the rows of the same call without that filter, for the chosen categories
+expect_response_filter_rows <- function(out, out_all, response_levels, info = NULL) {
+  expected <- out_all[as.character(out_all$Response) %in% response_levels, , drop = FALSE]
+  expect_identical(nrow(out), nrow(expected), info = info)
+  expect_gt(nrow(out), 0)
+  expect_identical(colnames(out), colnames(out_all), info = info)
+  for (column in colnames(expected)) {
+    if (is.numeric(expected[[column]])) {
+      expect_equal(
+        out[[column]],
+        expected[[column]],
+        tolerance = 1e-10,
+        info = paste(info, column)
+      )
+    } else {
+      expect_identical(
+        as.character(out[[column]]),
+        as.character(expected[[column]]),
+        info = paste(info, column)
+      )
+    }
+  }
+}
+
+
+ordinal_filter_models <- function() {
+  data(housing, package = "MASS", envir = environment())
+  data(wine, package = "ordinal", envir = environment())
+  wine_letters <- wine
+  levels(wine_letters$rating) <- letters[1:5]
+  list(
+    polr = list(
+      model = MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE),
+      factor = "Infl",
+      response = "Sat",
+      one = "Low",
+      two = c("Low", "High")
+    ),
+    clm = list(
+      model = ordinal::clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing),
+      factor = "Infl",
+      response = "Sat",
+      one = "Medium",
+      two = c("Low", "High")
+    ),
+    clmm = list(
+      model = ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine),
+      factor = "temp",
+      response = "rating",
+      one = "2",
+      two = c("1", "3")
+    ),
+    clmm_letters = list(
+      model = ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine_letters),
+      factor = "temp",
+      response = "rating",
+      one = "b",
+      two = c("a", "c")
+    )
+  )
+}
+
+
+response_filter_term <- function(response, response_levels) {
+  if (length(response_levels) == 1) {
+    paste0(response, "='", response_levels, "'")
+  } else {
+    paste0(response, "=c(", toString(paste0("'", response_levels, "'")), ")")
+  }
+}
+
+
+test_that("estimate_means, backend emmeans, ordinal, by filter on the response", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("ordinal")
+  models <- ordinal_filter_models()
+  for (name in names(models)) {
+    fit <- models[[name]]
+    out_all <- estimate_means(fit$model, fit$factor, backend = "emmeans")
+    for (response_levels in list(fit$one, fit$two)) {
+      term <- response_filter_term(fit$response, response_levels)
+      out <- estimate_means(fit$model, c(fit$factor, term), backend = "emmeans")
+      expect_response_filter_rows(out, out_all, response_levels, info = paste(name, term))
+    }
+  }
+})
+
+
+test_that("estimate_contrasts, backend emmeans, ordinal, by filter on the response", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("ordinal")
+  models <- ordinal_filter_models()
+  for (name in names(models)) {
+    fit <- models[[name]]
+    p_adjusts <- if (name %in% c("polr", "clm")) c("none", "holm") else "none"
+    for (p_adjust in p_adjusts) {
+      out_all <- estimate_contrasts(
+        fit$model,
+        fit$factor,
+        backend = "emmeans",
+        p_adjust = p_adjust
+      )
+      for (response_levels in list(fit$one, fit$two)) {
+        term <- response_filter_term(fit$response, response_levels)
+        out <- estimate_contrasts(
+          fit$model,
+          fit$factor,
+          by = term,
+          backend = "emmeans",
+          p_adjust = p_adjust
+        )
+        expect_response_filter_rows(
+          out,
+          out_all,
+          response_levels,
+          info = paste(name, term, p_adjust)
+        )
+      }
+    }
+  }
+})
+
+
+test_that("estimate_means, backend emmeans, ordinal, response named in by without values", {
+  skip_if_not_installed("emmeans")
+  data(housing, package = "MASS")
+  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
+  out_all <- estimate_means(m, "Infl", backend = "emmeans")
+  out <- estimate_means(m, c("Infl", "Sat"), backend = "emmeans")
+  expect_response_filter_rows(out, out_all, levels(housing$Sat))
+  expect_identical(nrow(out), 9L)
+})
+
+
+test_that("estimate_slopes, backend emmeans, ordinal, by filter on the response", {
+  skip_if_not_installed("emmeans")
+  data(housing, package = "MASS")
+  set.seed(1)
+  housing$x <- as.numeric(housing$Cont) + stats::rnorm(nrow(housing))
+  m <- MASS::polr(Sat ~ Infl * x, weights = Freq, data = housing, Hess = TRUE)
+
+  out_all <- estimate_slopes(m, "x", by = "Infl", backend = "emmeans")
+  out <- estimate_slopes(m, "x", by = c("Infl", "Sat='Low'"), backend = "emmeans")
+  expect_response_filter_rows(out, out_all, "Low")
+  expect_identical(nrow(out), 3L)
+
+  # contrasts of slopes
+  out_all <- estimate_contrasts(m, "x", by = "Infl", backend = "emmeans")
+  out <- estimate_contrasts(m, "x", by = c("Infl", "Sat='Low'"), backend = "emmeans")
+  expect_response_filter_rows(out, out_all, "Low")
+  expect_identical(nrow(out), 3L)
+
+  # slopes use `"cum.prob"` internally, but the mode check uses `predict`
+  err <- expect_error(estimate_slopes(
+    m,
+    "x",
+    by = c("Infl", "Sat='Low'"),
+    predict = "cum.prob",
+    backend = "emmeans"
+  ))
+  expect_match(
+    gsub("\\s+", " ", conditionMessage(err)),
+    "The response `Sat` can only be used in `by` with `predict = \"prob\"`",
+    fixed = TRUE
+  )
+})
+
+
+test_that("estimate_means and estimate_contrasts, backend emmeans, ordinal, errors for response filters", {
+  skip_if_not_installed("emmeans")
+  data(housing, package = "MASS")
+  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE)
+  # the messages are wrapped at spaces, so we compare them on one line
+  error_message <- function(expr) {
+    err <- expect_error(expr)
+    gsub("\\s+", " ", conditionMessage(err))
+  }
+
+  # the response categories can only be filtered for `predict = "prob"`
+  modes <- c("cum.prob", "exc.prob", "linear.predictor", "latent", "mean.class")
+  for (mode in modes) {
+    messages <- c(
+      error_message(estimate_means(
+        m,
+        c("Infl", "Sat='Low'"),
+        predict = mode,
+        backend = "emmeans"
+      )),
+      error_message(estimate_contrasts(
+        m,
+        "Infl",
+        by = "Sat='Low'",
+        predict = mode,
+        backend = "emmeans"
+      ))
+    )
+    for (msg in messages) {
+      expect_match(msg, "The response `Sat` can only be used in `by`", fixed = TRUE, info = mode)
+      expect_match(msg, "`predict = \"prob\"`", fixed = TRUE, info = mode)
+    }
+  }
+
+  # values that are no response categories
+  messages <- c(
+    error_message(estimate_means(m, c("Infl", "Sat='Lowest'"), backend = "emmeans")),
+    error_message(estimate_contrasts(m, "Infl", by = "Sat='Lowest'", backend = "emmeans"))
+  )
+  for (msg in messages) {
+    expect_match(msg, "The response `Sat` has no category `Lowest`.", fixed = TRUE)
+  }
+  msg <- error_message(estimate_means(
+    m,
+    c("Infl", "Sat=c('Lowest','Highest','Low')"),
+    backend = "emmeans"
+  ))
+  expect_match(msg, "has no categories `Lowest` and `Highest`.", fixed = TRUE)
+
+  # values that cannot be read as response categories
+  msg <- error_message(estimate_means(m, c("Infl", "Sat=[fivenum]"), backend = "emmeans"))
+  expect_match(msg, "Could not read the response categories in `Sat=[fivenum]`.", fixed = TRUE)
+
+  # a category of the data that the model does not have
+  skip_if_not_installed("ordinal")
+  data(wine, package = "ordinal")
+  wine$rating <- factor(wine$rating, levels = 1:6, ordered = TRUE)
+  m_clm <- ordinal::clm(rating ~ temp, data = wine)
+  msg <- error_message(estimate_means(m_clm, c("temp", "rating='6'"), backend = "emmeans"))
+  expect_match(msg, "The response `rating` has no category `6`.", fixed = TRUE)
+})
+
+
+test_that("estimate_means, backend emmeans, ordinal, several by terms on the response", {
+  skip_if_not_installed("emmeans")
+  data(housing, package = "MASS")
+  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
+  out_all <- estimate_means(m, "Infl", backend = "emmeans")
+
+  # the categories of all terms are kept
+  out <- estimate_means(m, c("Infl", "Sat='Low'", "Sat='High'"), backend = "emmeans")
+  expect_response_filter_rows(out, out_all, c("Low", "High"))
+  expect_identical(nrow(out), 6L)
+
+  # the response without values does not undo the filter
+  out <- estimate_means(m, c("Infl", "Sat", "Sat='Low'"), backend = "emmeans")
+  expect_response_filter_rows(out, out_all, "Low")
+  expect_identical(nrow(out), 3L)
+})