From b61b648aa472099d76a399db2d96779f25ba6d9a Mon Sep 17 00:00:00 2001
From: Daniel <mail@danielluedecke.de>
Date: Wed, 30 Sep 2026 12:53:02 +0200
Subject: [PATCH 1/4] Support ordinal (argument `mode` for emmeans backend

---
 DESCRIPTION                   |  2 +-
 NEWS.md                       |  3 +++
 R/get_emmeans.R               | 16 +++++++++++++++-
 tests/testthat/test-ordinal.R | 17 +++++++++++++++++
 4 files changed, 36 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 3317cfd32..50974ce38 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.7
+Version: 0.17.0.8
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
diff --git a/NEWS.md b/NEWS.md
index 227120259..77bbbc05c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -18,6 +18,9 @@
   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.
+
 ## Bug fixes
 
 * Fixed issue in `estimate_contrasts()` with wrong assignment of estimates in
diff --git a/R/get_emmeans.R b/R/get_emmeans.R
index 2bf09f01f..3f93d66fb 100644
--- a/R/get_emmeans.R
+++ b/R/get_emmeans.R
@@ -77,7 +77,11 @@ get_emmeans <- function(
     }
   } else {
     dpars <- FALSE
-    fun_args$type <- predict
+    if (predict %in% .emmeans_ordinal_types) {
+      fun_args$mode <- predict
+    } else {
+      fun_args$type <- predict
+    }
   }
 
   # add dots
@@ -122,6 +126,16 @@ get_emmeans <- function(
 # HELPERS (guess arguments) -----------------------------------------------
 # =========================================================================
 
+.emmeans_ordinal_types <- c(
+  "latent",
+  "linear.predictor",
+  "cum.prob",
+  "exc.prob",
+  "prob",
+  "mean.class"
+)
+
+
 #' @keywords internal
 .guess_emmeans_arguments <- function(model, by = NULL, verbose = TRUE, ...) {
   # Gather info
diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R
index cc0d0b321..7a45254a1 100644
--- a/tests/testthat/test-ordinal.R
+++ b/tests/testthat/test-ordinal.R
@@ -32,6 +32,23 @@ test_that("estimate_relation prints ordinal models correctly", {
 })
 
 
+test_that("estimate_means with backend emmeans work for ordinal", {
+  data(housing, package = "MASS")
+  m <- MASS::polr(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))
+
+  skip_if_not_installed("glmmTMB", minimum_version = "1.1.15.2")
+  m <- glmmTMB::glmmTMB(
+    Sat ~ Infl + Type + Cont,
+    data = housing,
+    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))
+})
+
+
 test_that("estimate_means, print bracl", {
   skip_if_not_installed("brglm2")
   # required for the penguins dataset, which was added in R 4.5.0

From ff69ccbad6a72251db15a17abb31f2bfb95e6031 Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Wed, 30 Sep 2026 16:43:08 -0500
Subject: [PATCH 2/4] Return per-category and per-threshold estimates for
 ordinal models with the emmeans backend (#670)

* Add tests for ordinal modes with the emmeans backend (failing)

* Return per-category and per-threshold estimates for ordinal models with the emmeans backend

* Update NEWS entry for ordinal models with the emmeans backend

* Fix ordinal emmeans output for transformed responses, plot names, and docs

- Name the response pseudo-factor after the formula's left-hand side, so
  that models with e.g. factor(y) as response no longer fail.
- Add Latent, Linear_predictor, and Mean_class to the known estimate names.
- Document the ordinal modes and the "prob" default for estimate_means().
- Check the number of matched rows in the glmmTMB test.
---
 NEWS.md                       |  11 +-
 R/estimate_means.R            |   6 ++
 R/get_emmeans.R               |  82 +++++++++++++++
 R/utils.R                     |   5 +-
 man/estimate_contrasts.Rd     |   6 ++
 man/estimate_means.Rd         |   6 ++
 man/estimate_slopes.Rd        |   6 ++
 man/get_emmeans.Rd            |   6 ++
 tests/testthat/test-ordinal.R | 188 +++++++++++++++++++++++++++++++++-
 9 files changed, 310 insertions(+), 6 deletions(-)

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)
 })
 
 

From 41a49b34c06cc69c7787e515196cc9b93961eaeb Mon Sep 17 00:00:00 2001
From: Daniel <mail@danielluedecke.de>
Date: Wed, 30 Sep 2026 23:44:25 +0200
Subject: [PATCH 3/4] Update R/get_emmeans.R

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 R/get_emmeans.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/get_emmeans.R b/R/get_emmeans.R
index 920efc5c7..ce0c5f07c 100644
--- a/R/get_emmeans.R
+++ b/R/get_emmeans.R
@@ -187,7 +187,8 @@ get_emmeans <- function(
 # 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,
+  switch(
+    mode,
     prob = insight::find_terms(model, verbose = FALSE)$response,
     cum.prob = ,
     exc.prob = ,

From 28fd61ec5e4169d7b8c1796aa9c71c33702c1912 Mon Sep 17 00:00:00 2001
From: Daniel <mail@danielluedecke.de>
Date: Thu, 1 Oct 2026 12:59:17 +0200
Subject: [PATCH 4/4] skip test

---
 tests/testthat/test-ordinal.R | 23 +++++++++++++++--------
 1 file changed, 15 insertions(+), 8 deletions(-)

diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R
index 0cebbd25c..af8ea05af 100644
--- a/tests/testthat/test-ordinal.R
+++ b/tests/testthat/test-ordinal.R
@@ -92,8 +92,12 @@ test_that("estimate_means, backend emmeans, ordinal, default predict", {
   expect_identical(
     paste(out$Level1, out$Level2),
     c(
-      "Apartment Tower", "Atrium Apartment", "Atrium Tower",
-      "Terrace Apartment", "Terrace Atrium", "Terrace Tower"
+      "Apartment Tower",
+      "Atrium Apartment",
+      "Atrium Tower",
+      "Terrace Apartment",
+      "Terrace Atrium",
+      "Terrace Tower"
     )
   )
   expect_equal(
@@ -115,9 +119,7 @@ test_that("estimate_means, backend emmeans, ordinal, mean.class and latent", {
     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)
-    ))
+    expected <- as.data.frame(suppressMessages(emmeans::emmeans(m, "Type", mode = mode)))
     estimate_column <- setdiff(
       colnames(expected),
       c("Type", "SE", "df", "asymp.LCL", "asymp.UCL")
@@ -157,9 +159,11 @@ test_that("estimate_means, backend emmeans, ordinal, threshold modes", {
     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)
-    ))
+    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")
@@ -196,7 +200,10 @@ test_that("estimate_means, backend emmeans, ordinal, clm and glmmTMB", {
   expect_identical(nrow(compared), 12L)
   expect_equal(compared$Probability.x, compared$Probability.y, tolerance = 1e-6)
 
+  # ordinal glmmTMB not supported in marginaleffects <= 1.0.0
   skip_if_not_installed("glmmTMB", minimum_version = "1.1.15.2")
+  skip_if_not_installed("marginaleffects", minimum_version = "1.0.1")
+
   m <- glmmTMB::glmmTMB(
     Sat ~ Infl + Type + Cont,
     data = housing,