Color plots of emmeans ordinal output by response category or threshold

NEWS.md

@@ -28,6 +28,15 @@
   now `predict = "prob"`, which returns the same probabilities as the
   `"marginaleffects"` backend.
 
+* For ordinal models, `plot()` of `estimate_means()` output with
+  `backend = "emmeans"` now maps the `Response` or `Threshold` column to the
+  color with one `by` variable and to the facets with more. The
+  `"marginaleffects"` backend already does this for `Response`.
+
+* With `backend = "emmeans"`, `plot()` now shows predicted probabilities on a
+  percent y-axis, `"cum.prob"` and `"exc.prob"` included. This applies to
+  binomial, beta, ordinal, and multinomial models.
+
 * `comparison = "ratio"` is a new shortcut to compute pairwise comparisons
   as ratios.
 

R/estimate_means.R

@@ -538,6 +538,12 @@ estimate_means <- function(
   # add attributes from workhorse function
   attributes(means) <- utils::modifyList(attributes(means), info[.info_elements()])
 
+  # output from the emmeans backend has no `model_info` attribute, which we
+  # need for plotting, e.g. to color the response categories of ordinal models
+  if (backend == "emmeans") {
+    attr(means, "model_info") <- insight::model_info(model, response = 1, verbose = FALSE)
+  }
+
   # Output
   class(means) <- unique(c("estimate_means", class(means)))
   means

R/visualisation_recipe_internal.R

@@ -37,6 +37,17 @@
     data$Response <- factor(data$Response, levels = unique(data$Response))
   }
 
