diff --git a/DESCRIPTION b/DESCRIPTION
index 3069ae7a8..b1c5da55f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.18.0.1
+Version: 0.18.0.2
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
diff --git a/NEWS.md b/NEWS.md
index 704d59924..cdd4a67c5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,17 @@
 
 ## Bug fixes
 
+* `estimate_contrasts()` now applies the filter values in `by` to contrasts of
+  slopes with `estimate = "average"` and with `backend = "emmeans"`, for
+  example `estimate_contrasts(m, "wt", by = c("cyl = '4'", "am"))`, where
+  `cyl` and `am` are factors. Before, these contrasts compared the slopes of
+  all levels and ignored the filter. With `backend = "emmeans"`, the filter
+  applies only to `by` variables that are all factors or character vectors.
+  If the filter leaves fewer than two slopes to compare,
+  `estimate_contrasts()` now stops with a clear error. Ordinal and
+  multinomial models, and a custom hypothesis in `comparison` (such as
+  `"b1 = b2"`, a function or a matrix), get no such check.
+
 * `p_adjust = "tukey"` and `p_adjust = "sidak"` now give the same p-values as
   *emmeans* in `estimate_contrasts()` and `estimate_slopes()` with
   `backend = "marginaleffects"`:
diff --git a/R/get_emcontrasts.R b/R/get_emcontrasts.R
index a0e301db2..4dc665a17 100644
--- a/R/get_emcontrasts.R
+++ b/R/get_emcontrasts.R
@@ -90,8 +90,34 @@ get_emcontrasts <- function(
         "Please specify the `by` argument to calculate contrasts of slopes."
       )
     }
-    # Run emmeans
+    # filter values in `by`, but not for the response categories of ordinal
+    # models, which are filtered later
+    slope_filter <- .parse_by_filter(by)
+    slope_filter <- slope_filter[setdiff(
+      intersect(names(slope_filter), my_args$by),
+      c(my_args$contrast, my_args$response_filter$response)
+    )]
+    # emmeans holds numeric variables in `by` at their mean, and `at` does not
+    # work for factors that are created in the model formula. So we only
+    # filter when all variables in `by` are factors or character vectors
+    numeric_by <- vapply(my_args$by, function(i) is.numeric(model_data[[i]]), logical(1))
+    if (any(numeric_by)) {
+      slope_filter <- list()
+    }
+    slope_filter <- .present_filter_values(model_data, slope_filter)
+    # sanity check - we need at least two slopes to compare
+    if (
+      (is.character(comparison) && !.is_custom_comparison(comparison)) ||
+        inherits(comparison, "formula")
+    ) {
+      .check_slope_contrast_filter(model, my_args$by, slope_filter)
+    }
+    # Run emmeans. Filters in `by` go into `at`, so only those slopes are
+    # contrasted
     fun_args <- c(fun_args, list(specs = my_args$by, var = my_args$contrast))
+    if (length(slope_filter)) {
+      fun_args$at <- slope_filter
+    }
     if (is.null(ordinal$mode)) {
       estimated <- suppressMessages(do.call(emmeans::emtrends, fun_args))
       emm_by <- NULL
diff --git a/R/get_marginalcontrasts.R b/R/get_marginalcontrasts.R
index 70ae023df..ffb0a3770 100644
--- a/R/get_marginalcontrasts.R
+++ b/R/get_marginalcontrasts.R
@@ -128,22 +128,74 @@ get_marginalcontrasts <- function(
         "Please specify the `by` argument to calculate contrasts of slopes."
       )
     }
+    # filter values in `by`. For `estimate = "average"`, they are already
+    # removed from `by` and saved in `by_filter`
+    if (identical(estimate, "average")) {
+      slope_filter <- my_args$by_filter
+    } else {
+      slope_filter <- .parse_by_filter(my_args$by)
+    }
+    slope_filter <- .present_filter_values(insight::get_data(model, verbose = FALSE), slope_filter)
+    # sanity check - we need at least two slopes to compare. We use the
+    # original `by`, because a group in a `comparison` formula replaces it
+    if (
+      !.is_custom_comparison(comparison) &&
+        !is.function(comparison) &&
+        !is.matrix(comparison)
+    ) {
+      .check_slope_contrast_filter(model, by, slope_filter, estimate)
+    }
     # call slopes with hypothesis argument
