diff --git a/NEWS.md b/NEWS.md
index 77bbbc05c..71e2f42b3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -18,8 +18,15 @@
   meaningful predictions, slopes, and interaction contrasts using modelbased
   post-estimation functions.
 
-* Backend `"emmeans"` now supports ordinal models by passing the `predict`
-  argument to `mode` in the `emmeans()` call.
+* `estimate_means()` with `backend = "emmeans"` now supports ordinal models,
+  e.g. from `MASS::polr()` or `ordinal::clm()`. `predict` accepts the ordinal
+  modes of *emmeans* (`"prob"`, `"cum.prob"`, `"exc.prob"`,
+  `"linear.predictor"`, `"latent"` and `"mean.class"`). For `"prob"`, the
+  output has one row per response category, in a `Response` column. For
+  `"cum.prob"`, `"exc.prob"` and `"linear.predictor"`, the output has one row
+  per threshold, in a `Threshold` column. For ordinal models, the default is
+  now `predict = "prob"`, which returns the same probabilities as the
+  `"marginaleffects"` backend.
 
 ## Bug fixes
 
diff --git a/R/estimate_means.R b/R/estimate_means.R
index c8c067ca6..e4914f684 100644
--- a/R/estimate_means.R
+++ b/R/estimate_means.R
@@ -38,6 +38,12 @@
 #'   `"unlink"`, or `"log"`. If `predict = NULL` (default), the most appropriate
 #'   transformation is selected (which usually is `"response"`). See also
 #'   [this vignette](https://CRAN.R-project.org/package=emmeans/vignettes/transformations.html).
+#'   For ordinal models (e.g. from `MASS::polr()` or `ordinal::clm()`),
+#'   `estimate_means()` also accepts the ordinal modes of *emmeans*:
+#'   `"prob"`, `"cum.prob"`, `"exc.prob"`, `"linear.predictor"`, `"latent"`
+#'   and `"mean.class"` (see `vignette("models", package = "emmeans")`). For
+#'   these models, the default is `"prob"`, which returns one probability per
+#'   response category.
 #'
 #' See also section _Predictions on different scales_.
 #'
