From b62054688cf2c69c8856a49bcb1fc373e20129a8 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 20 Jul 2026 07:29:00 +0000
Subject: [PATCH 1/9] Initial plan


From 764fc3ed94ea2b579ec271d78b365accadb036b5 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 20 Jul 2026 07:47:24 +0000
Subject: [PATCH 2/9] Split Response levels from Level1/Level2 in
 estimate_contrasts for response-category models

---
 DESCRIPTION                              |  2 +-
 NEWS.md                                  |  6 ++
 R/format.R                               | 87 +++++++++++++++++++++++-
 tests/testthat/test-estimate_contrasts.R | 23 +++++++
 4 files changed, 116 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 80e160487..ab5ce2f98 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.16.0.10
+Version: 0.16.0.11
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
diff --git a/NEWS.md b/NEWS.md
index 09cfb430a..ffd28fbeb 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -8,6 +8,12 @@
 
 ## Changes
 
+* For contrasts of models where predictions vary by response category
+  (e.g., categorical, ordinal/cumulative or multinomial models, as well as
+  multivariate response models), `estimate_contrasts()` now returns separate
+  `Response1` and `Response2` columns instead of merging the response
+  category label into the `Level1` and `Level2` columns (#646).
+
 * The `iterations` argument can now also be used for Bayesian model in
   `estimate_means()`, `estimate_slopes()`, and `estimate_contrasts()`, which
   is then passed to the `ndraws` argument and controls how many samples are
diff --git a/R/format.R b/R/format.R
index df5f987ee..f60f42f6c 100644
--- a/R/format.R
+++ b/R/format.R
@@ -512,9 +512,21 @@ format.marginaleffects_contrasts <- function(
         }
       }
 
+      # for models where predictions are made for each level of the response
+      # (e.g., categorical, ordinal/cumulative or multinomial models, as well
+      # as multivariate response models), marginaleffects merges the response
+      # level and the contrasted factor level into the same string (e.g.,
+      # "setosa 3"). We split those into separate "Response1" and "Response2"
+      # columns, so contrasted factor levels and response levels are not
+      # confounded
+      params <- .split_contrast_response_levels(params, model)
+
       # we need to update these variables, because these are the new column
       # names for contrasts and focal terms
-      contrast <- c("Level1", "Level2")
+      contrast <- intersect(
+        c("Level1", "Response1", "Level2", "Response2"),
+        colnames(params)
+      )
 
       # remove old column
       x$Parameter <- NULL
