diff --git a/NEWS.md b/NEWS.md
index b0f5ba58e..f99fb75ad 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -116,6 +116,11 @@
   x-axis, e.g. for predictions from intercept-only models, or for
   `estimate_prediction()` without the `by` argument.
 
+* Fixed an "invalid regular expression" error in `estimate_contrasts()` for
+  multivariate models (for example from *brms*) with many observations.
+  `estimate_contrasts()` now uses the names of the response variables, not
+  their observed values, to find the response in the contrast labels.
+
 * The `grid` argument of `visualisation_recipe()` and `plot()` now changes the
   facet grid of plots with four or more `by` variables. Before, these functions
   applied the `facet` argument to the facet grid instead. If `facet` was not
diff --git a/R/format.R b/R/format.R
index 17a253616..6d130db5c 100644
--- a/R/format.R
+++ b/R/format.R
@@ -624,7 +624,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))
-  if (is.factor(response)) {
+  # for multivariate models, the response is a data frame or matrix, and the
+  # response levels are its column names, not the observed values
+  if (is.matrix(response) || is.data.frame(response)) {
+    response_levels <- colnames(response)
+  } else if (is.factor(response)) {
     response_levels <- levels(response)
   } else {
     response_levels <- unique(as.character(response))
diff --git a/tests/testthat/test-brms.R b/tests/testthat/test-brms.R
index 8fe9dca99..681bb014e 100644
--- a/tests/testthat/test-brms.R
+++ b/tests/testthat/test-brms.R
@@ -87,6 +87,36 @@ test_that("estimate_means - brms, multivariate", {
 })
 
 
+test_that("estimate_contrasts - brms, multivariate, many observations", {
+  # the response values must not be used as response levels, else the
+  # regular expressions built from them grow too long and fail
+  set.seed(123)
+  d <- data.frame(
+    y1 = rnorm(600),
+    y2 = rnorm(600),
+    g = factor(sample(c("a", "b", "c"), 600, replace = TRUE))
+  )
+  m <- suppressWarnings(brms::brm(
+    brms::mvbind(y1, y2) ~ g,
+    data = d,
+    refresh = 0,
+    iter = 500,
+    chains = 1,
+    seed = 123
+  ))
+
+  out <- estimate_contrasts(m, "g", resp = "y1")
+  expect_identical(nrow(out), 3L)
+  expect_identical(as.character(out$Level1), c("b", "c", "c"))
+  expect_identical(as.character(out$Level2), c("a", "a", "b"))
+
+  out <- estimate_contrasts(m, "g")
+  expect_identical(nrow(out), 15L)
+  expect_true(all(c("Response1", "Response2") %in% colnames(out)))
+  expect_identical(levels(out$Response1), c("y1", "y2"))
+})
+
+
 test_that("estimate_means - brms, Wiener", {
   skip_if_not_installed("curl")
   skip_if_offline()