diff --git a/DESCRIPTION b/DESCRIPTION
index a26d86c11..8b4a377ff 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.14.0.9
+Version: 0.14.0.10
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index 62f0a21d1..e383f61f7 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -170,6 +170,8 @@
 #' between-effects, or the difference of average slopes) at different levels
 #' of another variable are required, add that variable to the `contrast`
 #' argument instead, e.g. `contrast = c("x_within", "x_between", "factor")`.
+#' These contrasts can additionally be stratified by another variable using `by`
+#' again, e.g. `contrast = c("x_within", "x_between", "factor1"), by = "factor2"`.
 #' Note that when average slopes are contrasted, the `comparison` argument has
 #' no effect and is always set to `"pairwise"`. See also 'Examples'.
 #'
diff --git a/R/get_contexteffects.R b/R/get_contexteffects.R
index 98b9aa978..b36037e19 100644
--- a/R/get_contexteffects.R
+++ b/R/get_contexteffects.R
@@ -13,8 +13,16 @@
   # contrasts at each group level
   if (is.null(my_args$by)) {
     comparison <- stats::as.formula("~I(diff(x))")
+  } else if (length(my_args$by) > 2) {
+    # it is not possible to have more than two by-variables for now
+    insight::format_error(
+      "It is not possible to have more than two variables in `by` when calculating contrasts of slopes."
+    )
   } else {
-    comparison <- stats::as.formula(paste("~I(diff(x)) |", my_args$by))
+    comparison <- stats::as.formula(paste(
+      "~I(diff(x)) |",
+      paste(my_args$by, collapse = "+")
+    ))
   }
 
   # prepare arguments