+  # ordinal models, emmeans backend: cumulative or exceedance probabilities and
+  # the linear predictor are returned for each threshold. We skip this when
+  # a focal predictor is named "Threshold"
+  if (has_response_levels && "Threshold" %in% colnames(data) && !"Threshold" %in% by) {
+    by <- c(by, "Threshold")
+    # emmeans returns a factor with the thresholds in model order
+    if (!is.factor(data$Threshold)) {
+      data$Threshold <- factor(data$Threshold, levels = unique(data$Threshold))
+    }
+  }
+
   # mixture models? if so, we need one more stratification in "by"
   if (isTRUE(model_info$is_mixture) && "Class" %in% colnames(data)) {
     by <- c(by, "Class")
@@ -502,7 +513,10 @@
   if (
     !is.null(response_scale) &&
       response_scale %in%
-        c("response", "expectation", "invlink(link)", "prob", "probs") &&
+        c(
+          "response", "expectation", "invlink(link)", "prob", "probs",
+          "cum.prob", "exc.prob"
+        ) &&
       isTRUE(
         model_info$is_logit |
           model_info$is_binomial |

tests/testthat/test-attributes_estimatefun.R

@@ -12,7 +12,7 @@ test_that("attributes_means, means", {
     c(
       "names", "row.names", "class", "at", "by", "table_title", "table_footer",
       "model", "response", "ci", "backend", "coef_name", "call", "focal_terms",
-      "predict", "transform", "keep_iterations"
+      "predict", "transform", "keep_iterations", "model_info"
     )
   )
   estim <- suppressMessages(estimate_means(model, "Species", backend = "marginaleffects"))

tests/testthat/test-keep_iterations.R

@@ -167,7 +167,7 @@ test_that("estimate_means() - posterior draws, emmeans", {
     c(
       "names", "class", "row.names", "table_title", "table_footer",
       "model", "response", "ci", "backend", "coef_name", "call", "at", "by",
-      "focal_terms", "predict", "transform", "keep_iterations"
+      "focal_terms", "predict", "transform", "keep_iterations", "model_info"
     )
   )
   # fmt: skip
@@ -187,7 +187,7 @@ test_that("estimate_means() - posterior draws, emmeans", {
     c(
       "names", "class", "row.names", "table_title", "table_footer",
       "model", "response", "ci", "backend", "coef_name", "call", "at", "by",
-      "focal_terms", "predict", "transform", "keep_iterations"
+      "focal_terms", "predict", "transform", "keep_iterations", "model_info"
     )
   )
   expect_identical(dim(out), c(10L, 4005L))

tests/testthat/test-visualisation_recipe.R

@@ -208,6 +208,232 @@ test_that("visualization_recipe", {
   # plot(modelbased:::.visualization_recipe(x))
 })
 
+
+# emmeans backend: ordinal and other models ----------------------------------
+
+.recipe_mapping <- function(x, ...) {
+  r <- visualisation_recipe(x, ...)
+  geoms <- vapply(r, function(l) l$geom, character(1))
+  main <- r[[which(geoms %in% c("pointrange", "line"))[1]]]
+  facet <- r[geoms == "facet_wrap"]
+  list(
+    geoms = unname(geoms),
+    x = main$aes$x,
+    color = main$aes$color,
+    facet = if (length(facet)) deparse(facet[[1]]$facets) else NULL,
+    data = attributes(r)$data
+  )
+}
+
+# printed table as a list of character columns
+.printed_table <- function(x) {
+  lapply(as.list(format(x)), as.character)
+}
+
+# the printed table must not depend on the model information
+.expect_same_print <- function(x) {
+  y <- x
+  attr(y, "model_info") <- NULL
+  expect_identical(.printed_table(x), .printed_table(y))
+}
+
+
+test_that("visualization_recipe, emmeans ordinal, color by Response", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("MASS")
+  skip_if_not_installed("ordinal")
+  skip_if_not_installed("marginaleffects", minimum_version = "0.29.0")
+
+  data(housing, package = "MASS")
+  models <- list(
+    polr = MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE),
+    clm = ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing)
+  )
+  for (m in models) {
+    for (backend in c("emmeans", "marginaleffects")) {
+      x <- estimate_means(m, by = "Infl", backend = backend)
+      out <- .recipe_mapping(x)
+      expect_identical(out$x, "Infl")
+      expect_identical(out$color, "Response")
+      expect_null(out$facet)
+      expect_true("scale_y_continuous" %in% out$geoms)
+      expect_identical(attributes(x)$focal_terms, "Infl")
+      if (backend == "emmeans") .expect_same_print(x)
+
+      x <- estimate_means(m, by = c("Infl", "Type"), backend = backend)
+      out <- .recipe_mapping(x)
+      expect_identical(out$x, "Infl")
+      expect_identical(out$color, "Type")
+      expect_identical(out$facet, "~Response")
+      expect_identical(attributes(x)$focal_terms, c("Infl", "Type"))
+      if (backend == "emmeans") .expect_same_print(x)
+    }
+  }
+
+  # printed rows keep their order: the response varies within predictor levels
+  x <- estimate_means(models$polr, by = c("Infl", "Type"), backend = "emmeans")
+  printed <- .printed_table(x)
+  expect_identical(
+    printed$Infl[1:6],
+    c("Low", "Low", "Low", "Medium", "Medium", "Medium")
+  )
+  expect_identical(printed$Type[1:6], rep("Tower", 6))
+  expect_identical(
+    printed$Response[1:6],
+    c("Low", "Medium", "High", "Low", "Medium", "High")
+  )
+
+  # mixed models from package ordinal
+  data(wine, package = "ordinal")
+  m <- ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine)
+  x <- estimate_means(m, by = "temp", backend = "emmeans")
+  out <- .recipe_mapping(x)
+  expect_identical(out$x, "temp")
+  expect_identical(out$color, "Response")
+  .expect_same_print(x)
+})
+
+
+test_that("visualization_recipe, emmeans ordinal, color by Threshold", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("MASS")
+  skip_if_not_installed("ordinal")
+
+  data(housing, package = "MASS")
+  models <- list(
+    polr = MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE),
+    clm = ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing)
+  )
+  for (m in models) {
+    for (p in c("cum.prob", "exc.prob", "linear.predictor")) {
+      x <- estimate_means(m, by = "Infl", backend = "emmeans", predict = p)
+      out <- .recipe_mapping(x)
+      expect_identical(out$x, "Infl")
+      expect_identical(out$color, "Threshold")
+      expect_null(out$facet)
+      expect_identical(
+        "scale_y_continuous" %in% out$geoms,
+        p %in% c("cum.prob", "exc.prob")
+      )
+      expect_identical(attributes(x)$focal_terms, "Infl")
+      .expect_same_print(x)
+
+      x <- estimate_means(m, by = c("Infl", "Type"), backend = "emmeans", predict = p)
+      out <- .recipe_mapping(x)
+      expect_identical(out$x, "Infl")
+      expect_identical(out$color, "Type")
+      expect_identical(out$facet, "~Threshold")
+      .expect_same_print(x)
+    }
+  }
+
+  # thresholds in model order, not in alphabetical order
+  data(wine, package = "ordinal")
+  rating_levels <- c("very low", "low", "medium", "high", "very high")
+  wine$rating <- factor(
+    rating_levels[as.integer(wine$rating)],
+    levels = rating_levels,
+    ordered = TRUE
+  )
+  m <- MASS::polr(rating ~ temp + contact, data = wine, Hess = TRUE)
+  x <- estimate_means(m, by = "temp", backend = "emmeans", predict = "cum.prob")
+  out <- .recipe_mapping(x)
+  expect_identical(out$color, "Threshold")
+  expect_identical(
+    levels(out$data$Threshold),
+    c("very low|low", "low|medium", "medium|high", "high|very high")
+  )
+  expect_identical(
+    as.character(unique(out$data$Threshold)),
+    c("very low|low", "low|medium", "medium|high", "high|very high")
+  )
+  .expect_same_print(x)
+})
+
+
+test_that("visualization_recipe, ordinal, focal predictor named Threshold", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("MASS")
+  skip_if_not_installed("marginaleffects", minimum_version = "0.29.0")
+
+  data(housing, package = "MASS")
+  d <- housing
+  d$Threshold <- d$Infl
+  m <- MASS::polr(Sat ~ Threshold + Type, weights = Freq, data = d, Hess = TRUE)
+  for (backend in c("emmeans", "marginaleffects")) {
+    x <- estimate_means(m, by = "Threshold", backend = backend)
+    out <- .recipe_mapping(x)
+    expect_identical(out$x, "Threshold")
+    expect_identical(out$color, "Response")
+    expect_null(out$facet)
+  }
+})
+
+
+test_that("visualization_recipe, emmeans, percent scale for logistic models", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("marginaleffects", minimum_version = "0.29.0")
+
+  d <- mtcars
+  d$cyl <- factor(d$cyl)
+  m <- glm(am ~ cyl, family = binomial, data = d)
+  for (backend in c("emmeans", "marginaleffects")) {
+    x <- estimate_means(m, by = "cyl", backend = backend)
+    out <- .recipe_mapping(x)
+    expect_identical(out$geoms, c("pointrange", "labs", "scale_y_continuous"))
+  }
+  .expect_same_print(estimate_means(m, by = "cyl", backend = "emmeans"))
+})
+
+
+test_that("visualization_recipe, emmeans, recipes that do not change", {
+  skip_if_not_installed("emmeans")
+  skip_if_not_installed("MASS")
+
+  m <- lm(Sepal.Length ~ Species, data = iris)
+  x <- estimate_means(m, by = "Species", backend = "emmeans")
+  out <- .recipe_mapping(x, show_data = TRUE)
+  expect_identical(out$geoms, c("jitter", "pointrange", "labs"))
+  expect_identical(out$x, "Species")
+  expect_null(out$color)
+  .expect_same_print(x)
+
+  m <- lm(log(Sepal.Length) ~ Species, data = iris)
+  x <- estimate_means(m, by = "Species", backend = "emmeans")
+  out <- .recipe_mapping(x, show_data = TRUE)
+  expect_identical(out$geoms, c("jitter", "pointrange", "labs"))
+  expect_identical(out$x, "Species")
+  expect_null(out$color)
+  .expect_same_print(x)
+
+  d <- mtcars
+  d$cyl <- factor(d$cyl)
+  m <- glm(carb ~ cyl, family = poisson, data = d)
+  x <- estimate_means(m, by = "cyl", backend = "emmeans")
+  out <- .recipe_mapping(x)
+  expect_identical(out$geoms, c("pointrange", "labs"))
+  expect_identical(out$x, "cyl")
+  expect_null(out$color)
+  .expect_same_print(x)
+
+  m <- glm(am ~ cyl, family = binomial, data = d)
+  x <- estimate_means(m, by = "cyl", backend = "emmeans", predict = "link")
+  out <- .recipe_mapping(x)
+  expect_identical(out$geoms, c("pointrange", "labs"))
+  expect_identical(out$x, "cyl")
+  expect_null(out$color)
+  .expect_same_print(x)
+
+  data(housing, package = "MASS")
+  m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE)
+  for (p in c("mean.class", "latent")) {
+    x <- estimate_means(m, by = "Infl", backend = "emmeans", predict = p)
+    out <- .recipe_mapping(x)
+    expect_identical(out$geoms, c("pointrange", "labs"))
+    expect_identical(out$x, "Infl")
+    expect_null(out$color)
+    .expect_same_print(x)
+  }
 test_that("visualization_recipe, informative error when nothing goes on x-axis", {
   # intercept-only model
   model <- lm(qsec ~ 1, data = mtcars)