diff --git a/DESCRIPTION b/DESCRIPTION
index 0eaf096c9..81465f10d 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.2
+Version: 0.14.0.3
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
diff --git a/NEWS.md b/NEWS.md
index 8f00f3bd3..374d1c669 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -9,6 +9,10 @@
 
 * Support for models of class `nestedLogit`.
 
+* Added option  `comparison = "slope"` to `estimate_contrast()`, to allow
+  calculating contrasts of average slopes. This can also be useful when
+  estimating the context effects of within- and between-effects in a model.
+
 # modelbased 0.14.0
 
 ## Changes
diff --git a/R/clean_names.R b/R/clean_names.R
index 94b494691..b0850b6f1 100644
--- a/R/clean_names.R
+++ b/R/clean_names.R
@@ -1,10 +1,9 @@
 # Clean names -------------------------------------------------------------
 
-
 #' @keywords internal
 .clean_names_frequentist <- function(means, predict = NULL, info = NULL) {
-  names(means)[names(means) == "emmean"] <- .guess_estimate_name(predict, info)
-  names(means)[names(means) == "response"] <- .guess_estimate_name(predict, info)
+  names(means)[names(means) == "emmean"] <- .guess_estimate_name(predict, info = info)
+  names(means)[names(means) == "response"] <- .guess_estimate_name(predict, info = info)
   names(means)[names(means) == "prob"] <- "Probability"
   names(means)[names(means) == "estimate"] <- "Difference"
   names(means)[names(means) == "odds.ratio"] <- "Odds_ratio"
diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index b7af196ea..6e8941eed 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -33,8 +33,12 @@
 #'   [this website](https://marginaleffects.com/bonus/hypothesis.html) and
 #'   section _Comparison options_ below.
 #'   * String: One of `"pairwise"`, `"reference"`, `"sequential"`, `"meandev"`
-#'     `"meanotherdev"`, `"poly"`, `"helmert"`, or `"trt_vs_ctrl"`. To test
-#'     multiple hypotheses jointly (usually used for factorial designs),
+#'     `"meanotherdev"`, `"poly"`, `"helmert"`, `"slope"` or `"trt_vs_ctrl"`.
+#'     The `"slope"` option calculates contrasts between average slopes and can
+#'     also be used to calculate "context" effects, which is the difference of
+#'     within- and between-effects (see
+#'     https://statisticalhorizons.com/between-within-contextual-effects/). To
+#'     test multiple hypotheses jointly (usually used for factorial designs),
 #'     `comparison` can also be `"joint"`. In this case, use the `test` argument
 #'     to specify which test should be conducted: `"F"` (default) or `"Chi2"`.
 #'   * String: Special string options are `"inequality"`, `"inequality_ratio"`,
diff --git a/R/format.R b/R/format.R
index dcf9e25ce..350067b86 100644
--- a/R/format.R
+++ b/R/format.R
@@ -14,9 +14,16 @@ format.estimate_contrasts <- function(
     select <- NULL
     include_grid <- FALSE
   }
+  # change parameter name for context effects
+  if (isTRUE(attributes(x)$context_effects)) {
+    x$Parameter <- "Average slope"
+  }
+
   # don't print columns of adjusted_for variables
   adjusted_for <- attr(x, "adjusted_for", exact = TRUE)
-  if (!is.null(adjusted_for) && all(adjusted_for %in% colnames(x)) && !isTRUE(include_grid)) {
+  if (
+    !is.null(adjusted_for) && all(adjusted_for %in% colnames(x)) && !isTRUE(include_grid)
+  ) {
     # remove non-focal terms from data frame
     x[adjusted_for] <- NULL
   } else if (isTRUE(include_grid)) {
@@ -84,7 +91,13 @@ format.estimate_contrasts <- function(
   }
 
   if (!is.null(format) && format %in% c("md", "markdown", "html")) {
-    insight::format_table(x, ci_brackets = c("(", ")"), select = select, format = format, ...)
+    insight::format_table(
+      x,
+      ci_brackets = c("(", ")"),
+      select = select,
+      format = format,
+      ...
+    )
   } else {
     insight::format_table(x, select = select, ...)
   }
@@ -106,8 +119,12 @@ format.estimate_grouplevel <- format.estimate_contrasts
 #' @export
 format.estimate_smooth <- function(x, ...) {
   # Colnames
-  if ("Size" %in% names(x)) x$Size <- ifelse(x$Size < 1, paste0(insight::format_value(x$Size * 100), "%"), "100%")
-  if ("Part" %in% names(x)) x$Part <- insight::format_value(x$Part, protect_integers = TRUE)
+  if ("Size" %in% names(x)) {
+    x$Size <- ifelse(x$Size < 1, paste0(insight::format_value(x$Size * 100), "%"), "100%")
+  }
+  if ("Part" %in% names(x)) {
+    x$Part <- insight::format_value(x$Part, protect_integers = TRUE)
+  }
 
   insight::format_table(x, ...)
 }
@@ -137,6 +154,7 @@ format.marginaleffects_means <- function(x, model, ci = 0.95, ...) {
   }
   non_focal <- setdiff(colnames(model_data), attr(x, "focal_terms"))
   predict_type <- attributes(x)$predict
+  transform <- attributes(x)$transform
 
   # special attributes we get from "get_marginalcontrasts()"
   comparison <- list(...)$hypothesis
@@ -161,7 +179,7 @@ format.marginaleffects_means <- function(x, model, ci = 0.95, ...) {
     # for simple means, we don't want p-values
     remove_columns <- c(remove_columns, "p")
     # estimate name
-    estimate_name <- .guess_estimate_name(predict_type, info)
+    estimate_name <- .guess_estimate_name(predict_type, transform, info)
   }
 
   # reshape and format columns
@@ -190,7 +208,14 @@ format.marginaleffects_slopes <- function(x, model, ci = 0.95, ...) {
   }
   model_data <- insight::get_data(model, verbose = FALSE)
   # define all columns that should be removed
-  remove_columns <- c("Predicted", "s.value", "S", "CI", "rowid_dedup", equivalence_columns)
+  remove_columns <- c(
+    "Predicted",
+    "s.value",
+    "S",
+    "CI",
+    "rowid_dedup",
+    equivalence_columns
+  )
   # for contrasting slope, we need to keep the "Parameter" column
   # however, for estimating trends/slope, the "Parameter" column is usually
   # redundant. Since we cannot check for class-attributes, we simply check if
@@ -377,7 +402,11 @@ format.marginaleffects_contrasts <- function(
       # replace all comparison levels with tokens
       params[] <- lapply(params, function(comparison_pair) {
         for (j in seq_along(all_num_levels)) {
-          comparison_pair <- sub(all_num_levels[j], replace_num_levels[j], comparison_pair)
+          comparison_pair <- sub(
+            all_num_levels[j],
+            replace_num_levels[j],
+            comparison_pair
+          )
         }
         for (j in seq_along(all_levels)) {
           comparison_pair <- sub(
@@ -489,7 +518,10 @@ format.marginaleffects_contrasts <- function(
       if (!is.null(contrast_filter)) {
         # make sure we also have all levels for non-filtered variables
         contrast_filter <- insight::compact_list(c(
-          lapply(dgrid[setdiff(focal_terms, unique(c(by, names(contrast_filter))))], unique),
+          lapply(
+            dgrid[setdiff(focal_terms, unique(c(by, names(contrast_filter))))],
+            unique
+          ),
           contrast_filter
         ))
         # now create combinations of all filter variables
@@ -516,7 +548,6 @@ format.marginaleffects_contrasts <- function(
 # Helper ----------------------------------------------------------------------
 # -----------------------------------------------------------------------------
 
-
 # since we combine levels from different factors, we have to make
 # sure levels are unique across different terms. If not, paste
 # variable names to levels. We first find the intersection of all
@@ -553,15 +584,17 @@ equivalence_columns <- c(
 # outputs from {marginaleffects}
 
 #' @keywords internal
-.standardize_marginaleffects_columns <- function(x,
-                                                 remove_columns,
-                                                 model,
-                                                 model_data,
-                                                 info,
-                                                 ci = 0.95,
-                                                 estimate_name = NULL,
-                                                 is_contrast_analysis = FALSE,
-                                                 ...) {
+.standardize_marginaleffects_columns <- function(
+  x,
+  remove_columns,
+  model,
+  model_data,
+  info,
+  ci = 0.95,
+  estimate_name = NULL,
+  is_contrast_analysis = FALSE,
+  ...
+) {
   # tidy output - we want to tidy the output, using `model_parameters()` or
   # `describe_posterior()` for Bayesian models. We also need to know how the
   # coefficient column is named, because we replace that column name with an
@@ -578,7 +611,12 @@ equivalence_columns <- c(
   # column names for their "coefficient". We now extract the relevant one.
   possible_colnames <- c(
     attributes(params)$coefficient_name,
-    "Coefficient", "Slope", "Predicted", "Median", "Mean", "MAP"
+    "Coefficient",
+    "Slope",
+    "Predicted",
+    "Median",
+    "Mean",
+    "MAP"
   )
   coefficient_name <- intersect(possible_colnames, colnames(params))[1]
   # we need to remove some more columns
@@ -665,9 +703,18 @@ equivalence_columns <- c(
   if (.is_inequality_comparison(comparison_hypothesis)) {
     # fix for pairwise inequality labels - these are named like "(b1) - (b2)" etc.
     # but we want the original labels instead of b1, b2 etc.
-    if(comparison_hypothesis %in% c("inequality_pairwise", "inequality_ratio_pairwise") && !is.null(by_terms)) {
+    if (
+      comparison_hypothesis %in%
+        c("inequality_pairwise", "inequality_ratio_pairwise") &&
+        !is.null(by_terms)
+    ) {
       # clean parameter names
-      parameter_names <- gsub(")", "", gsub("(", "", params$Parameter, fixed = TRUE), fixed = TRUE)
+      parameter_names <- gsub(
+        ")",
+        "",
+        gsub("(", "", params$Parameter, fixed = TRUE),
+        fixed = TRUE
+      )
       # extract data for by-variable
       by_var <- model_data[[by_terms]]
       # make sure we have a factor
@@ -693,7 +740,11 @@ equivalence_columns <- c(
     }
 
     # fix labels for inequality analysis for slopes
-    if (comparison_hypothesis %in% c("inequality", "inequality_ratio") && isTRUE(attributes(x)$compute_slopes)) {
+    if (
+      comparison_hypothesis %in%
+        c("inequality", "inequality_ratio") &&
+        isTRUE(attributes(x)$compute_slopes)
+    ) {
       # for slopes, we either have the trend variable, or only the grouping,
       # but not the "inequality" variabe (the first in "by"). Update labels,
       # so users know by which variables slopes are averaged and grouped
@@ -752,7 +803,9 @@ equivalence_columns <- c(
 
 #' @keywords internal
 .add_contrasts_ci <- function(is_contrast_analysis, params) {
-  if (is_contrast_analysis && !"CI_low" %in% colnames(params) && "SE" %in% colnames(params)) {
+  if (
+    is_contrast_analysis && !"CI_low" %in% colnames(params) && "SE" %in% colnames(params)
+  ) {
     # extract ci-level
     if ("CI" %in% colnames(params)) {
       ci <- params[["CI"]][1]
@@ -782,17 +835,33 @@ equivalence_columns <- c(
 # based on on which scale predictions were requested
 
 #' @keywords internal
-.guess_estimate_name <- function(predict_type, info) {
+.guess_estimate_name <- function(predict_type, transform = NULL, info) {
   # estimate name
   if (is.null(predict_type)) {
     estimate_name <- "Mean"
-  } else if (!is.null(predict_type) && tolower(predict_type) %in% .brms_aux_elements()) {
+  } else if (tolower(predict_type) %in% .brms_aux_elements()) {
     # for Bayesian models with distributional parameter
     estimate_name <- tools::toTitleCase(predict_type)
-  } else if (!predict_type %in% c("none", "link") && (info$is_binomial || info$is_bernoulli || info$is_multinomial)) {
+  } else if (
+    !predict_type %in% c("none", "link") &&
+      (info$is_binomial || info$is_bernoulli || info$is_multinomial)
+  ) {
     # here we add all models that model the probability of an outcome, such as
     # binomial, multinomial, or Bernoulli models
     estimate_name <- "Probability"
+  } else if (
+    predict_type %in%
+      c("none", "link") &&
+      identical(transform, "exp") &&
+      (info$is_binomial || info$is_bernoulli || info$is_multinomial)
+  ) {
+    # here we add all models that have odds ratios as exponentiated coefficients
+    estimate_name <- "Odds_Ratio"
+  } else if (
+    predict_type %in% c("none", "link") && identical(transform, "exp") && (info$is_count)
+  ) {
+    # here we add all models that have IRRs as exponentiated coefficients
+    estimate_name <- "IRR"
   } else if (predict_type == "survival" && info$is_survival) {
     # this is for survival models, where we want to predict the survival probability
     estimate_name <- "Probability"
@@ -800,7 +869,11 @@ equivalence_columns <- c(
     # this is for zero-inflated models, where we want to predict the probability
     # of a zero-inflated outcome
     estimate_name <- "Probability"
-  } else if (predict_type %in% c("response", "invlink(link)") && (info$is_beta || info$is_orderedbeta)) {
+  } else if (
+    predict_type %in%
+      c("response", "invlink(link)") &&
+      (info$is_beta || info$is_orderedbeta)
+  ) {
     # this is for beta regression models, where we want to predict the mean
     # value of the outcome, which is a proportion
     estimate_name <- "Proportion"
@@ -834,7 +907,11 @@ equivalence_columns <- c(
       if (substring(input_string, match_positions[i], match_positions[i]) == "-") {
         inside_parentheses <- FALSE
         for (j in seq_along(match_positions)) {
-          if (i != j && match_positions[i] > match_positions[j] && match_positions[i] < (match_positions[j] + match_lengths[j])) {
+          if (
+            i != j &&
+              match_positions[i] > match_positions[j] &&
+              match_positions[i] < (match_positions[j] + match_lengths[j])
+          ) {
             inside_parentheses <- TRUE
             break
           }
@@ -850,11 +927,7 @@ equivalence_columns <- c(
     for (i in 1:(length(split_positions) - 1)) {
       parts <- c(
         parts,
-        substring(
-          input_string,
-          split_positions[i] + 1,
-          split_positions[i + 1] - 1
-        )
+        substring(input_string, split_positions[i] + 1, split_positions[i + 1] - 1)
       )
     }
   }
diff --git a/R/get_contexteffects.R b/R/get_contexteffects.R
new file mode 100644
index 000000000..09a46dc73
--- /dev/null
+++ b/R/get_contexteffects.R
@@ -0,0 +1,43 @@
+# special contrasts: context effects ----------------------------------------
+# ---------------------------------------------------------------------------
+
+.get_contexteffects <- function(
+  model,
+  my_args,
+  predict = NULL,
+  transform = NULL,
+  model_info,
+  ...
+) {
+  if (model_info$is_linear) {
+    out <- marginaleffects::avg_comparisons(
+      model,
+      variables = my_args$contrast,
+      hypothesis = my_args$comparison,
+      ...
+    )
+  } else {
+    dots <- list(...)
+    fun_args <- list(model, variables = my_args$contrast, hypothesis = my_args$comparison)
+    # set default for "type" argument, if not provided
+    if (is.null(predict)) {
+      fun_args$type <- "link"
+      # if "type" was not provided, also change transform argument. we do
+      # this only when user did not provide "type", else - if user provided
+      # "type" - we keep the default NULL
+      if (is.null(transform)) {
+        fun_args$transform <- "exp"
+      }
+    } else {
+      fun_args$type <- predict
+      fun_args$transform <- transform
+    }
+    out <- do.call(marginaleffects::avg_comparisons, c(fun_args, dots))
+  }
+  # save some labels for printing
+  attr(out, "by") <- my_args$by
+  attr(out, "contrast") <- my_args$contrast
+  attr(out, "context_effects") <- TRUE
+  class(out) <- unique(c("marginaleffects_means", class(out)))
+  out
+}
diff --git a/R/get_inequalitycontrasts.R b/R/get_inequalitycontrasts.R
index 9f534fab1..a2203cef7 100644
--- a/R/get_inequalitycontrasts.R
+++ b/R/get_inequalitycontrasts.R
@@ -1,7 +1,7 @@
 # special contrasts: inequality ---------------------------------------------
 # ---------------------------------------------------------------------------
 
-get_inequalitycontrasts <- function(
+.get_inequalitycontrasts <- function(
   model,
   model_data,
   my_args,
@@ -93,7 +93,10 @@ get_inequalitycontrasts <- function(
 
     # to calculate marginal effects inequalities, all contrast predictors
     # must be factors
-    check_factors <- .safe(vapply(model_data[my_args$contrast], is.factor, logical(1)), NULL)
+    check_factors <- .safe(
+      vapply(model_data[my_args$contrast], is.factor, logical(1)),
+      NULL
+    )
     if (is.null(check_factors) || !all(check_factors)) {
       insight::format_error(
         "All variables specified in `contrast` must be factors for `comparison = \"inequality\"`."
@@ -183,7 +186,10 @@ get_inequalitycontrasts <- function(
     inequality_pairwise = ,
     inequality = stats::as.formula(paste(c("~ pairwise", group), collapse = " | ")),
     inequality_ratio_pairwise = ,
-    inequality_ratio = stats::as.formula(paste(c("ratio ~ pairwise", group), collapse = " | ")),
+    inequality_ratio = stats::as.formula(paste(
+      c("ratio ~ pairwise", group),
+      collapse = " | "
+    )),
   )
   f2 <- stats::as.formula(paste(c("~ I(mean(abs(x)))", group), collapse = " | "))
   list(f1 = f1, f2 = f2)
@@ -200,7 +206,9 @@ get_inequalitycontrasts <- function(
   f <- unlist(strsplit(insight::safe_deparse(comparison), "|", fixed = TRUE))
   # check parts left and right of the bar "|"
   if (length(f) != 2) {
-    insight::format_error("Formula must contain exactly one `|` character separating two parts, e.g. `~ inequality | group`.")
+    insight::format_error(
+      "Formula must contain exactly one `|` character separating two parts, e.g. `~ inequality | group`."
+    )
   }
   # check parts left and right of the bar "|"
   left_part <- insight::trim_ws(f[[1]])
@@ -252,6 +260,10 @@ get_inequalitycontrasts <- function(
       return(out)
     }
   }
+  # handle special value: contrasting average slopes (context effects)
+  if (identical(comparison, "slope")) {
+    comparison <- "context"
+  }
   comparison
 }
 
@@ -268,7 +280,11 @@ get_inequalitycontrasts <- function(
   )
 
   if (!is.null(comparison)) {
-    if (length(comparison) == 1 && is.character(comparison) && comparison %in% inequality_comparisons) {
+    if (
+      length(comparison) == 1 &&
+        is.character(comparison) &&
+        comparison %in% inequality_comparisons
+    ) {
       return(TRUE)
     }
     # if we have a formula, we check whether it starts with "inequality". We
@@ -287,7 +303,7 @@ get_inequalitycontrasts <- function(
 
 # return the valid inequality comparison value
 .inequality_type <- function(comparison) {
-  if (!.is_inequality_comparison(comparison)){
+  if (!.is_inequality_comparison(comparison)) {
     return(NULL)
   }
   comparison
diff --git a/R/get_marginalcontrasts.R b/R/get_marginalcontrasts.R
index d97c73ad7..b727be51a 100644
--- a/R/get_marginalcontrasts.R
+++ b/R/get_marginalcontrasts.R
@@ -90,7 +90,7 @@ get_marginalcontrasts <- function(
     # Inequality and Total Effect Summary Measures for Nominal and Ordinal Variables
     # Sociological Science February 5, 10.15195/v12.a7
     # this requires a special handling, because we can only use it with avg_comparisons
-    out <- get_inequalitycontrasts(
+    out <- .get_inequalitycontrasts(
       model,
       model_data,
       my_args,
@@ -101,11 +101,20 @@ get_marginalcontrasts <- function(
       ...
     )
     predict <- "response"
+  } else if (isTRUE(my_args$context_effects)) {
+    out <- .get_contexteffects(model, my_args, predict, transform, model_info, ...)
+    # set defaults, for proper printing
+    if (is.null(predict)) {
+      predict <- "link"
+      if (is.null(transform)) {
+        transform <- "exp"
+      }
+    }
   } else if (compute_slopes) {
     # sanity check - contrast for slopes only makes sense when we have a "by" argument
     if (is.null(my_args$by)) {
       insight::format_error(
-        "Please specify the `by` argument to calculate contrasts of slopes."
+        "Please specify the `by` argument to calculate contrasts of slopes. If you want to calculate the average contrast between two slopes, use `comparison = \"slope\"` instead."
       )
     }
     # call slopes with hypothesis argument
@@ -156,10 +165,12 @@ get_marginalcontrasts <- function(
     info = list(
       contrast = my_args$contrast,
       predict = predict,
+      transform = transform,
       comparison = my_args$comparison,
       estimate = estimate,
       p_adjust = p_adjust,
       contrast_filter = my_args$contrast_filter,
+      context_effects = my_args$context_effects,
       keep_iterations = keep_iterations
     )
   )
@@ -170,11 +181,15 @@ get_marginalcontrasts <- function(
     out <- .p_adjust(model, out, p_adjust, verbose, ...)
   }
 
-  # remove "estimate_means" class attribute
-  class(out) <- setdiff(
-    unique(c("marginaleffects_contrasts", class(out))),
-    "estimate_means"
-  )
+  # no extra class attribute for context effects, because we don't want
+  # the regular contrast formatting here.
+  if (!isTRUE(my_args$context_effects)) {
+    # remove "estimate_means" class attribute
+    class(out) <- setdiff(
+      unique(c("marginaleffects_contrasts", class(out))),
+      "estimate_means"
+    )
+  }
   out
 }
 
@@ -232,6 +247,12 @@ get_marginalcontrasts <- function(
   # init
   comparison_slopes <- by_filter <- contrast_filter <- by_token <- NULL
   joint_test <- FALSE
+  context_effects <- FALSE
+  # overwrite "comparison" when it's set to "context".
+  if (identical(comparison, "context")) {
+    comparison <- "b1 - b2 = 0"
+    context_effects <- TRUE
+  }
   # save original `by`
   original_by <- my_args$by
   original_comparison <- comparison
@@ -404,6 +425,8 @@ get_marginalcontrasts <- function(
       contrast_filter = insight::compact_list(contrast_filter),
       # in case we have a joint/omnibus test
       joint_test = joint_test,
+      # remember if we want to calculate context effects
+      context_effects = context_effects,
       # cleaned `by` and `contrast`, without filtering information
       cleaned_by = gsub("=.*", "\\1", my_args$by),
       cleaned_contrast = gsub("=.*", "\\1", my_args$contrast)
diff --git a/R/get_marginalmeans.R b/R/get_marginalmeans.R
index d98a1eacb..238c02ac6 100644
--- a/R/get_marginalmeans.R
+++ b/R/get_marginalmeans.R
@@ -287,11 +287,19 @@ get_marginalmeans <- function(
 #   - `dots`: The `...` arguments, with arguments consumed by this function removed.
 .get_datagrid_means <- function(model, by, estimate, dots) {
   dg_factors <- switch(estimate, specific = "reference", "all")
-  dg_args <- list(model, by = by, factors = dg_factors, include_random = TRUE, verbose = FALSE)
+  dg_args <- list(
+    model,
+    by = by,
+    factors = dg_factors,
+    include_random = TRUE,
+    verbose = FALSE
+  )
   # did user request weights? These are not supported for data-grid
   # marginalization types
   if (
-    estimate %in% c("specific", "typical") && (!is.null(dots$weights) || !is.null(dots$wts))
+    estimate %in%
+      c("specific", "typical") &&
+      (!is.null(dots$weights) || !is.null(dots$wts))
   ) {
     insight::format_warning(
       "Using weights is not possible when `estimate` is set to \"typical\" or \"specific\". Use `estimate = \"average\"` to include weights for marginal means or contrasts."
@@ -322,7 +330,10 @@ get_marginalmeans <- function(
   # restore data types -  if we have defined numbers in `by`, like
   # `by = "predictor = 5"`, and `predictor` was a factor, it is returned as
   # numeric in the data grid. Fix this here, else marginal effects will fail
-  datagrid <- datawizard::data_restoretype(datagrid, insight::get_data(model, verbose = FALSE))
+  datagrid <- datawizard::data_restoretype(
+    datagrid,
+    insight::get_data(model, verbose = FALSE)
+  )
 
   list(datagrid = datagrid, datagrid_info = datagrid_info, dots = dots)
 }
@@ -420,7 +431,8 @@ get_marginalmeans <- function(
   # pass data grid in such situations - but we still created the data grid based
   # on the `by` variables, for internal use, for example filtering at this point
   if (
-    identical(estimate, "average") && all(datagrid_info$at_specs$varname %in% colnames(means))
+    identical(estimate, "average") &&
+      all(datagrid_info$at_specs$varname %in% colnames(means))
   ) {
     # sanity check - are all filter values from the data grid in the marginaleffects
     # object? For `estimate_average()`, predictions are based on the data, not
@@ -519,7 +531,7 @@ get_marginalmeans <- function(
     "at", "by", "focal_terms", "adjusted_for", "predict", "trend", "comparison",
     "contrast", "estimate", "p_adjust", "transform", "datagrid", "preserve_range",
     "coef_name", "slope", "ci", "model_info", "contrast_filter",
-    "keep_iterations", "joint_test", "vcov", "equivalence"
+    "keep_iterations", "joint_test", "vcov", "equivalence", "context_effects"
   )
 }
 
@@ -531,9 +543,16 @@ get_marginalmeans <- function(
   model,
   by = NULL,
   contrast = NULL,
+  comparison = NULL,
   verbose = TRUE,
   ...
 ) {
+  # special case: calculating context effects for models with within-between
+  # effects. In this case, we don't want any further checks
+  if (identical(comparison, "context")) {
+    return(list(by = by, contrast = contrast))
+  }
+
   # Gather info and data from model
   model_data <- insight::get_data(model, verbose = FALSE)
 
@@ -597,7 +616,10 @@ get_marginalmeans <- function(
         "Model contains an offset-term and you average predictions over the distribution of that offset. If you want to fix the offset to a specific value, for instance `1`, use `offset = 1` and set `estimate = \"typical\"`."
       )
       # if offset term is log-transformed, tell user. offset should be fixed then
-      log_offset <- insight::find_transformation(insight::find_offset(model, as_term = TRUE))
+      log_offset <- insight::find_transformation(insight::find_offset(
+        model,
+        as_term = TRUE
+      ))
       if (!is.null(log_offset) && startsWith(log_offset, "log")) {
         msg <- c(
           msg,
diff --git a/R/table_footer.R b/R/table_footer.R
index 255e60786..ce1aa7da8 100644
--- a/R/table_footer.R
+++ b/R/table_footer.R
@@ -22,14 +22,23 @@
   if (isTRUE(info$joint_test)) {
     table_footer <- NULL
   } else {
-    table_footer <- paste0("\nVariable predicted: ", toString(insight::find_response(model)))
+    table_footer <- paste0(
+      "\nVariable predicted: ",
+      toString(insight::find_response(model))
+    )
   }
 
   # modulated predictors (focal terms) ---------------------------------------
 
   if (!is.null(by) && !isTRUE(info$joint_test)) {
     modulate_string <- switch(type, inequality = , contrasts = "contrasted", "modulated")
-    table_footer <- paste0(table_footer, "\nPredictors ", modulate_string, ": ", toString(by))
+    table_footer <- paste0(
+      table_footer,
+      "\nPredictors ",
+      modulate_string,
+      ": ",
+      toString(by)
+    )
   }
 
   # predictors controlled (non-focal terms) ----------------------------------
@@ -49,7 +58,11 @@
       # over the list, because we may have different types of data
       for (av in seq_along(adjusted_values)) {
         if (is.numeric(adjusted_values[[av]])) {
-          adjusted_for[av] <- sprintf("%s (%.2g)", adjusted_for[av], adjusted_values[[av]])
+          adjusted_for[av] <- sprintf(
+            "%s (%.2g)",
+            adjusted_for[av],
+            adjusted_values[[av]]
+          )
         } else if (identical(type, "predictions")) {
           adjusted_for[av] <- sprintf("%s (%s)", adjusted_for[av], adjusted_values[[av]])
         }
@@ -87,7 +100,12 @@
 
   # tell user about scale of predictions / contrasts -------------------------
 
-  result_type <- switch(type, inequality = "Differences", contrasts = "Contrasts", "Predictions")
+  result_type <- switch(
+    type,
+    inequality = "Differences",
+    contrasts = "Contrasts",
+    "Predictions"
+  )
 
   if (!is.null(predict) && isFALSE(model_info$is_linear)) {
     # exceptions
@@ -141,7 +159,11 @@
       hypothesis_labels <- unlist(
         lapply(parameter_names, function(i) {
           rows <- as.numeric(sub(".", "", i))
-          paste0(i, " = ", toString(paste0(info$focal_terms, " [", transposed_dg[, rows], "]")))
+          paste0(
+            i,
+            " = ",
+            toString(paste0(info$focal_terms, " [", transposed_dg[, rows], "]"))
+          )
         }),
         use.names = FALSE
       )
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index b7611889f..48551a78b 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -122,8 +122,12 @@ described below, see documentation of \link[marginaleffects:comparisons]{margina
 section \emph{Comparison options} below.
 \itemize{
 \item String: One of \code{"pairwise"}, \code{"reference"}, \code{"sequential"}, \code{"meandev"}
-\code{"meanotherdev"}, \code{"poly"}, \code{"helmert"}, or \code{"trt_vs_ctrl"}. To test
-multiple hypotheses jointly (usually used for factorial designs),
+\code{"meanotherdev"}, \code{"poly"}, \code{"helmert"}, \code{"slope"} or \code{"trt_vs_ctrl"}.
+The \code{"slope"} option calculates contrasts between average slopes and can
+also be used to calculate "context" effects, which is the difference of
+within- and between-effects (see
+https://statisticalhorizons.com/between-within-contextual-effects/). To
+test multiple hypotheses jointly (usually used for factorial designs),
 \code{comparison} can also be \code{"joint"}. In this case, use the \code{test} argument
 to specify which test should be conducted: \code{"F"} (default) or \code{"Chi2"}.
 \item String: Special string options are \code{"inequality"}, \code{"inequality_ratio"},
diff --git a/man/get_emmeans.Rd b/man/get_emmeans.Rd
index cee4e34bc..db81057d1 100644
--- a/man/get_emmeans.Rd
+++ b/man/get_emmeans.Rd
@@ -140,8 +140,12 @@ described below, see documentation of \link[marginaleffects:comparisons]{margina
 section \emph{Comparison options} below.
 \itemize{
 \item String: One of \code{"pairwise"}, \code{"reference"}, \code{"sequential"}, \code{"meandev"}
-\code{"meanotherdev"}, \code{"poly"}, \code{"helmert"}, or \code{"trt_vs_ctrl"}. To test
-multiple hypotheses jointly (usually used for factorial designs),
+\code{"meanotherdev"}, \code{"poly"}, \code{"helmert"}, \code{"slope"} or \code{"trt_vs_ctrl"}.
+The \code{"slope"} option calculates contrasts between average slopes and can
+also be used to calculate "context" effects, which is the difference of
+within- and between-effects (see
+https://statisticalhorizons.com/between-within-contextual-effects/). To
+test multiple hypotheses jointly (usually used for factorial designs),
 \code{comparison} can also be \code{"joint"}. In this case, use the \code{test} argument
 to specify which test should be conducted: \code{"F"} (default) or \code{"Chi2"}.
 \item String: Special string options are \code{"inequality"}, \code{"inequality_ratio"},
diff --git a/tests/testthat/test-attributes_estimatefun.R b/tests/testthat/test-attributes_estimatefun.R
index 10a2d2c41..487d15143 100644
--- a/tests/testthat/test-attributes_estimatefun.R
+++ b/tests/testthat/test-attributes_estimatefun.R
@@ -43,7 +43,11 @@ test_that("attributes_means, contrasts", {
       "comparison", "contrast", "transform", "keep_iterations", "joint_test"
     )
   )
-  estim <- suppressMessages(estimate_contrasts(model, "Species", backend = "marginaleffects"))
+  estim <- suppressMessages(estimate_contrasts(
+    model,
+    "Species",
+    backend = "marginaleffects"
+  ))
   # fmt: skip
   expect_named(
     attributes(estim),
@@ -52,7 +56,7 @@ test_that("attributes_means, contrasts", {
       "model", "response", "ci", "p_adjust", "backend", "call", "focal_terms",
       "adjusted_for", "predict", "comparison", "contrast", "estimate",
       "transform", "datagrid", "preserve_range", "coef_name", "model_info",
-      "keep_iterations", "joint_test", "vcov"
+      "keep_iterations", "joint_test", "vcov", "context_effects"
     )
   )
   estim <- suppressMessages(estimate_contrasts(
@@ -68,7 +72,7 @@ test_that("attributes_means, contrasts", {
       "model", "response", "ci", "p_adjust", "backend", "call", "focal_terms",
       "adjusted_for", "predict", "comparison", "contrast", "estimate",
       "transform", "datagrid", "preserve_range", "coef_name", "model_info",
-      "keep_iterations", "joint_test", "vcov"
+      "keep_iterations", "joint_test", "vcov", "context_effects"
     )
   )
   estim <- suppressMessages(estimate_contrasts(
@@ -85,7 +89,7 @@ test_that("attributes_means, contrasts", {
       "model", "response", "ci", "p_adjust", "backend", "call", "focal_terms",
       "adjusted_for", "predict", "comparison", "contrast", "estimate",
       "transform", "datagrid", "preserve_range", "coef_name", "model_info",
-      "contrast_filter", "keep_iterations", "joint_test","vcov"
+      "contrast_filter", "keep_iterations", "joint_test","vcov", "context_effects"
     )
   )
 })
@@ -105,7 +109,11 @@ test_that("attributes_means, slopes", {
       "keep_iterations"
     )
   )
-  estim <- suppressMessages(estimate_slopes(model, "Sepal.Width", backend = "marginaleffects"))
+  estim <- suppressMessages(estimate_slopes(
+    model,
+    "Sepal.Width",
+    backend = "marginaleffects"
+  ))
   # fmt: skip
   expect_named(
     attributes(estim),
diff --git a/tests/testthat/test-estimate_contrasts_context.R b/tests/testthat/test-estimate_contrasts_context.R
new file mode 100644
index 000000000..9722283c7
--- /dev/null
+++ b/tests/testthat/test-estimate_contrasts_context.R
@@ -0,0 +1,54 @@
+skip_on_cran()
+skip_if_not_installed("marginaleffects", minimum_version = "0.29.0")
+skip_on_os("mac")
+skip_if(getRversion() < "4.5.0")
+skip_if_not_installed("datawizard")
+
+test_that("estimate_contrast, context effects, linear", {
+  data(penguins, package = "datasets")
+  d <- datawizard::demean(penguins, "bill_len", by = "species")
+  m <- lm(bill_dep ~ bill_len_between + bill_len_within, data = d)
+
+  b <- coef(summary(m))[2:3, 1]
+  se <- coef(summary(m))[2:3, 2]
+
+  out <- modelbased::estimate_contrasts(
+    m,
+    c("bill_len_between", "bill_len_within"),
+    comparison = "context"
+  )
+  expect_equal(out$Mean, b[1] - b[2], tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$SE, sqrt((se[1]^2 + se[2]^2)), tolerance = 1e-4, ignore_attr = TRUE)
+})
+
+test_that("estimate_contrast, context effects, glm", {
+  data(penguins, package = "datasets")
+  d <- datawizard::demean(penguins, "bill_len", by = "species")
+  d$out <- datawizard::categorize(d$flipper_len) - 1
+  m <- glm(out ~ bill_len_between + bill_len_within, data = d, family = binomial())
+
+  b <- coef(summary(m))[2:3, 1]
+  se <- coef(summary(m))[2:3, 2]
+
+  out <- modelbased::estimate_contrasts(
+    m,
+    c("bill_len_between", "bill_len_within"),
+    comparison = "context"
+  )
+  expect_equal(out$Odds_Ratio, exp(b[1] - b[2]), tolerance = 1e-4, ignore_attr = TRUE)
+
+  out <- modelbased::estimate_contrasts(
+    m,
+    c("bill_len_between", "bill_len_within"),
+    comparison = "slope"
+  )
+  expect_equal(out$Odds_Ratio, exp(b[1] - b[2]), tolerance = 1e-4, ignore_attr = TRUE)
+
+  out <- modelbased::estimate_contrasts(
+    m,
+    c("bill_len_between", "bill_len_within"),
+    comparison = "context",
+    predict = "response"
+  )
+  expect_equal(out$Probability, 0.01784138, tolerance = 1e-4, ignore_attr = TRUE)
+})
diff --git a/tests/testthat/test-keep_iterations.R b/tests/testthat/test-keep_iterations.R
index 261266138..28589247e 100644
--- a/tests/testthat/test-keep_iterations.R
+++ b/tests/testthat/test-keep_iterations.R
@@ -69,7 +69,7 @@ test_that("estimate_contrasts() - posterior draws", {
       "model", "response", "ci", "p_adjust", "backend", "call", "focal_terms",
       "adjusted_for", "predict", "comparison", "contrast", "estimate",
       "transform", "datagrid", "preserve_range", "coef_name", "model_info",
-      "keep_iterations", "joint_test"
+      "keep_iterations", "joint_test", "context_effects"
     )
   )
   # fmt: skip
@@ -95,7 +95,7 @@ test_that("estimate_contrasts() - posterior draws", {
       "model", "response", "ci", "p_adjust", "backend", "call", "focal_terms",
       "adjusted_for", "predict", "comparison", "contrast", "estimate",
       "transform", "datagrid", "preserve_range", "coef_name", "model_info",
-      "keep_iterations", "joint_test"
+      "keep_iterations", "joint_test", "context_effects"
     )
   )
   # fmt: skip
@@ -196,7 +196,12 @@ test_that("estimate_means() - posterior draws, emmeans", {
 test_that("estimate_contrasts() - posterior draws, emmeans", {
   m <- insight::download_model("brms_1")
   skip_if(is.null(m))
-  out <- estimate_contrasts(m, by = "wt=c(3,4,5)", keep_iterations = 5, backend = "emmeans")
+  out <- estimate_contrasts(
+    m,
+    by = "wt=c(3,4,5)",
+    keep_iterations = 5,
+    backend = "emmeans"
+  )
   # fmt: skip
   expect_named(
     attributes(out),
@@ -217,7 +222,12 @@ test_that("estimate_contrasts() - posterior draws, emmeans", {
   )
   expect_identical(dim(out), c(3L, 12L))
 
-  out <- estimate_contrasts(m, by = "wt=c(3,4,5)", keep_iterations = TRUE, backend = "emmeans")
+  out <- estimate_contrasts(
+    m,
+    by = "wt=c(3,4,5)",
+    keep_iterations = TRUE,
+    backend = "emmeans"
+  )
   # fmt: skip
   expect_named(
     attributes(out),