@@ -570,6 +582,79 @@ format.marginaleffects_contrasts <- function(
 # variable names to levels. We first find the intersection of all
 # levels from all current contrast terms
 
+# for models where predictions are made separately for each level of the
+# response (categorical, ordinal/cumulative, multinomial or multivariate
+# models), marginaleffects merges the response level and the contrasted
+# factor level into the same string, e.g. "setosa 3". This function splits
+# "Level1" and "Level2" into "Response1"/"Level1" and "Response2"/"Level2",
+# so response levels and contrasted factor levels are not confounded. If
+# splitting isn't possible (e.g., not a response-level model, or the response
+# levels cannot be unambiguously detected), `params` is returned unchanged.
+
+.split_contrast_response_levels <- function(params, model) {
+  if (is.null(model) || !all(c("Level1", "Level2") %in% colnames(params))) {
+    return(params)
+  }
+
+  model_info <- .safe(insight::model_info(model, response = 1, verbose = FALSE))
+  is_response_model <- isTRUE(model_info$is_multinomial) ||
+    isTRUE(model_info$is_categorical) ||
+    isTRUE(model_info$is_ordinal) ||
+    isTRUE(model_info$is_cumulative) ||
+    isTRUE(insight::is_multivariate(model))
+
+  if (!is_response_model) {
+    return(params)
+  }
+
+  # response levels, used to detect and split off the response-category
+  # label from the contrasted factor level
+  response_levels <- .safe(unique(as.character(insight::get_response(model, verbose = FALSE))))
+  if (isTRUE(insight::is_multivariate(model))) {
+    response_levels <- unique(c(response_levels, insight::find_response(model)))
+  }
+  if (!length(response_levels)) {
+    return(params)
+  }
+
+  # sort by length, descending, to avoid partial matches for levels that are
+  # substrings of other levels
+  response_levels <- response_levels[order(nchar(response_levels), decreasing = TRUE)]
+
+  split1 <- .split_response_prefix(params$Level1, response_levels)
+  split2 <- .split_response_prefix(params$Level2, response_levels)
+
+  # sanity check - only split if we found a response label for *all* rows, in
+  # both columns, otherwise keep the original (unsplit) columns
+  if (all(!is.na(split1$Response)) && all(!is.na(split2$Response))) {
+    params$Level1 <- factor(split1$Level, levels = unique(split1$Level))
+    params$Response1 <- factor(split1$Response, levels = unique(split1$Response))
+    params$Level2 <- factor(split2$Level, levels = unique(split2$Level))
+    params$Response2 <- factor(split2$Response, levels = unique(split2$Response))
+  }
+
+  params
+}
+
+
+# split off a known "response level" prefix from a vector of strings, e.g.
+# "setosa 3" -> Response = "setosa", Level = "3"
+
+.split_response_prefix <- function(x, response_levels) {
+  x_chr <- as.character(x)
+  out <- data.frame(Response = NA_character_, Level = x_chr, stringsAsFactors = FALSE)
+  for (i in response_levels) {
+    pattern <- paste0("^", i, "[[:space:]]+")
+    matches <- is.na(out$Response) & grepl(pattern, x_chr)
+    if (any(matches)) {
+      out$Response[matches] <- i
+      out$Level[matches] <- insight::trim_ws(sub(pattern, "", x_chr[matches]))
+    }
+  }
+  out
+}
+
+
 .fix_duplicated_contrastlevels <- function(params, contrast_names) {
   multiple_levels <- Reduce(
     function(i, j) intersect(i, j),
diff --git a/tests/testthat/test-estimate_contrasts.R b/tests/testthat/test-estimate_contrasts.R
index 081befc19..0c8e52e5a 100644
--- a/tests/testthat/test-estimate_contrasts.R
+++ b/tests/testthat/test-estimate_contrasts.R
@@ -1815,3 +1815,26 @@ test_that("estimate_contrast, p-adjust tukey works for contrasting slopes", {
   # here to avoid fragile tests while still ensuring close agreement.
   expect_equal(out1$p.value, out2$p, tolerance = 1e-2)
 })
+
+
+test_that("estimate_contrast, categorical/multinomial response models split off Response levels", {
+  skip_if_not_installed("nnet")
+
+  data(iris)
+  m <- nnet::multinom(Species ~ Sepal.Width, data = iris, trace = FALSE)
+  out <- estimate_contrasts(m, contrast = "Sepal.Width = c(2, 3)", backend = "marginaleffects")
+
+  # Level1/Response1/Level2/Response2 columns, in that order
+  expect_named(
+    out,
+    c(
+      "Level1", "Response1", "Level2", "Response2", "Difference", "SE",
+      "CI_low", "CI_high", "t", "df", "p"
+    )
+  )
+  # contrasted levels no longer contain the response category label
+  expect_true(all(as.character(out$Level1) %in% c("2", "3")))
+  expect_true(all(as.character(out$Level2) %in% c("2", "3")))
+  expect_true(all(as.character(out$Response1) %in% levels(iris$Species)))
+  expect_true(all(as.character(out$Response2) %in% levels(iris$Species)))
+})

From 0ad97eff376aeac1389cd550a4665f1ee91e9233 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 20 Jul 2026 07:50:21 +0000
Subject: [PATCH 3/9] Address code review feedback: refine model_info lookup
 and response level ordering

---
 R/format.R                               | 46 +++++++++++++++---------
 tests/testthat/test-estimate_contrasts.R | 16 +++++++++
 2 files changed, 46 insertions(+), 16 deletions(-)

diff --git a/R/format.R b/R/format.R
index f60f42f6c..9482795d9 100644
--- a/R/format.R
+++ b/R/format.R
@@ -519,10 +519,12 @@ format.marginaleffects_contrasts <- function(
       # "setosa 3"). We split those into separate "Response1" and "Response2"
       # columns, so contrasted factor levels and response levels are not
       # confounded
-      params <- .split_contrast_response_levels(params, model)
+      params <- .split_contrast_response_levels(params, model, x)
 
       # we need to update these variables, because these are the new column
-      # names for contrasts and focal terms
+      # names for contrasts and focal terms. We use `intersect()` here because
+      # "Response1"/"Response2" are only present if the contrasted factor and
+      # response category levels were successfully split above
       contrast <- intersect(
         c("Level1", "Response1", "Level2", "Response2"),
         colnames(params)
@@ -591,26 +593,36 @@ format.marginaleffects_contrasts <- function(
 # splitting isn't possible (e.g., not a response-level model, or the response
 # levels cannot be unambiguously detected), `params` is returned unchanged.
 
-.split_contrast_response_levels <- function(params, model) {
+.split_contrast_response_levels <- function(params, model, x = NULL) {
   if (is.null(model) || !all(c("Level1", "Level2") %in% colnames(params))) {
     return(params)
   }
 
-  model_info <- .safe(insight::model_info(model, response = 1, verbose = FALSE))
+  # prefer the model info already stored as attribute (consistent with how
+  # other format-methods retrieve model information), and only compute it
+  # ourselves as a fallback
+  model_info <- attr(x, "model_info", exact = TRUE)
+  if (is.null(model_info)) {
+    model_info <- .safe(insight::model_info(model, verbose = FALSE))
+  }
+  is_multivariate <- isTRUE(.safe(insight::is_multivariate(model)))
   is_response_model <- isTRUE(model_info$is_multinomial) ||
     isTRUE(model_info$is_categorical) ||
     isTRUE(model_info$is_ordinal) ||
     isTRUE(model_info$is_cumulative) ||
-    isTRUE(insight::is_multivariate(model))
+    is_multivariate
 
   if (!is_response_model) {
     return(params)
   }
 
   # response levels, used to detect and split off the response-category
-  # label from the contrasted factor level
-  response_levels <- .safe(unique(as.character(insight::get_response(model, verbose = FALSE))))
-  if (isTRUE(insight::is_multivariate(model))) {
+  # label from the contrasted factor level. We keep the original order of
+  # the response levels (e.g., factor level order), so that "Response1" and
+  # "Response2" use a meaningful, consistent level order
+  response <- .safe(insight::get_response(model, verbose = FALSE))
+  response_levels <- if (is.factor(response)) levels(response) else unique(as.character(response))
+  if (is_multivariate) {
     response_levels <- unique(c(response_levels, insight::find_response(model)))
   }
   if (!length(response_levels)) {
@@ -618,19 +630,21 @@ format.marginaleffects_contrasts <- function(
   }
 
   # sort by length, descending, to avoid partial matches for levels that are
-  # substrings of other levels
-  response_levels <- response_levels[order(nchar(response_levels), decreasing = TRUE)]
+  # substrings of other levels (used for pattern matching only)
+  search_levels <- response_levels[order(nchar(response_levels), decreasing = TRUE)]
 
-  split1 <- .split_response_prefix(params$Level1, response_levels)
-  split2 <- .split_response_prefix(params$Level2, response_levels)
+  split1 <- .split_response_prefix(params$Level1, search_levels)
+  split2 <- .split_response_prefix(params$Level2, search_levels)
 
   # sanity check - only split if we found a response label for *all* rows, in
   # both columns, otherwise keep the original (unsplit) columns
   if (all(!is.na(split1$Response)) && all(!is.na(split2$Response))) {
-    params$Level1 <- factor(split1$Level, levels = unique(split1$Level))
-    params$Response1 <- factor(split1$Response, levels = unique(split1$Response))
-    params$Level2 <- factor(split2$Level, levels = unique(split2$Level))
-    params$Response2 <- factor(split2$Response, levels = unique(split2$Response))
+    level_order <- unique(c(split1$Level, split2$Level))
+    response_order <- intersect(response_levels, unique(c(split1$Response, split2$Response)))
+    params$Level1 <- factor(split1$Level, levels = level_order)
+    params$Response1 <- factor(split1$Response, levels = response_order)
+    params$Level2 <- factor(split2$Level, levels = level_order)
+    params$Response2 <- factor(split2$Response, levels = response_order)
   }
 
   params
diff --git a/tests/testthat/test-estimate_contrasts.R b/tests/testthat/test-estimate_contrasts.R
index 0c8e52e5a..28dc82f2c 100644
--- a/tests/testthat/test-estimate_contrasts.R
+++ b/tests/testthat/test-estimate_contrasts.R
@@ -1837,4 +1837,20 @@ test_that("estimate_contrast, categorical/multinomial response models split off
   expect_true(all(as.character(out$Level2) %in% c("2", "3")))
   expect_true(all(as.character(out$Response1) %in% levels(iris$Species)))
   expect_true(all(as.character(out$Response2) %in% levels(iris$Species)))
+
+  # contrasts within the same response category should match the difference
+  # between the corresponding marginal means for that category
+  means <- estimate_means(m, by = "Sepal.Width = c(2, 3)", backend = "marginaleffects")
+  for (resp in levels(iris$Species)) {
+    row <- out[
+      as.character(out$Level1) == "3" &
+        as.character(out$Level2) == "2" &
+        as.character(out$Response1) == resp &
+        as.character(out$Response2) == resp,
+    ]
+    expect_identical(nrow(row), 1L)
+    mean_3 <- means$Probability[means$Sepal.Width == "3" & means$Response == resp]
+    mean_2 <- means$Probability[means$Sepal.Width == "2" & means$Response == resp]
+    expect_equal(row$Difference, mean_3 - mean_2, tolerance = 1e-4)
+  }
 })

From b9ffc2b541764db5e78e71424dcada16ce76c600 Mon Sep 17 00:00:00 2001
From: Dominique Makowski <dom.mak19@gmail.com>
Date: Mon, 20 Jul 2026 09:10:37 +0100
Subject: [PATCH 4/9] Apply suggestions from code review

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

diff --git a/R/format.R b/R/format.R
index 9482795d9..79ba4d0a2 100644
--- a/R/format.R
+++ b/R/format.R
@@ -621,7 +621,11 @@ format.marginaleffects_contrasts <- function(
   # the response levels (e.g., factor level order), so that "Response1" and
   # "Response2" use a meaningful, consistent level order
   response <- .safe(insight::get_response(model, verbose = FALSE))
-  response_levels <- if (is.factor(response)) levels(response) else unique(as.character(response))
+  response_levels <- if (is.factor(response)) {
+    levels(response)
+  } else {
+    unique(as.character(response))
+  }
   if (is_multivariate) {
     response_levels <- unique(c(response_levels, insight::find_response(model)))
   }

From 158f2688e28102dad778170f07b620e27ef14ab7 Mon Sep 17 00:00:00 2001
From: Dominique Makowski <dom.mak19@gmail.com>
Date: Mon, 20 Jul 2026 09:24:52 +0100
Subject: [PATCH 5/9] Potential fix for pull request finding

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
---
 R/format.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/format.R b/R/format.R
index 79ba4d0a2..5d109ff35 100644
--- a/R/format.R
+++ b/R/format.R
@@ -662,7 +662,8 @@ format.marginaleffects_contrasts <- function(
   x_chr <- as.character(x)
   out <- data.frame(Response = NA_character_, Level = x_chr, stringsAsFactors = FALSE)
   for (i in response_levels) {
-    pattern <- paste0("^", i, "[[:space:]]+")
+    i_esc <- gsub("([][{}()+*^$|\\\\?.])", "\\\\\\1", i, perl = TRUE)
+    pattern <- paste0("^", i_esc, "[[:space:]]+")
     matches <- is.na(out$Response) & grepl(pattern, x_chr)
     if (any(matches)) {
       out$Response[matches] <- i

From f0fbcda157c63922b758a623c6b2f4d6f5055402 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 20 Jul 2026 08:27:18 +0000
Subject: [PATCH 6/9] Only coerce non-factor contrast columns to factor

---
 R/format.R | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/R/format.R b/R/format.R
index 5d109ff35..eab247495 100644
--- a/R/format.R
+++ b/R/format.R
@@ -538,7 +538,9 @@ format.marginaleffects_contrasts <- function(
 
       # make sure terms are factors, for data_arrange later
       for (i in contrast) {
-        x[[i]] <- factor(x[[i]], levels = unique(x[[i]]))
+        if (!is.factor(x[[i]])) {
+          x[[i]] <- factor(x[[i]], levels = unique(x[[i]]))
+        }
       }
 
       # filter contrast-predictor levels, if requested (e.g., `contrast = "x=c('a', 'b')"`)

From aa43c336c927bf9689aafae1e9c90e19181f6610 Mon Sep 17 00:00:00 2001
From: Daniel <mail@danielluedecke.de>
Date: Mon, 20 Jul 2026 10:33:14 +0200
Subject: [PATCH 7/9] air formatting

---
 R/format.R | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/R/format.R b/R/format.R
index eab247495..3680324e3 100644
--- a/R/format.R
+++ b/R/format.R
@@ -646,7 +646,10 @@ format.marginaleffects_contrasts <- function(
   # both columns, otherwise keep the original (unsplit) columns
   if (all(!is.na(split1$Response)) && all(!is.na(split2$Response))) {
     level_order <- unique(c(split1$Level, split2$Level))
-    response_order <- intersect(response_levels, unique(c(split1$Response, split2$Response)))
+    response_order <- intersect(
+      response_levels,
+      unique(c(split1$Response, split2$Response))
+    )
     params$Level1 <- factor(split1$Level, levels = level_order)
     params$Response1 <- factor(split1$Response, levels = response_order)
     params$Level2 <- factor(split2$Level, levels = level_order)

From 8f23382d25bb5a7ddd02883aea7e02272873566e Mon Sep 17 00:00:00 2001
From: Daniel <mail@danielluedecke.de>
Date: Mon, 20 Jul 2026 10:36:07 +0200
Subject: [PATCH 8/9] code style

---
 R/format.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/format.R b/R/format.R
index 3680324e3..d27b5cf83 100644
--- a/R/format.R
+++ b/R/format.R
@@ -623,10 +623,10 @@ format.marginaleffects_contrasts <- function(
   # the response levels (e.g., factor level order), so that "Response1" and
   # "Response2" use a meaningful, consistent level order
   response <- .safe(insight::get_response(model, verbose = FALSE))
-  response_levels <- if (is.factor(response)) {
-    levels(response)
+  if (is.factor(response)) {
+    response_levels <- levels(response)
   } else {
-    unique(as.character(response))
+    response_levels <- unique(as.character(response))
   }
   if (is_multivariate) {
     response_levels <- unique(c(response_levels, insight::find_response(model)))

From 42b58bbd9ed3e69487ddb53b0013452e7fd27d2a Mon Sep 17 00:00:00 2001
From: Daniel <mail@danielluedecke.de>
Date: Mon, 20 Jul 2026 19:12:51 +0200
Subject: [PATCH 9/9] validate against marginaleffects

---
 tests/testthat/test-estimate_contrasts.R | 42 ++++++++++++++++++++++--
 1 file changed, 39 insertions(+), 3 deletions(-)

diff --git a/tests/testthat/test-estimate_contrasts.R b/tests/testthat/test-estimate_contrasts.R
index 28dc82f2c..8f01f1e15 100644
--- a/tests/testthat/test-estimate_contrasts.R
+++ b/tests/testthat/test-estimate_contrasts.R
@@ -1822,14 +1822,27 @@ test_that("estimate_contrast, categorical/multinomial response models split off
 
   data(iris)
   m <- nnet::multinom(Species ~ Sepal.Width, data = iris, trace = FALSE)
-  out <- estimate_contrasts(m, contrast = "Sepal.Width = c(2, 3)", backend = "marginaleffects")
+  out <- estimate_contrasts(
+    m,
+    contrast = "Sepal.Width = c(2, 3)",
+    backend = "marginaleffects"
+  )
 
   # Level1/Response1/Level2/Response2 columns, in that order
   expect_named(
     out,
     c(
-      "Level1", "Response1", "Level2", "Response2", "Difference", "SE",
-      "CI_low", "CI_high", "t", "df", "p"
+      "Level1",
+      "Response1",
+      "Level2",
+      "Response2",
+      "Difference",
+      "SE",
+      "CI_low",
+      "CI_high",
+      "t",
+      "df",
+      "p"
     )
   )
   # contrasted levels no longer contain the response category label
@@ -1838,6 +1851,29 @@ test_that("estimate_contrast, categorical/multinomial response models split off
   expect_true(all(as.character(out$Response1) %in% levels(iris$Species)))
   expect_true(all(as.character(out$Response2) %in% levels(iris$Species)))
 
+  out2 <- marginaleffects::avg_predictions(
+    m,
+    by = "Sepal.Width",
+    newdata = data.frame(Sepal.Width = c(2, 3)),
+    hypothesis = ~pairwise
+  )
+
+  expect_identical(
+    paste0(
+      "(",
+      out$Response1,
+      " ",
+      out$Level1,
+      ") - (",
+      out$Response2,
+      " ",
+      out$Level2,
+      ")"
+    ),
+    out2$hypothesis
+  )
+  expect_equal(out$Difference, out2$estimate, tolerance = 1e-4)
+
   # contrasts within the same response category should match the difference
   # between the corresponding marginal means for that category
   means <- estimate_means(m, by = "Sepal.Width = c(2, 3)", backend = "marginaleffects")