diff --git a/DESCRIPTION b/DESCRIPTION
index 057c7ab23..3069ae7a8 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
+Version: 0.18.0.1
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
diff --git a/NEWS.md b/NEWS.md
index d497a6dca..704d59924 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,37 @@
+# modelbased (devel)
+
+## Bug fixes
+
+* `p_adjust = "tukey"` and `p_adjust = "sidak"` now give the same p-values as
+  *emmeans* in `estimate_contrasts()` and `estimate_slopes()` with
+  `backend = "marginaleffects"`:
+  * For comparisons that are not pairwise, for example
+    `comparison = "reference"` or `"sequential"`, `"tukey"` is replaced by
+    `"sidak"` with a message. Before, the Tukey adjustment was applied.
+  * For slopes tested against a value other than 0, for example
+    `hypothesis = 1`, and for slopes that are not compared pairwise, `"tukey"`
+    is replaced by `"sidak"` with a message. Slopes compared pairwise, for
+    example with `comparison = ~pairwise` or `post_process = ~pairwise`, keep
+    the Tukey adjustment.
+  * The Tukey adjustment of slopes compared pairwise with `estimate_slopes()`
+    now counts the slopes, not the comparisons. Before, 4 slopes were counted
+    as 6, the number of comparisons, which made the p-values too large. With a
+    group term, for example `comparison = ~pairwise | group`, the slopes and
+    the comparisons are counted in each group.
+  * The Sidak adjustment of contrasts now counts the comparisons in each group
+    of the `by` variables. Before, it counted the levels of the contrast
+    variable, so for 4 levels and 6 pairwise comparisons, the p-values were
+    too small. A joint test, or another output with one test in each group,
+    is therefore no longer changed by `"sidak"`.
+  * A partly matched `hypothesis` argument, for example `hypo = ~pairwise`, is
+    now treated as `hypothesis`.
+
+* `estimate_slopes()` with `backend = "emmeans"` now applies `p_adjust` for the
+  `stats::p.adjust()` methods, `"tukey"` and `"sidak"`. Before, it returned the
+  unadjusted p-values for every method. As in `emmeans::test()`, `"tukey"` is
+  replaced by `"sidak"` with a message, because the slopes are tested
+  against 0.
+
 # modelbased 0.18.0
 
 ## Breaking Changes
diff --git a/R/estimate_contrast_methods.R b/R/estimate_contrast_methods.R
index d71f7e8e6..8e16877e1 100644
--- a/R/estimate_contrast_methods.R
+++ b/R/estimate_contrast_methods.R
@@ -152,7 +152,9 @@ estimate_contrasts.estimate_predicted <- function(
     }
   }
 
-  # p-value adjustment?
+  # p-value adjustment? Formatting drops the "by" attribute, which defines the
+  # families of the Sidak adjustment
+  attr(out, "by") <- by
   out <- .p_adjust(model, out, p_adjust, verbose, ...)
 
   # Table formatting
diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index 38c9bb143..836c2fbda 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -26,7 +26,15 @@
 #' Plagborg-Møller, 2019). Details for the other options can be found in the
 #' p-value adjustment section of the `emmeans::test` documentation or
 #' `?stats::p.adjust`. Note that certain options provided by the **emmeans**
-#' package are only available if you set `backend = "emmeans"`.
+#' package are only available if you set `backend = "emmeans"`. The Tukey
+#' adjustment is only appropriate for pairwise comparisons. With
+#' `backend = "marginaleffects"`, `"tukey"` is therefore replaced by `"sidak"`
+#' (with a message) for comparisons that are not pairwise, for example
+#' `comparison = "reference"` or `"sequential"`, as **emmeans** does. The Sidak
+#' adjustment counts the comparisons within each group of the `by` variables
+#' that are columns of the output. When the output has no such column, for
+#' example for contrasts of slopes without a group term in `comparison`, it
+#' counts all comparisons.
 #' @param comparison Specify the type of contrasts or tests that should be
 #' carried out. See also section *Comparison options* below for details.
 #' * When `backend = "emmeans"`, can be one of `"pairwise"`, `"revpairwise"`,
diff --git a/R/estimate_slopes.R b/R/estimate_slopes.R
index 12b1e7646..97c34c602 100644
--- a/R/estimate_slopes.R
+++ b/R/estimate_slopes.R
@@ -26,19 +26,25 @@
 #' `slope` from the `by` argument. `trend` is an alias for `slope`, for backward
 #' compatibility. See 'Examples'.
 #' @param p_adjust The p-values adjustment method for frequentist multiple