-    out <- estimate_slopes(
-      model = model,
-      trend = my_args$contrast,
-      by = my_args$by,
-      ci = ci,
-      estimate = estimate,
-      hypothesis = my_args$comparison_slopes,
-      backend = "marginaleffects",
-      transform = transform,
-      post_process = post_process,
-      iterations = iterations,
-      keep_iterations = keep_iterations,
-      verbose = verbose,
-      ...
-    )
+    .contrast_slopes <- function(slope_by = my_args$by, ...) {
+      estimate_slopes(
+        model = model,
+        trend = my_args$contrast,
+        by = slope_by,
+        ci = ci,
+        estimate = estimate,
+        hypothesis = my_args$comparison_slopes,
+        backend = "marginaleffects",
+        transform = transform,
+        post_process = post_process,
+        iterations = iterations,
+        keep_iterations = keep_iterations,
+        verbose = verbose,
+        ...
+      )
+    }
+    # for `estimate = "average"`, there is no data grid that we can filter, and
+    # the hypothesis compares all slopes, so we cannot filter the contrasts
+    # later. Instead, we only keep the data rows that pass the filter. The
+    # slopes are averaged within each group of `by`, so they don't change.
+    dots <- list(...)
+    if (
+      identical(estimate, "average") &&
+        length(slope_filter) &&
+        is.null(dots$newdata) &&
+        is.null(dots$data)
+    ) {
+      filtered_data <- .filter_data_by(
+        insight::get_data(model, verbose = FALSE),
+        slope_filter
+      )
+      # we add the filter back to `by`, so the data grid that is used to
+      # format the contrast labels only holds the filtered values
+      slope_by <- vapply(
+        my_args$by,
+        function(i) {
+          if (i %in% names(slope_filter)) {
+            paste0(i, "=", insight::safe_deparse(slope_filter[[i]]))
+          } else {
+            i
+          }
+        },
+        character(1),
+        USE.NAMES = FALSE
+      )
+      out <- .contrast_slopes(slope_by = slope_by, newdata = filtered_data, ...)
+    } else {
+      out <- .contrast_slopes(...)
+    }
   } else {
     # for contrasts of categorical predictors, we call avg_predictions
     # -------------------------------------------------------------------------
@@ -260,6 +312,106 @@ get_marginalcontrasts <- function(
 }
 
 
+# filter values in `by` for contrasts of slopes -------------------------------
+
+# returns a named list with the filter values for each variable in `by` that
+# has a filter, like `by = c("Species='setosa'", "x")`
+.parse_by_filter <- function(by) {
+  by_filter <- list()
+  for (i in by[grepl("=", by, fixed = TRUE)]) {
+    filter_value <- insight::trim_ws(unlist(strsplit(i, "=", fixed = TRUE), use.names = FALSE))
+    if (length(filter_value) > 1) {
+      value <- .safe(eval(str2lang(filter_value[2])))
+      if (!is.null(value)) {
+        by_filter[[filter_value[1]]] <- value
+      }
+    }
+  }
+  by_filter
+}
+
+
+# for factors and character vectors, keep only the filter values that are
+# levels in the data
+.present_filter_values <- function(data, by_filter) {
+  for (i in names(by_filter)) {
+    if (!is.null(data[[i]]) && !is.numeric(data[[i]])) {
+      values <- as.character(by_filter[[i]])
+      by_filter[[i]] <- values[values %in% unique(as.character(data[[i]]))]
+    }
+  }
+  by_filter
+}
+
+
+# keep only the data rows whose values pass the filters in `by`
+.filter_data_by <- function(data, by_filter) {
+  keep <- rep(TRUE, nrow(data))
+  for (i in names(by_filter)) {
+    if (!is.null(data[[i]])) {
+      keep <- keep & as.character(data[[i]]) %in% as.character(by_filter[[i]])
+    }
+  }
+  data[keep, , drop = FALSE]
+}
+
+
+# stop when the filters in `by` leave fewer than two slopes to compare
+.check_slope_contrast_filter <- function(model, by, by_filter, estimate = NULL) {
+  if (!length(by_filter) || !length(by)) {
+    return(invisible(NULL))
+  }
+  # models with response categories have one slope per category, so we
+  # cannot count the slopes from `by` alone
+  model_info <- insight::model_info(model, response = 1, verbose = FALSE)
+  if (isTRUE(model_info$is_ordinal) || isTRUE(model_info$is_multinomial) || isTRUE(model_info$is_categorical)) {
+    return(invisible(NULL))
+  }
+  model_data <- insight::get_data(model, verbose = FALSE)
+  by <- insight::trim_ws(gsub("=.*", "", by))
+  # we cannot count the slopes for variables that are not in the data
+  if (!all(by %in% colnames(model_data))) {
+    return(invisible(NULL))
+  }
+  # for `estimate = "average"`, slopes are averaged over the data, so we count
+  # the combinations of the `by` variables in the filtered data
+  if (identical(estimate, "average")) {
+    filtered_data <- .filter_data_by(model_data, by_filter)
+    n_slopes <- nrow(unique(filtered_data[by]))
+  } else {
+    n_slopes <- .count_grid_slopes(model_data, by, by_filter)
+  }
+  if (n_slopes < 2) {
+    insight::format_error(
+      "The filter in `by` leaves fewer than two slopes to compare. Contrasts of slopes need at least two slopes. Remove the filter, or add more values to it."
+    )
+  }
+  invisible(NULL)
+}
+
+
+# number of slopes in the data grid, which has all combinations of the levels
+# or values of the `by` variables
+.count_grid_slopes <- function(model_data, by, by_filter) {
+  n_slopes <- 1
+  for (i in by) {
+    x <- model_data[[i]]
+    values <- by_filter[[i]]
+    if (is.null(values)) {
+      # numeric variables without filter have more than one value in the grid
+      n <- if (is.numeric(x)) 2 else insight::n_unique(x)
+    } else if (is.numeric(x)) {
+      # the data grid uses all filter values of numeric variables
+      n <- insight::n_unique(values)
+    } else {
+      n <- sum(unique(as.character(values)) %in% unique(as.character(x)))
+    }
+    n_slopes <- n_slopes * n
+  }
+  n_slopes
+}
+
+
 # make "comparison" argument compatible -----------------------------------
 # -------------------------------------------------------------------------
 
diff --git a/tests/testthat/test-estimate_contrasts_slopes_filter.R b/tests/testthat/test-estimate_contrasts_slopes_filter.R
new file mode 100644
index 000000000..d61c04998
--- /dev/null
+++ b/tests/testthat/test-estimate_contrasts_slopes_filter.R
@@ -0,0 +1,200 @@
+skip_if_not_installed("emmeans")
+skip_if_not_installed("marginaleffects", minimum_version = "0.29.0")
+
+mt <- mtcars
+mt$g <- factor(mt$cyl)
+mt$a <- factor(mt$am)
+m2 <- lm(mpg ~ wt * g, data = mt)
+m3 <- lm(mpg ~ wt * g * a, data = mt)
+m6 <- lm(mpg ~ wt * g * a + hp, data = mt)
+
+# filters for `by`, the matching `at` for emmeans, and the levels that each
+# filter allows
+filters <- list(
+  F1 = list(
+    by = c("g='4'", "a"),
+    at = list(g = "4"),
+    allowed = list(g = "4", a = c("0", "1"))
+  ),
+  F2 = list(
+    by = c("g=c('4','6')", "a='1'"),
+    at = list(g = c("4", "6"), a = "1"),
+    allowed = list(g = c("4", "6"), a = "1")
+  ),
+  F3 = list(
+    by = c("g", "a='1'"),
+    at = list(a = "1"),
+    allowed = list(g = c("4", "6", "8"), a = "1")
+  ),
+  F4 = list(
+    by = "g=c('8','4')",
+    at = list(g = c("8", "4")),
+    allowed = list(g = c("8", "4"))
+  )
+)
+
+# identify the slope in a label like "4, 1", "6" or "g4 a1"; variables that
+# the label leaves out take the one value their filter allows
+.slope_key <- function(label, allowed) {
+  vars <- names(allowed)
+  tokens <- strsplit(trimws(label), ",\\s*|\\s+")[[1]]
+  prefixed <- vapply(
+    tokens,
+    function(tk) {
+      any(vapply(vars, function(v) startsWith(tk, v) && substring(tk, nchar(v) + 1) %in% allowed[[v]], logical(1)))
+    },
+    logical(1)
+  )
+  values <- list()
+  if (all(prefixed)) {
+    for (tk in tokens) {
+      v <- vars[vapply(vars, function(v) startsWith(tk, v), logical(1))][1]
+      values[[v]] <- substring(tk, nchar(v) + 1)
+    }
+  } else if (length(tokens) == length(vars)) {
+    values <- stats::setNames(as.list(tokens), vars)
+  } else {
+    values <- stats::setNames(as.list(tokens), vars[lengths(allowed) > 1])
+  }
+  for (v in vars) {
+    if (is.null(values[[v]])) {
+      values[[v]] <- allowed[[v]]
+    }
+  }
+  paste(vars, unlist(values[vars]), sep = "=", collapse = ";")
+}
+
+# check that `out` and the emmeans contrasts in `oracle` pair one to one,
+# with the same two slopes, the same difference and the same p-value
+expect_matches_oracle <- function(out, oracle, allowed) {
+  oracle <- as.data.frame(oracle)
+  expect_identical(nrow(out), nrow(oracle))
+  oracle_pairs <- strsplit(as.character(oracle$contrast), " - ", fixed = TRUE)
+  oracle_1 <- vapply(oracle_pairs, function(x) .slope_key(x[1], allowed), character(1))
+  oracle_2 <- vapply(oracle_pairs, function(x) .slope_key(x[2], allowed), character(1))
+  used <- logical(nrow(oracle))
+  for (i in seq_len(nrow(out))) {
+    key_1 <- .slope_key(as.character(out$Level1[i]), allowed)
+    key_2 <- .slope_key(as.character(out$Level2[i]), allowed)
+    same <- which(oracle_1 == key_1 & oracle_2 == key_2)
+    flipped <- which(oracle_1 == key_2 & oracle_2 == key_1)
+    expect_identical(length(same) + length(flipped), 1L)
+    if (length(same) + length(flipped) != 1L) {
+      next
+    }
+    j <- c(same, flipped)
+    sign <- if (length(same)) 1 else -1
+    expect_false(used[j])
+    used[j] <- TRUE
+    expect_equal(out$Difference[i], sign * oracle$estimate[j], tolerance = 1e-5, ignore_attr = TRUE)
+    expect_equal(out$p[i], oracle$p.value[j], tolerance = 1e-5, ignore_attr = TRUE)
+  }
+  expect_true(all(used))
+}
+
+.oracle <- function(model, filter) {
+  specs <- if (length(filter$allowed) == 1) ~g else ~ g * a
+  pairs(emmeans::emtrends(model, specs, "wt", at = filter$at), adjust = "none")
+}
+
+probes <- list(
+  list(model = m3, filter = filters$F1),
+  list(model = m3, filter = filters$F2),
+  list(model = m3, filter = filters$F3),
+  list(model = m6, filter = filters$F1),
+  list(model = m6, filter = filters$F2),
+  list(model = m6, filter = filters$F3),
+  list(model = m2, filter = filters$F4)
+)
+
+
+test_that("estimate_contrasts, slopes, by-filters, marginaleffects", {
+  for (probe in probes) {
+    oracle <- .oracle(probe$model, probe$filter)
+    for (estimate in c("typical", "average")) {
+      out <- estimate_contrasts(probe$model, "wt", by = probe$filter$by, estimate = estimate)
+      expect_matches_oracle(out, oracle, probe$filter$allowed)
+    }
+  }
+})
+
+
+test_that("estimate_contrasts, slopes, by-filters, emmeans", {
+  for (probe in probes) {
+    oracle <- .oracle(probe$model, probe$filter)
+    out <- estimate_contrasts(probe$model, "wt", by = probe$filter$by, backend = "emmeans")
+    expect_matches_oracle(out, oracle, probe$filter$allowed)
+  }
+})
+
+
+test_that("estimate_contrasts, slopes, by-filters leave fewer than two slopes", {
+  calls <- list(
+    list(model = m2, by = "g='4'"),
+    list(model = m3, by = c("g='4'", "a='1'"))
+  )
+  for (cl in calls) {
+    expect_error(
+      estimate_contrasts(cl$model, "wt", by = cl$by, estimate = "typical"),
+      regexp = "at least two"
+    )
+    expect_error(
+      estimate_contrasts(cl$model, "wt", by = cl$by, estimate = "average"),
+      regexp = "at least two"
+    )
+    expect_error(
+      estimate_contrasts(cl$model, "wt", by = cl$by, backend = "emmeans"),
+      regexp = "at least two"
+    )
+  }
+  # spaces around "=", and `estimate = "specific"`
+  expect_error(
+    estimate_contrasts(m3, "wt", by = c("g = '4'", "a = '1'"), estimate = "specific"),
+    regexp = "at least two"
+  )
+  out <- estimate_contrasts(m3, "wt", by = c("g = '4'", "a"), estimate = "average")
+  expect_matches_oracle(out, .oracle(m3, filters$F1), filters$F1$allowed)
+  # with `estimate = "average"`, only slopes with data count: for 5 gears,
+  # all cars have am = 1
+  mt$gear <- factor(mt$gear)
+  m7 <- lm(mpg ~ wt * a + gear, data = mt)
+  expect_error(
+    estimate_contrasts(m7, "wt", by = c("gear='5'", "a"), estimate = "average"),
+    regexp = "at least two"
+  )
+  expect_identical(nrow(estimate_contrasts(m7, "wt", by = c("gear='5'", "a"))), 1L)
+})
+
+
+test_that("estimate_contrasts, slopes, by-filter values that are not levels", {
+  oracle <- .oracle(m3, filters$F1)
+  out <- estimate_contrasts(m3, "wt", by = c("g=c('4','5')", "a"), estimate = "average")
+  expect_matches_oracle(out, oracle, filters$F1$allowed)
+  out <- estimate_contrasts(m3, "wt", by = c("g=c('4','5')", "a"), backend = "emmeans")
+  expect_matches_oracle(out, oracle, filters$F1$allowed)
+})
+
+
+test_that("estimate_contrasts, slopes, emmeans, numeric by or factors from the formula", {
+  # emmeans does not filter these, but still returns all contrasts
+  m_factor <- lm(mpg ~ wt * factor(cyl), data = mtcars)
+  out <- estimate_contrasts(m_factor, "wt", by = "cyl=c(4,6)", backend = "emmeans")
+  expect_identical(nrow(out), 3L)
+  expect_true(all(is.finite(out$Difference)))
+  m_numeric <- lm(mpg ~ wt * hp * g, data = mt)
+  out <- estimate_contrasts(m_numeric, "wt", by = c("g='4'", "hp"), backend = "emmeans")
+  expect_identical(nrow(out), 3L)
+  expect_true(all(is.finite(out$Difference)))
+})
+
+
+test_that("estimate_contrasts, slopes, no by-filter", {
+  oracle <- pairs(emmeans::emtrends(m3, ~ g * a, "wt"), adjust = "none")
+  allowed <- list(g = c("4", "6", "8"), a = c("0", "1"))
+  for (estimate in c("typical", "average")) {
+    out <- estimate_contrasts(m3, "wt", by = c("g", "a"), estimate = estimate)
+    expect_matches_oracle(out, oracle, allowed)
+  }
+  out <- estimate_contrasts(m3, "wt", by = c("g", "a"), backend = "emmeans")
+  expect_matches_oracle(out, oracle, allowed)
+})