@@ -47,10 +55,20 @@
 
   # pairwise comparison of context effects?
   if (identical(my_args$comparison, "context_pairwise")) {
-    # save original levels for formatting
-    original_levels <- .safe(out[[my_args$by]])
+    # if we have more than one by-variable, use the last one for stratification
+    stratify_by <- my_args$by[length(my_args$by)]
+    # save original levels from the contrast-variable for formatting
+    original_levels <- .safe(out[[my_args$by[1]]])
+    if (length(my_args$by) > 1) {
+      comparison <- stats::as.formula(paste(
+        "~pairwise |",
+        paste(stratify_by, collapse = "+")
+      ))
+    } else {
+      comparison <- ~pairwise
+    }
     # calculate pairwise comparisons
-    out <- marginaleffects::hypotheses(out, hypothesis = ~pairwise)
+    out <- marginaleffects::hypotheses(out, hypothesis = comparison)
     # format comparison levels. we first split the column with "b" parameter
     # names, like "(b1) - (b2)", into two columns. Then we remove the "b" and
     # just keep the number, which indicates the row.
@@ -68,8 +86,8 @@
       params$Level2 <- original_levels[params$Level2]
       # save attributes
       att <- attributes(out)
-      # remove old "by" column and bind new one back to output
-      out[[my_args$by]] <- NULL
+      # remove old "by" column and bind new one with the contrast-levels back to output
+      out[[my_args$by[1]]] <- NULL
       out <- cbind(params, out)
       # add back original attributes
       cn <- colnames(out)
diff --git a/R/get_marginalcontrasts.R b/R/get_marginalcontrasts.R
index 1735ce916..4f4dc04f7 100644
--- a/R/get_marginalcontrasts.R
+++ b/R/get_marginalcontrasts.R
@@ -358,11 +358,14 @@ get_marginalcontrasts <- function(
     }
     # overwrite some of the previous arguments
     context_effects <- TRUE
-    # if we have no "by" variable, user doesn't want to stratify, so set to
-    # pairwise comparisons of categorical variable
-    if (is.null(my_args$by) && length(my_args$contrast) > 2) {
+    if (length(my_args$contrast) > 2) {
+      # if we have more than just the slopes in "contrasts", the user wants to
+      # contrast by a second grouping variable - we want pairwise comparisons
       comparison <- "context_pairwise"
-      my_args$by <- my_args$contrast[3:length(my_args$contrast)]
+      original_by <- my_args$by <- c(
+        my_args$contrast[3:length(my_args$contrast)],
+        my_args$by
+      )
       my_args$contrast <- my_args$contrast[1:2]
     } else {
       comparison <- "context"
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index 71a322ac7..9ff31fd6d 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -420,6 +420,8 @@ pairwise comparisons of the difference of between within- and
 between-effects, or the difference of average slopes) at different levels
 of another variable are required, add that variable to the \code{contrast}
 argument instead, e.g. \code{contrast = c("x_within", "x_between", "factor")}.
+These contrasts can additionally be stratified by another variable using \code{by}
+again, e.g. \verb{contrast = c("x_within", "x_between", "factor1"), by = "factor2"}.
 Note that when average slopes are contrasted, the \code{comparison} argument has
 no effect and is always set to \code{"pairwise"}. See also 'Examples'.
 }
diff --git a/tests/testthat/test-estimate_contrasts_context.R b/tests/testthat/test-estimate_contrasts_context.R
index ec3455989..3d05a5ab3 100644
--- a/tests/testthat/test-estimate_contrasts_context.R
+++ b/tests/testthat/test-estimate_contrasts_context.R
@@ -47,6 +47,54 @@ test_that("estimate_contrast, context effects, linear", {
     tolerance = 1e-4,
     ignore_attr = TRUE
   )
+
+  skip_on_os(c("mac", "linux"))
+  m <- lm(bill_dep ~ sex * year * (bill_len_between + bill_len_within), data = d)
+  out <- estimate_contrasts(
+    m,
+    c("bill_len_between", "bill_len_within"),
+    by = c("sex", "year")
+  )
+  expect_identical(
+    capture.output(out),
+    c(
+      "Marginal Contrasts Analysis",
+      "",
+      "sex    | year | Difference |   SE |        95% CI |    z |      p",
+      "-----------------------------------------------------------------",
+      "female | 2007 |       0.28 | 0.08 | [ 0.12, 0.45] | 3.41 | < .001",
+      "female | 2008 |       0.26 | 0.06 | [ 0.15, 0.37] | 4.47 | < .001",
+      "female | 2009 |       0.24 | 0.10 | [ 0.05, 0.42] | 2.48 |  0.013",
+      "male   | 2007 |       0.46 | 0.10 | [ 0.26, 0.65] | 4.60 | < .001",
+      "male   | 2008 |       0.28 | 0.06 | [ 0.15, 0.40] | 4.41 | < .001",
+      "male   | 2009 |       0.10 | 0.10 | [-0.11, 0.30] | 0.93 |  0.355",
+      "",
+      "Variable predicted: bill_dep",
+      "Predictors contrasted: bill_len_between, bill_len_within",
+      "p-values are uncorrected."
+    )
+  )
+  out <- estimate_contrasts(
+    m,
+    c("bill_len_between", "bill_len_within", "sex"),
+    by = "year"
+  )
+  expect_identical(
+    capture.output(out),
+    c(
+      "Marginal Contrasts Analysis",
+      "",
+      "Level1 | Level2 | year | Difference |   SE |        95% CI |     z |     p",
+      "--------------------------------------------------------------------------",
+      "male   | female | 2007 |       0.17 | 0.13 | [-0.08, 0.43] |  1.33 | 0.183",
+      "male   | female | 2008 |       0.02 | 0.09 | [-0.15, 0.18] |  0.18 | 0.853",
+      "male   | female | 2009 |      -0.14 | 0.14 | [-0.42, 0.14] | -0.99 | 0.320",
+      "",
+      "Variable predicted: bill_dep",
+      "Predictors contrasted: bill_len_between, bill_len_within",
+      "p-values are uncorrected."
+    )
+  )
 })
 
 test_that("estimate_contrast, context effects, glm", {