-#' comparisons. For `estimate_slopes()`, multiple comparison only occurs for
-#' Johnson-Neyman intervals, i.e. in case of interactions with two numeric
-#' predictors (one specified in `slope`, one in `by`). In this case, the
-#' `"esarey"` or `"sup-t"` options are recommended, but `p_adjust` can also be
-#' one of `"none"` (default), `"hochberg"`, `"hommel"`, `"bonferroni"`, `"BH"`,
+#' comparisons, when several slopes are tested, or when slopes are compared.
+#' For Johnson-Neyman intervals, i.e. in case of interactions with two numeric
+#' predictors (one specified in `slope`, one in `by`), the `"esarey"` or
+#' `"sup-t"` options are recommended, but `p_adjust` can also be one of
+#' `"none"` (default), `"hochberg"`, `"hommel"`, `"bonferroni"`, `"BH"`,
 #' `"BY"`, `"fdr"`, `"tukey"`, `"sidak"`, or `"holm"`. `"sup-t"` computes
 #' simultaneous confidence bands, also called sup-t confidence band (Montiel
 #' Olea & Plagborg-Møller, 2019). The Tukey adjustment is only appropriate for
-#' pairwise comparisons. When **marginaleffects** computes slopes that are each
-#' tested against 0, `"tukey"` is therefore replaced by `"sidak"` (with a
-#' message), as `emmeans::test()` does. This replacement is not made when
-#' slopes are compared, for example with the `comparison` or `post_process`
-#' argument.
+#' pairwise comparisons. It is used when slopes are compared pairwise, for
+#' example with `comparison = ~pairwise` or `post_process = ~pairwise`. For
+#' slopes that are each tested against 0 or against another value (for
+#' example `hypothesis = 1`), and for slopes that are not compared pairwise,
+#' `"tukey"` is replaced by `"sidak"` (with a message), as `emmeans::test()`
+#' does. The Sidak adjustment counts the tests (rows) and the Tukey adjustment
+#' the compared slopes, in the whole output or, for comparisons with a group
+#' term (for example `comparison = ~pairwise | group`), in each group. In
+#' `estimate_slopes()` with `backend = "emmeans"`, only the `stats::p.adjust()`
+#' methods, `"tukey"` and `"sidak"` are applied to the slopes, and other
+#' methods are ignored.
 #' @inheritParams estimate_means
 #' @inheritParams parameters::model_parameters.default
 #'
@@ -190,6 +196,31 @@ estimate_slopes <- function(
       verbose = verbose,
       ...
     )
+    # p-value adjustment, done by `parameters::model_parameters()` for the
+    # "adjust" setting of the emmeans object. Not for Bayesian models, which
+    # have no p-values, and not for ordinal models.
+    emmeans_adjust <- c(stats::p.adjust.methods, "tukey", "sidak")
+    emmeans_method <- emmeans_adjust[match(tolower(p_adjust), tolower(emmeans_adjust))]
+    if (
+      length(emmeans_method) == 1 &&
+        !is.na(emmeans_method) &&
+        emmeans_method != "none" &&
+        is.null(attributes(estimated)$ordinal_mode) &&
+        !insight::model_info(model, response = 1)$is_bayesian
+    ) {
+      # slopes are tested against 0, so we use Sidak instead of Tukey, as
+      # `emmeans::test()` does
+      if (emmeans_method == "tukey") {
+        if (verbose && nrow(estimated@grid) > 1) {
+          insight::format_alert(
+            "The Tukey adjustment is only appropriate for pairwise comparisons. Slopes are tested against 0, so p-values were adjusted with the Sidak method instead."
+          )
+        }
+        emmeans_method <- "sidak"
+      }
+      estimated <- stats::update(estimated, adjust = emmeans_method)
+      attr(estimated, "p_adjust") <- emmeans_method
+    }
     trends <- .format_emmeans_slopes(model, estimated, ci, ...)
   } else {
     # Marginaleffects --------------------------------------------------------
diff --git a/R/p_adjust.R b/R/p_adjust.R
index a8f409f72..ccfc82b65 100644
--- a/R/p_adjust.R
+++ b/R/p_adjust.R
@@ -56,24 +56,106 @@
   }
 
   # slopes from `get_marginaltrends()` are each tested against 0, unless
-  # `comparison` (or `hypothesis`) or `post_process` is used. In that case,
-  # the rows are comparisons of slopes, and we keep Tukey. Else, we use Sidak,
-  # as emmeans does, because the Tukey adjustment is only for one set of
-  # pairwise comparisons. With one slope, or without p-values, the code below
-  # says that no adjustment was made.
+  # `comparison` (or `hypothesis`, also when partly matched, like `hypo`) or
+  # `post_process` is used. The Tukey adjustment is only for one set of
+  # pairwise comparisons, so we keep Tukey only for pairwise comparisons of
+  # slopes. Else, we use Sidak, as emmeans does. With one slope, or without
+  # p-values, the code below says that no adjustment was made.
   dots <- list(...)