diff --git a/R/get_emmeans.R b/R/get_emmeans.R
index 3f93d66fb..920efc5c7 100644
--- a/R/get_emmeans.R
+++ b/R/get_emmeans.R
@@ -62,6 +62,9 @@ get_emmeans <- function(
   # setup arguments
   fun_args <- list(model, specs = my_args$emmeans_specs, at = my_args$emmeans_at)
 
+  # ordinal models: prediction mode, used for column names when formatting
+  ordinal_mode <- NULL
+
   # handle distributional parameters
   if (inherits(model, "brmsfit")) {
     if (identical(predict, "response")) {
@@ -77,8 +80,28 @@ get_emmeans <- function(
     }
   } else {
     dpars <- FALSE
+    # for ordinal models, emmeans ignores `type = "response"`, and we need the
+    # `mode` argument instead. Probabilities are the default.
+    ordinal_model <- .is_emmeans_ordinal(model)
+    if (ordinal_model && identical(predict, "response")) {
+      predict <- "prob"
+    }
     if (predict %in% .emmeans_ordinal_types) {
       fun_args$mode <- predict
+      # some modes add a pseudo-factor to the reference grid (the response
+      # categories, or the thresholds). We add it to `specs`, else emmeans
+      # averages over it
+      if (ordinal_model) {
+        ordinal_mode <- predict
+        pseudo_factor <- .emmeans_ordinal_pseudo_factor(model, predict)
+        if (!is.null(pseudo_factor)) {
+          if (inherits(fun_args$specs, "formula")) {
+            fun_args$specs <- pseudo_factor
+          } else {
+            fun_args$specs <- c(fun_args$specs, pseudo_factor)
+          }
+        }
+      }
     } else {
       fun_args$type <- predict
     }
@@ -117,6 +140,7 @@ get_emmeans <- function(
   attr(estimated, "focal_terms") <- my_args$emmeans_specs
   attr(estimated, "transform") <- TRUE
   attr(estimated, "keep_iterations") <- keep_iterations
+  attr(estimated, "ordinal_mode") <- ordinal_mode
 
   estimated
 }
@@ -136,6 +160,43 @@ get_emmeans <- function(
 )
 
 
+# column names for the estimates of each ordinal mode
+.emmeans_ordinal_estimate_names <- c(
+  latent = "Latent",
+  linear.predictor = "Linear_predictor",
+  cum.prob = "Probability",
+  exc.prob = "Probability",
+  prob = "Probability",
+  mean.class = "Mean_class"
+)
+
+
+# frequentist ordinal models, which support the `mode` argument in emmeans
+.is_emmeans_ordinal <- function(model) {
+  if (inherits(model, "brmsfit")) {
+    return(FALSE)
+  }
+  m_info <- insight::model_info(model, response = 1, verbose = FALSE)
+  isTRUE(m_info$is_ordinal) && !isTRUE(m_info$is_bayesian)
+}
+
+
+# name of the pseudo-factor that emmeans adds to the reference grid for
+# ordinal models: the response categories for `mode = "prob"`, and the
+# thresholds for `mode = "cum.prob"`, `"exc.prob"` and `"linear.predictor"`.
+# emmeans names the response pseudo-factor after the left-hand side of the
+# formula, e.g. `factor(y)`, so we need the response term, not the variable
+.emmeans_ordinal_pseudo_factor <- function(model, mode) {
+  switch(mode,
+    prob = insight::find_terms(model, verbose = FALSE)$response,
+    cum.prob = ,
+    exc.prob = ,
+    linear.predictor = "cut",
+    NULL
+  )
+}
+
+
 #' @keywords internal
 .guess_emmeans_arguments <- function(model, by = NULL, verbose = TRUE, ...) {
   # Gather info
@@ -208,6 +269,7 @@ get_emmeans <- function(
   } else {
     means <- as.data.frame(stats::confint(x, level = ci))
     means$df <- NULL
+    means <- .clean_names_emmeans_ordinal(means, x, model)
     means <- .clean_names_frequentist(means, predict, m_info)
   }
 
@@ -226,6 +288,26 @@ get_emmeans <- function(
 }
 
 
+# for ordinal models, renames the estimate column and the pseudo-factor
+# columns (response categories or thresholds)
+.clean_names_emmeans_ordinal <- function(means, x, model) {
+  ordinal_mode <- attributes(x)$ordinal_mode
+  if (is.null(ordinal_mode)) {
+    return(means)
+  }
+  estimate_name <- x@misc$estName
+  if (!is.null(estimate_name) && estimate_name %in% colnames(means)) {
+    colnames(means)[colnames(means) == estimate_name] <- .emmeans_ordinal_estimate_names[[ordinal_mode]]
+  }
+  pseudo_factor <- .emmeans_ordinal_pseudo_factor(model, ordinal_mode)
+  if (!is.null(pseudo_factor) && pseudo_factor %in% colnames(means)) {
+    new_name <- if (ordinal_mode == "prob") "Response" else "Threshold"
+    colnames(means)[colnames(means) == pseudo_factor] <- new_name
+  }
+  means
+}
+
+
 # adds posterior draws to output for emmeans objects
 .add_posterior_draws_emmeans <- function(info, estimated) {
   # add posterior draws?
diff --git a/R/utils.R b/R/utils.R
index 60d7b8795..6053416b8 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -39,7 +39,10 @@
     "Median",
     "MAP",
     "Coefficient",
-    "Odds_ratio"
+    "Odds_ratio",
+    "Latent",
+    "Linear_predictor",
+    "Mean_class"
   )
   dpars <- insight::find_auxiliary(model, verbose = FALSE)
   if (!is.null(dpars)) {
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index 9c7d69be4..1067ec14c 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -112,6 +112,12 @@ response scale.
 \code{"unlink"}, or \code{"log"}. If \code{predict = NULL} (default), the most appropriate
 transformation is selected (which usually is \code{"response"}). See also
 \href{https://CRAN.R-project.org/package=emmeans/vignettes/transformations.html}{this vignette}.
+For ordinal models (e.g. from \code{MASS::polr()} or \code{ordinal::clm()}),
+\code{estimate_means()} also accepts the ordinal modes of \emph{emmeans}:
+\code{"prob"}, \code{"cum.prob"}, \code{"exc.prob"}, \code{"linear.predictor"}, \code{"latent"}
+and \code{"mean.class"} (see \code{vignette("models", package = "emmeans")}). For
+these models, the default is \code{"prob"}, which returns one probability per
+response category.
 }
 
 See also section \emph{Predictions on different scales}.}
diff --git a/man/estimate_means.Rd b/man/estimate_means.Rd
index ff5e03c9e..b97165164 100644
--- a/man/estimate_means.Rd
+++ b/man/estimate_means.Rd
@@ -52,6 +52,12 @@ response scale.
 \code{"unlink"}, or \code{"log"}. If \code{predict = NULL} (default), the most appropriate
 transformation is selected (which usually is \code{"response"}). See also
 \href{https://CRAN.R-project.org/package=emmeans/vignettes/transformations.html}{this vignette}.
+For ordinal models (e.g. from \code{MASS::polr()} or \code{ordinal::clm()}),
+\code{estimate_means()} also accepts the ordinal modes of \emph{emmeans}:
+\code{"prob"}, \code{"cum.prob"}, \code{"exc.prob"}, \code{"linear.predictor"}, \code{"latent"}
+and \code{"mean.class"} (see \code{vignette("models", package = "emmeans")}). For
+these models, the default is \code{"prob"}, which returns one probability per
+response category.
 }
 
 See also section \emph{Predictions on different scales}.}
diff --git a/man/estimate_slopes.Rd b/man/estimate_slopes.Rd
index 1d3ac0ad8..55b3d97c2 100644
--- a/man/estimate_slopes.Rd
+++ b/man/estimate_slopes.Rd
@@ -66,6 +66,12 @@ response scale.
 \code{"unlink"}, or \code{"log"}. If \code{predict = NULL} (default), the most appropriate
 transformation is selected (which usually is \code{"response"}). See also
 \href{https://CRAN.R-project.org/package=emmeans/vignettes/transformations.html}{this vignette}.
+For ordinal models (e.g. from \code{MASS::polr()} or \code{ordinal::clm()}),
+\code{estimate_means()} also accepts the ordinal modes of \emph{emmeans}:
+\code{"prob"}, \code{"cum.prob"}, \code{"exc.prob"}, \code{"linear.predictor"}, \code{"latent"}
+and \code{"mean.class"} (see \code{vignette("models", package = "emmeans")}). For
+these models, the default is \code{"prob"}, which returns one probability per
+response category.
 }
 
 See also section \emph{Predictions on different scales}.}
diff --git a/man/get_emmeans.Rd b/man/get_emmeans.Rd
index 910d50433..310644bda 100644
--- a/man/get_emmeans.Rd
+++ b/man/get_emmeans.Rd
@@ -133,6 +133,12 @@ response scale.
 \code{"unlink"}, or \code{"log"}. If \code{predict = NULL} (default), the most appropriate
 transformation is selected (which usually is \code{"response"}). See also
 \href{https://CRAN.R-project.org/package=emmeans/vignettes/transformations.html}{this vignette}.
+For ordinal models (e.g. from \code{MASS::polr()} or \code{ordinal::clm()}),
+\code{estimate_means()} also accepts the ordinal modes of \emph{emmeans}:
+\code{"prob"}, \code{"cum.prob"}, \code{"exc.prob"}, \code{"linear.predictor"}, \code{"latent"}
+and \code{"mean.class"} (see \code{vignette("models", package = "emmeans")}). For
+these models, the default is \code{"prob"}, which returns one probability per
+response category.
 }
 
 See also section \emph{Predictions on different scales}.}
diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R
index 7a45254a1..0cebbd25c 100644
--- a/tests/testthat/test-ordinal.R
+++ b/tests/testthat/test-ordinal.R
@@ -32,11 +32,169 @@ test_that("estimate_relation prints ordinal models correctly", {
 })
 
 
-test_that("estimate_means with backend emmeans work for ordinal", {
+# compares probabilities from the emmeans and the marginaleffects backend,
+# matching rows by focal terms and response category
+.compare_ordinal_backends <- function(out_emmeans, out_marginaleffects, by) {
+  merge_by <- c(by, "Response")
+  out_emmeans <- as.data.frame(out_emmeans)[c(merge_by, "Probability")]
+  out_marginaleffects <- as.data.frame(out_marginaleffects)[c(merge_by, "Probability")]
+  for (i in merge_by) {
+    out_emmeans[[i]] <- as.character(out_emmeans[[i]])
+    out_marginaleffects[[i]] <- as.character(out_marginaleffects[[i]])
+  }
+  merge(out_emmeans, out_marginaleffects, by = merge_by)
+}
+
+
+test_that("estimate_means, backend emmeans, ordinal, predict = 'prob'", {
+  skip_if_not_installed("emmeans")
+  data(housing, package = "MASS")
+  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
+
+  # one row per focal level and response category
+  out <- estimate_means(m, "Type", predict = "prob", backend = "emmeans")
+  expect_identical(nrow(out), 12L)
+  expect_setequal(as.character(out$Response), levels(housing$Sat))
+  compared <- .compare_ordinal_backends(out, estimate_means(m, "Type"), "Type")
+  expect_identical(nrow(compared), 12L)
+  expect_equal(compared$Probability.x, compared$Probability.y, tolerance = 1e-6)
+
+  # two focal terms
+  out <- estimate_means(m, c("Type", "Infl"), predict = "prob", backend = "emmeans")
+  expect_identical(nrow(out), 36L)
+  compared <- .compare_ordinal_backends(
+    out,
+    estimate_means(m, c("Type", "Infl")),
+    c("Type", "Infl")
+  )
+  expect_identical(nrow(compared), 36L)
+  expect_equal(compared$Probability.x, compared$Probability.y, tolerance = 1e-6)
+})
+
+
+test_that("estimate_means, backend emmeans, ordinal, default predict", {
+  skip_if_not_installed("emmeans")
+  data(housing, package = "MASS")
+  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
+
+  # default returns probabilities, same as predict = "prob"
+  out_prob <- estimate_means(m, "Type", predict = "prob", backend = "emmeans")
+  out_default <- estimate_means(m, "Type", backend = "emmeans")
+  expect_identical(nrow(out_default), 12L)
+  expect_identical(
+    paste(out_default$Type, out_default$Response),
+    paste(out_prob$Type, out_prob$Response)
+  )
+  expect_equal(out_default$Probability, out_prob$Probability, tolerance = 1e-6)
+
+  # contrasts are not affected by the new default for means
+  out <- estimate_contrasts(m, "Type", backend = "emmeans")
+  expect_identical(
+    paste(out$Level1, out$Level2),
+    c(
+      "Apartment Tower", "Atrium Apartment", "Atrium Tower",
+      "Terrace Apartment", "Terrace Atrium", "Terrace Tower"
+    )
+  )
+  expect_equal(
+    out$Difference,
+    c(-0.5723501, 0.2061636, -0.3661866, -0.5186648, -0.7248283, -1.0910149),
+    tolerance = 1e-6
+  )
+})
+
+
+test_that("estimate_means, backend emmeans, ordinal, mean.class and latent", {
+  skip_if_not_installed("emmeans")
+  data(housing, package = "MASS")
+  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
+
+  modes <- c(mean.class = "Mean_class", latent = "Latent")
+  for (mode in names(modes)) {
+    out <- estimate_means(m, "Type", predict = mode, backend = "emmeans")
+    expect_identical(nrow(out), 4L)
+    expect_true(modes[[mode]] %in% colnames(out))
+    expect_identical(attributes(out)$coef_name, modes[[mode]])
+    expected <- as.data.frame(suppressMessages(
+      emmeans::emmeans(m, "Type", mode = mode)
+    ))
+    estimate_column <- setdiff(
+      colnames(expected),
+      c("Type", "SE", "df", "asymp.LCL", "asymp.UCL")
+    )
+    compared <- merge(
+      data.frame(
+        Type = as.character(out$Type),
+        x = out[[modes[[mode]]]],
+        stringsAsFactors = FALSE
+      ),
+      data.frame(
+        Type = as.character(expected$Type),
+        y = expected[[estimate_column]],
+        stringsAsFactors = FALSE
+      ),
+      by = "Type"
+    )
+    expect_identical(nrow(compared), 4L)
+    expect_equal(compared$x, compared$y, tolerance = 1e-6)
+  }
+})
+
+
+test_that("estimate_means, backend emmeans, ordinal, threshold modes", {
+  skip_if_not_installed("emmeans")
   data(housing, package = "MASS")
   m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
+
+  modes <- c(
+    cum.prob = "Probability",
+    exc.prob = "Probability",
+    linear.predictor = "Linear_predictor"
+  )
+  for (mode in names(modes)) {
+    out <- estimate_means(m, "Type", predict = mode, backend = "emmeans")
+    expect_identical(nrow(out), 8L)
+    expect_setequal(as.character(out$Threshold), c("Low|Medium", "Medium|High"))
+    expect_true(modes[[mode]] %in% colnames(out))
+    expect_identical(attributes(out)$coef_name, modes[[mode]])
+    expected <- as.data.frame(suppressMessages(
+      emmeans::emmeans(m, c("Type", "cut"), mode = mode)
+    ))
+    estimate_column <- setdiff(
+      colnames(expected),
+      c("Type", "cut", "SE", "df", "asymp.LCL", "asymp.UCL")
+    )
+    compared <- merge(
+      data.frame(
+        Type = as.character(out$Type),
+        Threshold = as.character(out$Threshold),
+        x = out[[modes[[mode]]]],
+        stringsAsFactors = FALSE
+      ),
+      data.frame(
+        Type = as.character(expected$Type),
+        Threshold = as.character(expected$cut),
+        y = expected[[estimate_column]],
+        stringsAsFactors = FALSE
+      ),
+      by = c("Type", "Threshold")
+    )
+    expect_identical(nrow(compared), 8L)
+    expect_equal(compared$x, compared$y, tolerance = 1e-6)
+  }
+})
+
+
+test_that("estimate_means, backend emmeans, ordinal, clm and glmmTMB", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("ordinal")
+  data(housing, package = "MASS")
+  m <- ordinal::clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
   out <- estimate_means(m, "Type", predict = "prob", backend = "emmeans")
-  expect_equal(out$Probability, c(0.33333, 0.33333, 0.33333, 0.33333))
+  expect_identical(nrow(out), 12L)
+  compared <- .compare_ordinal_backends(out, estimate_means(m, "Type"), "Type")
+  expect_identical(nrow(compared), 12L)
+  expect_equal(compared$Probability.x, compared$Probability.y, tolerance = 1e-6)
 
   skip_if_not_installed("glmmTMB", minimum_version = "1.1.15.2")
   m <- glmmTMB::glmmTMB(
@@ -45,7 +203,31 @@ test_that("estimate_means with backend emmeans work for ordinal", {
     family = glmmTMB::ordinal()
   )
   out <- estimate_means(m, "Type", predict = "prob", backend = "emmeans")
-  expect_equal(out$Probability, c(0.33333, 0.33333, 0.33333, 0.33333))
+  expect_identical(nrow(out), 12L)
+  compared <- .compare_ordinal_backends(out, estimate_means(m, "Type"), "Type")
+  expect_identical(nrow(compared), 12L)
+  expect_equal(compared$Probability.x, compared$Probability.y, tolerance = 1e-6)
+})
+
+
+test_that("estimate_means, backend emmeans, ordinal, transformed response", {
+  skip_if_not_installed("emmeans")
+  data(housing, package = "MASS")
+  housing$SatNum <- as.integer(housing$Sat)
+  # `Hess = TRUE`, else `vcov()` re-fits the model in an environment where
+  # `SatNum` does not exist
+  m <- MASS::polr(
+    factor(SatNum) ~ Infl + Type + Cont,
+    weights = Freq,
+    data = housing,
+    Hess = TRUE
+  )
+  out <- estimate_means(m, "Type", backend = "emmeans")
+  expect_identical(nrow(out), 12L)
+  expect_setequal(as.character(out$Response), c("1", "2", "3"))
+  compared <- .compare_ordinal_backends(out, estimate_means(m, "Type"), "Type")
+  expect_identical(nrow(compared), 12L)
+  expect_equal(compared$Probability.x, compared$Probability.y, tolerance = 1e-6)
 })