+  n_slopes <- NULL
+  slope_comparison <- NULL
+  # number of tests in the family of each row. Comparisons of slopes with a
+  # group term (like `~pairwise | group`) have one family per group, else
+  # all rows are one family
+  slope_family <- rep(nrow(params), nrow(params))
+  if (is_unformatted) {
+    dots_names <- names(dots)
+    if (is.null(dots_names)) {
+      dots_names <- character(0)
+    }
+    # names that are a prefix of "hypothesis", like `hypo`, which R matches
+    # to the `hypothesis` argument of `avg_slopes()`
+    hypothesis_arg <- dots_names[nzchar(dots_names) & startsWith("hypothesis", dots_names)]
+    # `post_process` is applied last, so it defines the rows of the output
+    slope_comparison <- c(
+      dots["post_process"],
+      dots["comparison"],
+      dots[hypothesis_arg]
+    )
+    slope_comparison <- Filter(Negate(is.null), slope_comparison)
+    if (length(slope_comparison)) {
+      slope_family <- .p_adjust_family_size(
+        params,
+        by = .comparison_group_terms(slope_comparison[[1]])
+      )
+    }
+  }
+  if (p_adjust == "tukey" && is_unformatted && "p.value" %in% colnames(params)) {
+    reason <- NULL
+    if (!length(slope_comparison)) {
+      if (nrow(params) > 1) {
+        reason <- "Slopes are tested against 0"
+      }
+    } else if (is.numeric(slope_comparison[[1]]) && is.null(dim(slope_comparison[[1]]))) {
+      reason <- paste0(
+        "Slopes are tested against ",
+        toString(insight::format_value(slope_comparison[[1]], protect_integers = TRUE))
+      )
+    } else if (.is_pairwise_comparison(slope_comparison[[1]])) {
+      # for k slopes, we have k * (k - 1) / 2 pairwise comparisons, in each
+      # group. Without the group columns, we cannot count them.
+      group_terms <- .comparison_group_terms(slope_comparison[[1]])
+      n_slopes <- (1 + sqrt(1 + 8 * slope_family)) / 2
+      if (
+        !all(group_terms %in% colnames(params)) ||
+          anyNA(n_slopes) ||
+          any(n_slopes != round(n_slopes))
+      ) {
+        n_slopes <- NULL
+        reason <- "The number of compared slopes could not be found"
+      }
+    } else {
+      reason <- "Slopes are not compared pairwise"
+    }
+    if (!is.null(reason)) {
+      # with one test per family, Sidak does not change the p-values, so we
+      # say nothing
+      if (verbose && max(slope_family) > 1) {
+        insight::format_alert(paste0(
+          "The Tukey adjustment is only appropriate for pairwise comparisons. ",
+          reason,
+          ", so p-values were adjusted with the Sidak method instead."
+        ))
+      }
+      p_adjust <- "sidak"
+      attr(params, "p_adjust") <- "sidak"
+    }
+  }
+
+  # contrasts (formatted output): the Tukey adjustment is only appropriate for
+  # pairwise comparisons. For other comparisons, like "reference" or
+  # "sequential", we use Sidak, as emmeans does.
+  comparison <- attributes(params)$comparison
+  # number of contrasts in the family of each row: the rows within each group
+  # of the `by` variables, or of the group term of `comparison`
+  if (!is_unformatted) {
+    contrast_family <- .p_adjust_family_size(
+      params,
+      by = c(attributes(params)$by, .comparison_group_terms(comparison))
+    )
+  }
   if (
     p_adjust == "tukey" &&
-      is_unformatted &&
-      nrow(params) > 1 &&
-      "p.value" %in% colnames(params) &&
-      is.null(dots[["comparison"]]) &&
-      is.null(dots[["hypothesis"]]) &&
-      is.null(dots[["post_process"]])
+      !is_unformatted &&
+      !is.null(comparison) &&
+      !.is_pairwise_comparison(comparison) &&
+      "p" %in% colnames(params)
   ) {
-    if (verbose) {
+    # with one contrast per family, Sidak does not change the p-values, so we
+    # say nothing
+    if (verbose && max(contrast_family) > 1) {
       insight::format_alert(
-        "The Tukey adjustment is only appropriate for pairwise comparisons. Slopes are tested against 0, so p-values were adjusted with the Sidak method instead."
+        "The Tukey adjustment is only appropriate for pairwise comparisons, so p-values were adjusted with the Sidak method instead."
       )
     }
     p_adjust <- "sidak"
@@ -92,12 +174,26 @@
 
   # needed for rank adjustment
   focal_terms <- .safe(datagrid[focal])
-  if (is_unformatted) {
+  if (!is.null(n_slopes)) {
+    # pairwise comparisons of slopes
+    rank_adjust <- n_slopes
+  } else if (is_unformatted) {
     rank_adjust <- nrow(params)
   } else if (is.null(focal_terms)) {
     rank_adjust <- 1
   } else {
     rank_adjust <- prod(vapply(focal_terms, insight::n_unique, numeric(1)))
+    # pairwise contrasts: k levels give k * (k - 1) / 2 contrasts in each
+    # family. This also counts the levels when a group term or filter makes
+    # the focal terms differ from the compared levels.
+    n_levels <- (1 + sqrt(1 + 8 * contrast_family)) / 2
+    if (
+      .is_pairwise_comparison(comparison) &&
+        !anyNA(n_levels) &&
+        all(n_levels == round(n_levels))
+    ) {
+      rank_adjust <- n_levels
+    }
   }
 
   p_column <- if (is_unformatted) "p.value" else "p"
@@ -125,7 +221,7 @@
       statistic <- NULL
     }
     if (!is.null(statistic)) {
-      if (rank_adjust < 2) {
+      if (min(rank_adjust) < 2) {
         if (verbose) {
           insight::format_alert(
             "Tukey adjustment requires at least 2 groups. P-values were not adjusted."
@@ -150,14 +246,67 @@
       insight::format_alert("No test-statistic found. P-values were not adjusted.")
     }
   } else if (p_adjust == "sidak") {
-    # sidak adjustment
-    params[[p_column]] <- 1 - (1 - params[[p_column]])^rank_adjust
+    # sidak adjustment. The family is the set of tests: for slopes, all rows,
+    # and for contrasts, the rows within each group of the `by` variables
+    # (like emmeans).
+    if (is_unformatted) {
+      family_size <- slope_family
+    } else {
+      family_size <- contrast_family
+    }
+    # 1 - (1 - p)^n, without losing precision for small p-values
+    params[[p_column]] <- -expm1(family_size * log1p(-params[[p_column]]))
   }
 
   params
 }
 
 
+# TRUE if a comparison is a formula whose right side is `pairwise` or
+# `revpairwise`, with or without a `| group` term, e.g. `difference ~ pairwise`
+.is_pairwise_comparison <- function(comparison) {
+  if (is.character(comparison)) {
+    return(length(comparison) == 1 && comparison %in% c("pairwise", "revpairwise"))
+  }
+  if (!inherits(comparison, "formula")) {
+    return(FALSE)
+  }
+  rhs <- comparison[[length(comparison)]]
+  if (is.call(rhs) && identical(rhs[[1]], as.name("|"))) {
+    rhs <- rhs[[2]]
+  }
+  is.name(rhs) && as.character(rhs) %in% c("pairwise", "revpairwise")
+}
+
+
+# variables of the group term of a comparison formula, like `h` in
+# `~pairwise | h`, or NULL
+.comparison_group_terms <- function(comparison) {
+  if (!inherits(comparison, "formula")) {
+    return(NULL)
+  }
+  rhs <- comparison[[length(comparison)]]
+  if (is.call(rhs) && identical(rhs[[1]], as.name("|"))) {
+    return(all.vars(rhs[[3]]))
+  }
+  NULL
+}
+
+
+# number of tests in the family of each row: the rows with the same values of
+# the `by` variables (by default, the `by` attribute)
+.p_adjust_family_size <- function(params, by = attributes(params)$by) {
+  by <- .safe(insight::trim_ws(gsub("=.*", "", by)))
+  by <- intersect(by, colnames(params))
+  if (!length(by)) {
+    return(rep(nrow(params), nrow(params)))
+  }
+  # missing values of a `by` variable form their own group
+  groups <- interaction(lapply(params[by], addNA), drop = TRUE)
+  as.vector(table(groups)[as.integer(groups)])
+}
+
+
 .p_adjust_supt <- function(params, model) {
   insight::check_if_installed("mvtnorm")
   # get correlation matrix, based on the covariance matrix
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index bf67bb755..0ed291730 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -280,7 +280,15 @@ confidence bands, also called sup-t confidence band (Montiel Olea &
 Plagborg-Møller, 2019). Details for the other options can be found in the
 p-value adjustment section of the \code{emmeans::test} documentation or
 \code{?stats::p.adjust}. Note that certain options provided by the \strong{emmeans}
-package are only available if you set \code{backend = "emmeans"}.}
+package are only available if you set \code{backend = "emmeans"}. The Tukey
+adjustment is only appropriate for pairwise comparisons. With
+\code{backend = "marginaleffects"}, \code{"tukey"} is therefore replaced by \code{"sidak"}
+(with a message) for comparisons that are not pairwise, for example
+\code{comparison = "reference"} or \code{"sequential"}, as \strong{emmeans} does. The Sidak
+adjustment counts the comparisons within each group of the \code{by} variables
+that are columns of the output. When the output has no such column, for
+example for contrasts of slopes without a group term in \code{comparison}, it
+counts all comparisons.}
 
 \item{transform}{A function applied to predictions and confidence intervals
 to (back-) transform results, which can be useful in case the regression
diff --git a/man/estimate_slopes.Rd b/man/estimate_slopes.Rd
index c766719c7..18c585d21 100644
--- a/man/estimate_slopes.Rd
+++ b/man/estimate_slopes.Rd
@@ -172,19 +172,25 @@ appropriate transformation-function. Note that no standard errors are returned
 when transformations are applied.}
 
 \item{p_adjust}{The p-values adjustment method for frequentist multiple
-comparisons. For \code{estimate_slopes()}, multiple comparison only occurs for
-Johnson-Neyman intervals, i.e. in case of interactions with two numeric
-predictors (one specified in \code{slope}, one in \code{by}). In this case, the
-\code{"esarey"} or \code{"sup-t"} options are recommended, but \code{p_adjust} can also be
-one of \code{"none"} (default), \code{"hochberg"}, \code{"hommel"}, \code{"bonferroni"}, \code{"BH"},
+comparisons, when several slopes are tested, or when slopes are compared.
+For Johnson-Neyman intervals, i.e. in case of interactions with two numeric
+predictors (one specified in \code{slope}, one in \code{by}), the \code{"esarey"} or
+\code{"sup-t"} options are recommended, but \code{p_adjust} can also be one of
+\code{"none"} (default), \code{"hochberg"}, \code{"hommel"}, \code{"bonferroni"}, \code{"BH"},
 \code{"BY"}, \code{"fdr"}, \code{"tukey"}, \code{"sidak"}, or \code{"holm"}. \code{"sup-t"} computes
 simultaneous confidence bands, also called sup-t confidence band (Montiel
 Olea & Plagborg-Møller, 2019). The Tukey adjustment is only appropriate for
-pairwise comparisons. When \strong{marginaleffects} computes slopes that are each
-tested against 0, \code{"tukey"} is therefore replaced by \code{"sidak"} (with a
-message), as \code{emmeans::test()} does. This replacement is not made when
-slopes are compared, for example with the \code{comparison} or \code{post_process}
-argument.}
+pairwise comparisons. It is used when slopes are compared pairwise, for
+example with \code{comparison = ~pairwise} or \code{post_process = ~pairwise}. For
+slopes that are each tested against 0 or against another value (for
+example \code{hypothesis = 1}), and for slopes that are not compared pairwise,
+\code{"tukey"} is replaced by \code{"sidak"} (with a message), as \code{emmeans::test()}
+does. The Sidak adjustment counts the tests (rows) and the Tukey adjustment
+the compared slopes, in the whole output or, for comparisons with a group
+term (for example \code{comparison = ~pairwise | group}), in each group. In
+\code{estimate_slopes()} with \code{backend = "emmeans"}, only the \code{stats::p.adjust()}
+methods, \code{"tukey"} and \code{"sidak"} are applied to the slopes, and other
+methods are ignored.}
 
 \item{keep_iterations}{If \code{TRUE}, will keep all iterations (draws) of
 bootstrapped or Bayesian models. They will be added as additional columns
diff --git a/man/get_emmeans.Rd b/man/get_emmeans.Rd
index fe43628cc..7fe429137 100644
--- a/man/get_emmeans.Rd
+++ b/man/get_emmeans.Rd
@@ -362,19 +362,25 @@ subsequent, multi-step comparisons. After the initial comparison in
 specified post-process tests. See 'Examples'.}
 
 \item{p_adjust}{The p-values adjustment method for frequentist multiple
-comparisons. For \code{estimate_slopes()}, multiple comparison only occurs for
-Johnson-Neyman intervals, i.e. in case of interactions with two numeric
-predictors (one specified in \code{slope}, one in \code{by}). In this case, the
-\code{"esarey"} or \code{"sup-t"} options are recommended, but \code{p_adjust} can also be
-one of \code{"none"} (default), \code{"hochberg"}, \code{"hommel"}, \code{"bonferroni"}, \code{"BH"},
+comparisons, when several slopes are tested, or when slopes are compared.
+For Johnson-Neyman intervals, i.e. in case of interactions with two numeric
+predictors (one specified in \code{slope}, one in \code{by}), the \code{"esarey"} or
+\code{"sup-t"} options are recommended, but \code{p_adjust} can also be one of
+\code{"none"} (default), \code{"hochberg"}, \code{"hommel"}, \code{"bonferroni"}, \code{"BH"},
 \code{"BY"}, \code{"fdr"}, \code{"tukey"}, \code{"sidak"}, or \code{"holm"}. \code{"sup-t"} computes
 simultaneous confidence bands, also called sup-t confidence band (Montiel
 Olea & Plagborg-Møller, 2019). The Tukey adjustment is only appropriate for
-pairwise comparisons. When \strong{marginaleffects} computes slopes that are each
-tested against 0, \code{"tukey"} is therefore replaced by \code{"sidak"} (with a
-message), as \code{emmeans::test()} does. This replacement is not made when
-slopes are compared, for example with the \code{comparison} or \code{post_process}
-argument.}
+pairwise comparisons. It is used when slopes are compared pairwise, for
+example with \code{comparison = ~pairwise} or \code{post_process = ~pairwise}. For
+slopes that are each tested against 0 or against another value (for
+example \code{hypothesis = 1}), and for slopes that are not compared pairwise,
+\code{"tukey"} is replaced by \code{"sidak"} (with a message), as \code{emmeans::test()}
+does. The Sidak adjustment counts the tests (rows) and the Tukey adjustment
+the compared slopes, in the whole output or, for comparisons with a group
+term (for example \code{comparison = ~pairwise | group}), in each group. In
+\code{estimate_slopes()} with \code{backend = "emmeans"}, only the \code{stats::p.adjust()}
+methods, \code{"tukey"} and \code{"sidak"} are applied to the slopes, and other
+methods are ignored.}
 
 \item{iterations}{This argument has two distinct effects, depending on the
 context. When \code{effectsize = "boot"}, it defines the number of bootstrap
diff --git a/tests/testthat/test-p_adjust.R b/tests/testthat/test-p_adjust.R
new file mode 100644
index 000000000..29b2da0c7
--- /dev/null
+++ b/tests/testthat/test-p_adjust.R
@@ -0,0 +1,318 @@
+skip_on_cran()
+skip_if_not_installed("emmeans")
+skip_if_not_installed("marginaleffects")
+
+# data and models for the p-value adjustment tests
+set.seed(12)
+n <- 200
+d_adjust <- data.frame(
+  g = factor(rep(letters[1:4], each = n / 4)),
+  h = factor(rep(c("u", "v"), n / 2)),
+  x = rnorm(n)
+)
+d_adjust$y <- 1 +
+  as.integer(d_adjust$g) * 0.3 * d_adjust$x +
+  0.4 * as.integer(d_adjust$g) +
+  rnorm(n)
+m1 <- lm(y ~ g * h, data = d_adjust)
+m2 <- lm(y ~ x * g, data = d_adjust)
+
+# p-values of emmeans, in the row order of modelbased. modelbased compares
+# "b - a" where emmeans compares "a - b", so rows are matched by the sorted
+# pair of levels, and by the levels of `by`, if any
+emmeans_p <- function(out, em, by = NULL) {
+  em <- as.data.frame(summary(em))
+  pairs_em <- strsplit(as.character(em$contrast), " - ", fixed = TRUE)
+  key_em <- vapply(pairs_em, function(i) paste(sort(gsub("[()]", "", i)), collapse = "|"), "")
+  if (all(c("Level1", "Level2") %in% colnames(out))) {
+    key_out <- vapply(
+      seq_len(nrow(out)),
+      function(i) paste(sort(c(as.character(out$Level1[i]), as.character(out$Level2[i]))), collapse = "|"),
+      ""
+    )
+  } else {
+    # compared slopes, labelled like "(b) - (a)"
+    pairs_out <- strsplit(as.character(out$Parameter), " - ", fixed = TRUE)
+    key_out <- vapply(pairs_out, function(i) paste(sort(gsub("[()]", "", i)), collapse = "|"), "")
+  }
+  if (!is.null(by)) {
+    key_em <- paste(key_em, em[[by]])
+    key_out <- paste(key_out, out[[by]])
+  }
+  em$p.value[match(key_out, key_em)]
+}
+
+# the message can wrap over lines
+sidak_message <- "(?s)(?=.*Tukey)(?=.*Sidak)"
+
+
+test_that("p_adjust tukey, pairwise contrasts, marginaleffects", {
+  out <- estimate_contrasts(m1, "g", p_adjust = "tukey")
+  em <- suppressMessages(emmeans::emmeans(m1, "g"))
+  expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "tukey")), tolerance = 1e-5)
+  expect_identical(attributes(out)$p_adjust, "tukey")
+
+  out <- estimate_contrasts(m1, "g", by = "h", p_adjust = "tukey")
+  em <- suppressMessages(emmeans::emmeans(m1, "g", by = "h"))
+  expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "tukey"), by = "h"), tolerance = 1e-5)
+})
+
+
+test_that("p_adjust tukey, non-pairwise contrasts use sidak, marginaleffects", {
+  em <- suppressMessages(emmeans::emmeans(m1, "g"))
+  comparisons <- list("reference", ~reference, "sequential", "meandev")
+  methods <- c("trt.vs.ctrl", "trt.vs.ctrl", "consec", "eff")
+  for (i in seq_along(comparisons)) {
+    expect_message(
+      {
+        out <- estimate_contrasts(
+          m1,
+          "g",
+          comparison = comparisons[[i]],
+          p_adjust = "tukey"
+        )
+      },
+      regexp = sidak_message,
+      perl = TRUE
+    )
+    em_p <- as.data.frame(summary(emmeans::contrast(em, methods[i], adjust = "sidak")))$p.value
+    expect_equal(out$p, em_p, tolerance = 1e-5)
+    expect_identical(attributes(out)$p_adjust, "sidak")
+    expect_no_message(
+      estimate_contrasts(
+        m1,
+        "g",
+        comparison = comparisons[[i]],
+        p_adjust = "tukey",
+        verbose = FALSE
+      ),
+      message = "Tukey|Sidak"
+    )
+  }
+
+  # the example of the issue
+  mt <- mtcars
+  mt$gear <- factor(mt$gear)
+  m5 <- lm(mpg ~ gear, data = mt)
+  out <- suppressMessages(estimate_contrasts(
+    m5,
+    "gear",
+    comparison = "reference",
+    p_adjust = "tukey"
+  ))
+  em_p <- as.data.frame(summary(emmeans::contrast(
+    emmeans::emmeans(m5, "gear"),
+    "trt.vs.ctrl",
+    adjust = "sidak"
+  )))$p.value
+  expect_equal(out$p, em_p, tolerance = 1e-5)
+})
+
+
+test_that("p_adjust sidak, contrasts, family size, marginaleffects", {
+  out <- estimate_contrasts(m1, "g", p_adjust = "sidak")
+  em <- suppressMessages(emmeans::emmeans(m1, "g"))
+  expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "sidak")), tolerance = 1e-5)
+
+  out <- estimate_contrasts(m1, "g", by = "h", p_adjust = "sidak")
+  em <- suppressMessages(emmeans::emmeans(m1, "g", by = "h"))
+  expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "sidak"), by = "h"), tolerance = 1e-5)
+
+  em <- suppressMessages(emmeans::emmeans(m1, "g"))
+  comparisons <- list("reference", ~reference, "sequential", "meandev")
+  methods <- c("trt.vs.ctrl", "trt.vs.ctrl", "consec", "eff")
+  for (i in seq_along(comparisons)) {
+    expect_no_message(
+      {
+        out <- estimate_contrasts(
+          m1,
+          "g",
+          comparison = comparisons[[i]],
+          p_adjust = "sidak"
+        )
+      },
+      message = "Sidak"
+    )
+    em_p <- as.data.frame(summary(emmeans::contrast(em, methods[i], adjust = "sidak")))$p.value
+    expect_equal(out$p, em_p, tolerance = 1e-5)
+  }
+})
+
+
+test_that("p_adjust tukey and sidak, compared slopes, marginaleffects", {
+  em <- emmeans::emtrends(m2, "g", var = "x")
+  for (p_adjust in c("tukey", "sidak")) {
+    em_pairs <- pairs(em, adjust = p_adjust)
+    expect_no_message(
+      {
+        out <- estimate_slopes(m2, "x", by = "g", comparison = ~pairwise, p_adjust = p_adjust)
+      },
+      message = "Sidak"
+    )
+    expect_equal(out$p, emmeans_p(out, em_pairs), tolerance = 1e-5)
+    expect_no_message(
+      {
+        out <- estimate_contrasts(m2, contrast = "x", by = "g", p_adjust = p_adjust)
+      },
+      message = "Sidak"
+    )
+    expect_equal(out$p, emmeans_p(out, em_pairs), tolerance = 1e-5)
+  }
+
+  # post-processing: unadjusted p-values are based on z, so we only check the
+  # number of slopes in the Tukey adjustment, which uses the model's df
+  out <- suppressMessages(estimate_slopes(
+    m2,
+    "x",
+    by = "g",
+    post_process = ~pairwise,
+    p_adjust = "tukey"
+  ))
+  z <- intersect(c("z", "t", "statistic"), colnames(out))[1]
+  expect_equal(
+    as.vector(out$p),
+    stats::ptukey(sqrt(2) * abs(as.vector(out[[z]])), 4, stats::df.residual(m2), lower.tail = FALSE),
+    tolerance = 1e-8
+  )
+})
+
+
+test_that("p_adjust tukey, slopes tested against a value, marginaleffects", {
+  expect_message(
+    {
+      out <- estimate_slopes(m2, "x", by = "g", hypothesis = 1, p_adjust = "tukey")
+    },
+    regexp = sidak_message,
+    perl = TRUE
+  )
+  em_p <- as.data.frame(summary(emmeans::test(
+    emmeans::emtrends(m2, "g", var = "x"),
+    null = 1,
+    adjust = "sidak"
+  )))$p.value
+  expect_equal(out$p, em_p, tolerance = 1e-5)
+  expect_identical(attributes(out)$p_adjust, "sidak")
+
+  # partly matched `hypothesis` argument
+  em_pairs <- pairs(emmeans::emtrends(m2, "g", var = "x"), adjust = "tukey")
+  expect_no_message(
+    {
+      out <- estimate_slopes(m2, "x", by = "g", hypo = ~pairwise, p_adjust = "tukey")
+    },
+    message = "Sidak"
+  )
+  expect_equal(out$p, emmeans_p(out, em_pairs), tolerance = 1e-5)
+})
+
+
+test_that("p_adjust, slopes, emmeans backend", {
+  mt <- mtcars
+  mt$cyl <- factor(mt$cyl)
+  m4 <- lm(mpg ~ wt * cyl, data = mt)
+  em <- emmeans::emtrends(m4, "cyl", var = "wt")
+  for (p_adjust in c("holm", "bonferroni", "sidak", "tukey", "none")) {
+    out <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = p_adjust)
+    em_p <- suppressMessages(as.data.frame(summary(emmeans::test(em, adjust = p_adjust)))$p.value)
+    expect_equal(out$p, em_p, tolerance = 1e-5)
+  }
+  p_none <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = "none")$p
+  p_holm <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = "holm")$p
+  expect_false(isTRUE(all.equal(p_none, p_holm)))
+})
+
+
+test_that("p_adjust tukey and sidak, slopes compared pairwise within groups", {
+  m3 <- lm(y ~ x * g * h, data = d_adjust)
+  em <- emmeans::emtrends(m3, ~ g | h, var = "x")
+  for (p_adjust in c("tukey", "sidak")) {
+    expect_no_message(
+      {
+        out <- estimate_slopes(
+          m3,
+          "x",
+          by = c("g", "h"),
+          comparison = ~ pairwise | h,
+          p_adjust = p_adjust
+        )
+      },
+      message = "Sidak"
+    )
+    expect_equal(out$p, emmeans_p(out, pairs(em, adjust = p_adjust), by = "h"), tolerance = 1e-5)
+  }
+})
+
+
+test_that("p_adjust, families of grouped contrasts and of contrasts of predictions", {
+  # contrasts of slopes with a group term, Tukey counts the slopes in each group
+  m3 <- lm(y ~ x * g * h, data = d_adjust)
+  out <- estimate_contrasts(
+    m3,
+    contrast = "x",
+    by = c("g", "h"),
+    comparison = ~ pairwise | h,
+    p_adjust = "tukey"
+  )
+  em <- emmeans::emtrends(m3, "g", by = "h", var = "x")
+  expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "tukey"), by = "h"), tolerance = 1e-5)
+
+  # contrasts of `estimate_relation()` output, Sidak counts within `by` groups
+  out <- estimate_contrasts(
+    estimate_relation(m1, by = c("g", "h")),
+    contrast = "g",
+    by = "h",
+    p_adjust = "sidak"
+  )
+  em <- suppressMessages(emmeans::emmeans(m1, "g", by = "h"))
+  expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "sidak"), by = "h"), tolerance = 1e-5)
+
+  # missing values of a `by` variable form their own family
+  expect_identical(
+    modelbased:::.p_adjust_family_size(data.frame(h = c("a", NA, "a", "b")), by = "h"),
+    c(2L, 1L, 2L, 1L)
+  )
+})
+
+
+test_that("p_adjust, no Sidak message when nothing is adjusted, and method names", {
+  # one contrast or one slope: Sidak does not change the p-value
+  expect_no_message(
+    estimate_contrasts(m1, "g", comparison = "b1 = b2", p_adjust = "tukey"),
+    message = "Sidak"
+  )
+  expect_no_message(
+    estimate_slopes(lm(mpg ~ wt, data = mtcars), "wt", hypothesis = 1, p_adjust = "tukey"),
+    message = "Sidak"
+  )
+  # a matrix hypothesis is not a test against a value
+  expect_message(
+    estimate_slopes(
+      m2,
+      "x",
+      by = "g",
+      hypothesis = matrix(c(-1, 1, 0, 0, 0, -1, 1, 0), 4),
+      p_adjust = "tukey"
+    ),
+    regexp = "not compared pairwise"
+  )
+
+  # emmeans backend: method names in any case, and the `p_adjust` attribute
+  mt <- mtcars
+  mt$cyl <- factor(mt$cyl)
+  m4 <- lm(mpg ~ wt * cyl, data = mt)
+  em <- emmeans::emtrends(m4, "cyl", var = "wt")
+  out <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = "bh")
+  expect_equal(
+    out$p,
+    as.data.frame(summary(emmeans::test(em, adjust = "BH")))$p.value,
+    tolerance = 1e-5
+  )
+  expect_identical(attributes(out)$p_adjust, "BH")
+  out <- suppressMessages(estimate_slopes(
+    m4,
+    "wt",
+    by = "cyl",
+    backend = "emmeans",
+    p_adjust = "tukey"
+  ))
+  expect_identical(attributes(out)$p_adjust, "sidak")
+})