From 5894e9ac332d978059c91639b5b18d17b5c1c22b Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 10:59:37 -0500
Subject: [PATCH 01/10] Use Sidak for p_adjust = "tukey" on contrasts that are
 not pairwise, and count the comparisons in each by group for Sidak

---
 R/p_adjust.R                   |  61 ++++++++++++++-
 tests/testthat/test-p_adjust.R | 133 +++++++++++++++++++++++++++++++++
 2 files changed, 192 insertions(+), 2 deletions(-)
 create mode 100644 tests/testthat/test-p_adjust.R

diff --git a/R/p_adjust.R b/R/p_adjust.R
index a8f409f72..cff64ed24 100644
--- a/R/p_adjust.R
+++ b/R/p_adjust.R
@@ -80,6 +80,26 @@
     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
+  if (
+    p_adjust == "tukey" &&
+      !is_unformatted &&
+      !is.null(comparison) &&
+      !.is_pairwise_comparison(comparison) &&
+      "p" %in% colnames(params)
+  ) {
+    if (verbose) {
+      insight::format_alert(
+        "The Tukey adjustment is only appropriate for pairwise comparisons, so p-values were adjusted with the Sidak method instead."
+      )
+    }
+    p_adjust <- "sidak"
+    attr(params, "p_adjust") <- "sidak"
+  }
+
   # esarey is specifically for johnson-neyman intervals
   if (p_adjust == "esarey") {
     return(.p_adjust_esarey(params))
@@ -150,14 +170,51 @@
       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 <- nrow(params)
+    } else {
+      family_size <- .p_adjust_family_size(params)
+    }
+    params[[p_column]] <- 1 - (1 - params[[p_column]])^family_size
   }
 
   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")
+}
+
+
+# number of contrasts in the family of each row: the rows with the same values
+# of the `by` variables
+.p_adjust_family_size <- function(params) {
+  by <- .safe(insight::trim_ws(gsub("=.*", "", attributes(params)$by)))
+  by <- intersect(by, colnames(params))
+  if (!length(by)) {
+    return(rep(nrow(params), nrow(params)))
+  }
+  groups <- interaction(params[by], drop = TRUE)
+  as.vector(table(groups)[as.character(groups)])
+}
+
+
 .p_adjust_supt <- function(params, model) {
   insight::check_if_installed("mvtnorm")
   # get correlation matrix, based on the covariance matrix
diff --git a/tests/testthat/test-p_adjust.R b/tests/testthat/test-p_adjust.R
new file mode 100644
index 000000000..3cd119822
--- /dev/null
+++ b/tests/testthat/test-p_adjust.R
@@ -0,0 +1,133 @@
+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 = "|"), "")
+  key_out <- vapply(
+    seq_len(nrow(out)),
+    function(i) paste(sort(c(as.character(out$Level1[i]), as.character(out$Level2[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)
+  }
+})

From a6762a8372bbd93aa6e05d85d67b8e508eb18dfd Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 11:13:50 -0500
Subject: [PATCH 02/10] Count the compared slopes for p_adjust = "tukey", and
 use Sidak for slopes tested against a value or not compared pairwise

---
 R/p_adjust.R                   | 71 ++++++++++++++++++++---------
 tests/testthat/test-p_adjust.R | 82 +++++++++++++++++++++++++++++++---
 2 files changed, 128 insertions(+), 25 deletions(-)

diff --git a/R/p_adjust.R b/R/p_adjust.R
index cff64ed24..f3cf3fe69 100644
--- a/R/p_adjust.R
+++ b/R/p_adjust.R
@@ -56,28 +56,56 @@
   }
 
   # 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(...)
-  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"]])
-  ) {
-    if (verbose) {
-      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."
+  n_slopes <- NULL
+  if (p_adjust == "tukey" && is_unformatted && "p.value" %in% colnames(params)) {
+    dots_names <- names(dots)
+    if (is.null(dots_names)) {
+      dots_names <- character(0)
+    }
+    hypothesis_arg <- dots_names[nzchar(dots_names) & startsWith("hypothesis", dots_names)]
+    slope_comparison <- c(
+      dots["comparison"],
+      dots[hypothesis_arg],
+      dots["post_process"]
+    )
+    slope_comparison <- Filter(Negate(is.null), slope_comparison)
+    reason <- NULL
+    if (!length(slope_comparison)) {
+      if (nrow(params) > 1) {
+        reason <- "Slopes are tested against 0"
+      }
+    } else if (is.numeric(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
+      n_slopes <- (1 + sqrt(1 + 8 * nrow(params))) / 2
+      if (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)) {
+      if (verbose) {
+        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"
     }
-    p_adjust <- "sidak"
-    attr(params, "p_adjust") <- "sidak"
   }
 
   # contrasts (formatted output): the Tukey adjustment is only appropriate for
@@ -112,7 +140,10 @@
 
   # 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
diff --git a/tests/testthat/test-p_adjust.R b/tests/testthat/test-p_adjust.R
index 3cd119822..173ca6fde 100644
--- a/tests/testthat/test-p_adjust.R
+++ b/tests/testthat/test-p_adjust.R
@@ -24,11 +24,17 @@ 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 = "|"), "")
-  key_out <- vapply(
-    seq_len(nrow(out)),
-    function(i) paste(sort(c(as.character(out$Level1[i]), as.character(out$Level2[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]])
@@ -131,3 +137,69 @@ test_that("p_adjust sidak, contrasts, family size, marginaleffects", {
     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)
+})

From ab2673235971b55b5124015d9bc53a3bceba6ec6 Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 11:18:34 -0500
Subject: [PATCH 03/10] Apply p_adjust to slopes from the emmeans backend

---
 R/estimate_slopes.R            | 22 ++++++++++++++++++++++
 tests/testthat/test-p_adjust.R | 16 ++++++++++++++++
 2 files changed, 38 insertions(+)

diff --git a/R/estimate_slopes.R b/R/estimate_slopes.R
index 12b1e7646..3a6870ab8 100644
--- a/R/estimate_slopes.R
+++ b/R/estimate_slopes.R
@@ -190,6 +190,28 @@ estimate_slopes <- function(
       verbose = verbose,
       ...
     )
+    # p-value adjustment, with the methods emmeans provides
+    emmeans_adjust <- c(
+      stats::p.adjust.methods,
+      "tukey",
+      "sidak",
+      "scheffe",
+      "mvt",
+      "dunnettx"
+    )
+    if (!is.null(p_adjust) && tolower(p_adjust) %in% tolower(emmeans_adjust)) {
+      # slopes are tested against 0, so we use Sidak instead of Tukey, as
+      # `emmeans::test()` does
+      if (tolower(p_adjust) == "tukey") {
+        if (verbose) {
+          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."
+          )
+        }
+        p_adjust <- "sidak"
+      }
+      estimated <- stats::update(estimated, adjust = p_adjust)
+    }
     trends <- .format_emmeans_slopes(model, estimated, ci, ...)
   } else {
     # Marginaleffects --------------------------------------------------------
diff --git a/tests/testthat/test-p_adjust.R b/tests/testthat/test-p_adjust.R
index 173ca6fde..a291b015f 100644
--- a/tests/testthat/test-p_adjust.R
+++ b/tests/testthat/test-p_adjust.R
@@ -203,3 +203,19 @@ test_that("p_adjust tukey, slopes tested against a value, marginaleffects", {
   )
   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)))
+})

From 3f626de8914ba5a2e031ae3b7cafb7a492708af8 Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 11:24:21 -0500
Subject: [PATCH 04/10] Document the Tukey and Sidak rules for p_adjust, add
 NEWS, bump version

---
 NEWS.md                   | 29 +++++++++++++++++++++++++++--
 R/estimate_contrasts.R    |  8 +++++++-
 R/estimate_slopes.R       | 35 ++++++++++++++++-------------------
 man/estimate_contrasts.Rd |  8 +++++++-
 man/estimate_slopes.Rd    | 23 +++++++++++++----------
 man/get_emmeans.Rd        | 23 +++++++++++++----------
 6 files changed, 83 insertions(+), 43 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index d497a6dca..b7b1b9126 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -155,8 +155,33 @@
   `get_marginaltrends()` now use the Sidak adjustment and say so in a message.
   The Tukey adjustment is only appropriate for pairwise comparisons, but these
   slopes are each tested against 0. `emmeans::test()` makes the same
-  substitution. Slopes that are compared, for example with the `comparison`
-  or `post_process` argument, keep the Tukey adjustment.
+  substitution. Slopes that are compared pairwise, for example with
+  `comparison = ~pairwise` or `post_process = ~pairwise`, keep the Tukey
+  adjustment.
+
+* `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.
+  * The Tukey adjustment of slopes compared pairwise with `estimate_slopes()`
+    now counts the slopes, not the comparisons. Before, 4 slopes were counted
+    as 6 groups, which made the p-values too large.
+  * 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 partly matched `hypothesis` argument, for example `hypo = ~pairwise`, is
+    now treated as `hypothesis`.
+
+* `estimate_slopes()` with `backend = "emmeans"` now applies `p_adjust`. 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.
 
 * `estimate_means()` with `backend = "emmeans"` now shows the response level
   names for models from `ordinal::clmm()`. Before, the `Response` column showed
diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index 38c9bb143..243affa40 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -26,7 +26,13 @@
 #' 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,
+#' and for contrasts of slopes, 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 3a6870ab8..85edb26d6 100644
--- a/R/estimate_slopes.R
+++ b/R/estimate_slopes.R
@@ -26,19 +26,22 @@
 #' `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 all tests (all rows) of the output. This
+#' also holds for `backend = "emmeans"`, which applies `p_adjust` to the
+#' slopes.
 #' @inheritParams estimate_means
 #' @inheritParams parameters::model_parameters.default
 #'
@@ -190,15 +193,9 @@ estimate_slopes <- function(
       verbose = verbose,
       ...
     )
-    # p-value adjustment, with the methods emmeans provides
-    emmeans_adjust <- c(
-      stats::p.adjust.methods,
-      "tukey",
-      "sidak",
-      "scheffe",
-      "mvt",
-      "dunnettx"
-    )
+    # p-value adjustment, done by `parameters::model_parameters()` for the
+    # "adjust" setting of the emmeans object
+    emmeans_adjust <- c(stats::p.adjust.methods, "tukey", "sidak")
     if (!is.null(p_adjust) && tolower(p_adjust) %in% tolower(emmeans_adjust)) {
       # slopes are tested against 0, so we use Sidak instead of Tukey, as
       # `emmeans::test()` does
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index bf67bb755..de718c29d 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -280,7 +280,13 @@ 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,
+and for contrasts of slopes, 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..a7c8eca86 100644
--- a/man/estimate_slopes.Rd
+++ b/man/estimate_slopes.Rd
@@ -172,19 +172,22 @@ 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 all tests (all rows) of the output. This
+also holds for \code{backend = "emmeans"}, which applies \code{p_adjust} to the
+slopes.}
 
 \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..7d68f422f 100644
--- a/man/get_emmeans.Rd
+++ b/man/get_emmeans.Rd
@@ -362,19 +362,22 @@ 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 all tests (all rows) of the output. This
+also holds for \code{backend = "emmeans"}, which applies \code{p_adjust} to the
+slopes.}
 
 \item{iterations}{This argument has two distinct effects, depending on the
 context. When \code{effectsize = "boot"}, it defines the number of bootstrap

From 2fecf9307f787ca40a92385726cc9e35a3f538f1 Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 11:31:15 -0500
Subject: [PATCH 05/10] Count Tukey and Sidak families in each group for slopes
 compared with a group term, and narrow the p_adjust docs

---
 NEWS.md                        | 13 +++++----
 R/estimate_contrasts.R         |  6 ++--
 R/estimate_slopes.R            |  8 +++--
 R/p_adjust.R                   | 53 +++++++++++++++++++++++++++-------
 man/estimate_contrasts.Rd      |  6 ++--
 man/estimate_slopes.Rd         |  8 +++--
 man/get_emmeans.Rd             |  8 +++--
 tests/testthat/test-p_adjust.R | 21 ++++++++++++++
 8 files changed, 95 insertions(+), 28 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index b7b1b9126..9e09b8aef 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -170,7 +170,9 @@
     is replaced by `"sidak"` with a message.
   * The Tukey adjustment of slopes compared pairwise with `estimate_slopes()`
     now counts the slopes, not the comparisons. Before, 4 slopes were counted
-    as 6 groups, which made the p-values too large.
+    as 6 groups, 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
@@ -178,10 +180,11 @@
   * A partly matched `hypothesis` argument, for example `hypo = ~pairwise`, is
     now treated as `hypothesis`.
 
-* `estimate_slopes()` with `backend = "emmeans"` now applies `p_adjust`. 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.
+* `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.
 
 * `estimate_means()` with `backend = "emmeans"` now shows the response level
   names for models from `ordinal::clmm()`. Before, the `Response` column showed
diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R
index 243affa40..836c2fbda 100644
--- a/R/estimate_contrasts.R
+++ b/R/estimate_contrasts.R
@@ -31,8 +31,10 @@
 #' `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,
-#' and for contrasts of slopes, all comparisons.
+#' 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 85edb26d6..ff6f2dbe6 100644
--- a/R/estimate_slopes.R
+++ b/R/estimate_slopes.R
@@ -39,9 +39,11 @@
 #' 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 all tests (all rows) of the output. This
-#' also holds for `backend = "emmeans"`, which applies `p_adjust` to the
-#' slopes.
+#' does. The Tukey and Sidak adjustments count all tests (all rows) of the
+#' output, or, for comparisons with a group term (for example
+#' `comparison = ~pairwise | group`), the tests in each group. 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
 #'
diff --git a/R/p_adjust.R b/R/p_adjust.R
index f3cf3fe69..0b23413a3 100644
--- a/R/p_adjust.R
+++ b/R/p_adjust.R
@@ -63,7 +63,12 @@
   # p-values, the code below says that no adjustment was made.
   dots <- list(...)
   n_slopes <- NULL
-  if (p_adjust == "tukey" && is_unformatted && "p.value" %in% colnames(params)) {
+  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)
@@ -75,6 +80,14 @@
       dots["post_process"]
     )
     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) {
@@ -86,9 +99,15 @@
         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
-      n_slopes <- (1 + sqrt(1 + 8 * nrow(params))) / 2
-      if (n_slopes != round(n_slopes)) {
+      # 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"
       }
@@ -176,7 +195,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."
@@ -205,7 +224,7 @@
     # and for contrasts, the rows within each group of the `by` variables
     # (like emmeans).
     if (is_unformatted) {
-      family_size <- nrow(params)
+      family_size <- slope_family
     } else {
       family_size <- .p_adjust_family_size(params)
     }
@@ -233,10 +252,24 @@
 }
 
 
-# number of contrasts in the family of each row: the rows with the same values
-# of the `by` variables
-.p_adjust_family_size <- function(params) {
-  by <- .safe(insight::trim_ws(gsub("=.*", "", attributes(params)$by)))
+# 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)))
diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd
index de718c29d..0ed291730 100644
--- a/man/estimate_contrasts.Rd
+++ b/man/estimate_contrasts.Rd
@@ -285,8 +285,10 @@ 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,
-and for contrasts of slopes, all comparisons.}
+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 a7c8eca86..ef441d041 100644
--- a/man/estimate_slopes.Rd
+++ b/man/estimate_slopes.Rd
@@ -185,9 +185,11 @@ example with \code{comparison = ~pairwise} or \code{post_process = ~pairwise}. F
 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 all tests (all rows) of the output. This
-also holds for \code{backend = "emmeans"}, which applies \code{p_adjust} to the
-slopes.}
+does. The Tukey and Sidak adjustments count all tests (all rows) of the
+output, or, for comparisons with a group term (for example
+\code{comparison = ~pairwise | group}), the tests in each group. 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 7d68f422f..57058a12b 100644
--- a/man/get_emmeans.Rd
+++ b/man/get_emmeans.Rd
@@ -375,9 +375,11 @@ example with \code{comparison = ~pairwise} or \code{post_process = ~pairwise}. F
 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 all tests (all rows) of the output. This
-also holds for \code{backend = "emmeans"}, which applies \code{p_adjust} to the
-slopes.}
+does. The Tukey and Sidak adjustments count all tests (all rows) of the
+output, or, for comparisons with a group term (for example
+\code{comparison = ~pairwise | group}), the tests in each group. 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
index a291b015f..7ebb9f43a 100644
--- a/tests/testthat/test-p_adjust.R
+++ b/tests/testthat/test-p_adjust.R
@@ -219,3 +219,24 @@ test_that("p_adjust, slopes, emmeans backend", {
   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)
+  }
+})

From fde4121952076f76ee416502ac1fcdf7b9426335 Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 11:31:59 -0500
Subject: [PATCH 06/10] Say what the Tukey and Sidak adjustments count for
 slopes

---
 R/estimate_slopes.R    | 6 +++---
 man/estimate_slopes.Rd | 6 +++---
 man/get_emmeans.Rd     | 6 +++---
 3 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/R/estimate_slopes.R b/R/estimate_slopes.R
index ff6f2dbe6..fc4caeb11 100644
--- a/R/estimate_slopes.R
+++ b/R/estimate_slopes.R
@@ -39,9 +39,9 @@
 #' 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 Tukey and Sidak adjustments count all tests (all rows) of the
-#' output, or, for comparisons with a group term (for example
-#' `comparison = ~pairwise | group`), the tests in each group. With
+#' 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. 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
diff --git a/man/estimate_slopes.Rd b/man/estimate_slopes.Rd
index ef441d041..1c80ca769 100644
--- a/man/estimate_slopes.Rd
+++ b/man/estimate_slopes.Rd
@@ -185,9 +185,9 @@ example with \code{comparison = ~pairwise} or \code{post_process = ~pairwise}. F
 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 Tukey and Sidak adjustments count all tests (all rows) of the
-output, or, for comparisons with a group term (for example
-\code{comparison = ~pairwise | group}), the tests in each group. With
+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. 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.}
 
diff --git a/man/get_emmeans.Rd b/man/get_emmeans.Rd
index 57058a12b..6c9f7b7e3 100644
--- a/man/get_emmeans.Rd
+++ b/man/get_emmeans.Rd
@@ -375,9 +375,9 @@ example with \code{comparison = ~pairwise} or \code{post_process = ~pairwise}. F
 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 Tukey and Sidak adjustments count all tests (all rows) of the
-output, or, for comparisons with a group term (for example
-\code{comparison = ~pairwise | group}), the tests in each group. With
+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. 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.}
 

From 2c9b6a1c39a6fe85e54c6920cd9b8f1f4e877623 Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 11:32:43 -0500
Subject: [PATCH 07/10] Bump version

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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",

From d507fb4512f93d8157c48c433127e04f60905dfb Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 11:52:42 -0500
Subject: [PATCH 08/10] Count p_adjust families for grouped contrasts of slopes
 and contrasts of predictions, accept any case of method names with the
 emmeans backend, and skip the Sidak message when nothing is adjusted

---
 NEWS.md                        |  5 ++-
 R/estimate_contrast_methods.R  |  4 +-
 R/estimate_slopes.R            | 28 +++++++++----
 R/p_adjust.R                   | 44 ++++++++++++++++----
 man/estimate_slopes.Rd         |  7 ++--
 man/get_emmeans.Rd             |  7 ++--
 tests/testthat/test-p_adjust.R | 76 ++++++++++++++++++++++++++++++++++
 7 files changed, 144 insertions(+), 27 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 9e09b8aef..f860b1c04 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -170,13 +170,14 @@
     is replaced by `"sidak"` with a message.
   * The Tukey adjustment of slopes compared pairwise with `estimate_slopes()`
     now counts the slopes, not the comparisons. Before, 4 slopes were counted
-    as 6 groups, which made the p-values too large. With a group term, for
+    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.
+    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`.
 
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_slopes.R b/R/estimate_slopes.R
index fc4caeb11..97c34c602 100644
--- a/R/estimate_slopes.R
+++ b/R/estimate_slopes.R
@@ -41,9 +41,10 @@
 #' `"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. With
-#' `backend = "emmeans"`, only the `stats::p.adjust()` methods, `"tukey"` and
-#' `"sidak"` are applied to the slopes, and other methods are ignored.
+#' 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
 #'
@@ -196,20 +197,29 @@ estimate_slopes <- function(
       ...
     )
     # p-value adjustment, done by `parameters::model_parameters()` for the
-    # "adjust" setting of the emmeans object
+    # "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")
-    if (!is.null(p_adjust) && tolower(p_adjust) %in% tolower(emmeans_adjust)) {
+    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 (tolower(p_adjust) == "tukey") {
-        if (verbose) {
+      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."
           )
         }
-        p_adjust <- "sidak"
+        emmeans_method <- "sidak"
       }
-      estimated <- stats::update(estimated, adjust = p_adjust)
+      estimated <- stats::update(estimated, adjust = emmeans_method)
+      attr(estimated, "p_adjust") <- emmeans_method
     }
     trends <- .format_emmeans_slopes(model, estimated, ci, ...)
   } else {
diff --git a/R/p_adjust.R b/R/p_adjust.R
index 0b23413a3..9c6c045d8 100644
--- a/R/p_adjust.R
+++ b/R/p_adjust.R
@@ -74,10 +74,11 @@
       dots_names <- character(0)
     }
     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],
-      dots["post_process"]
+      dots[hypothesis_arg]
     )
     slope_comparison <- Filter(Negate(is.null), slope_comparison)
     if (length(slope_comparison)) {
@@ -93,7 +94,7 @@
       if (nrow(params) > 1) {
         reason <- "Slopes are tested against 0"
       }
-    } else if (is.numeric(slope_comparison[[1]])) {
+    } 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))
@@ -115,7 +116,9 @@
       reason <- "Slopes are not compared pairwise"
     }
     if (!is.null(reason)) {
-      if (verbose) {
+      # 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,
@@ -131,6 +134,14 @@
   # 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 &&
@@ -138,7 +149,9 @@
       !.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, so p-values were adjusted with the Sidak method instead."
       )
@@ -168,6 +181,17 @@
     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"
@@ -226,9 +250,10 @@
     if (is_unformatted) {
       family_size <- slope_family
     } else {
-      family_size <- .p_adjust_family_size(params)
+      family_size <- contrast_family
     }
-    params[[p_column]] <- 1 - (1 - params[[p_column]])^family_size
+    # 1 - (1 - p)^n, without losing precision for small p-values
+    params[[p_column]] <- -expm1(family_size * log1p(-params[[p_column]]))
   }
 
   params
@@ -274,8 +299,9 @@
   if (!length(by)) {
     return(rep(nrow(params), nrow(params)))
   }
-  groups <- interaction(params[by], drop = TRUE)
-  as.vector(table(groups)[as.character(groups)])
+  # 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)])
 }
 
 
diff --git a/man/estimate_slopes.Rd b/man/estimate_slopes.Rd
index 1c80ca769..18c585d21 100644
--- a/man/estimate_slopes.Rd
+++ b/man/estimate_slopes.Rd
@@ -187,9 +187,10 @@ 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. 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.}
+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 6c9f7b7e3..7fe429137 100644
--- a/man/get_emmeans.Rd
+++ b/man/get_emmeans.Rd
@@ -377,9 +377,10 @@ 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. 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.}
+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
index 7ebb9f43a..07b9010f5 100644
--- a/tests/testthat/test-p_adjust.R
+++ b/tests/testthat/test-p_adjust.R
@@ -240,3 +240,79 @@ test_that("p_adjust tukey and sidak, slopes compared pairwise within groups", {
     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(
+    .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")
+})

From 8565c50cd52878f5fa0b3eeb9fe29ac962849cb5 Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 12:16:24 -0500
Subject: [PATCH 09/10] Explain the partial match of hypothesis, and call the
 internal helper with ::: in tests

---
 R/p_adjust.R                   | 2 ++
 tests/testthat/test-p_adjust.R | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/R/p_adjust.R b/R/p_adjust.R
index 9c6c045d8..ccfc82b65 100644
--- a/R/p_adjust.R
+++ b/R/p_adjust.R
@@ -73,6 +73,8 @@
     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(
diff --git a/tests/testthat/test-p_adjust.R b/tests/testthat/test-p_adjust.R
index 07b9010f5..29b2da0c7 100644
--- a/tests/testthat/test-p_adjust.R
+++ b/tests/testthat/test-p_adjust.R
@@ -267,7 +267,7 @@ test_that("p_adjust, families of grouped contrasts and of contrasts of predictio
 
   # missing values of a `by` variable form their own family
   expect_identical(
-    .p_adjust_family_size(data.frame(h = c("a", NA, "a", "b")), by = "h"),
+    modelbased:::.p_adjust_family_size(data.frame(h = c("a", NA, "a", "b")), by = "h"),
     c(2L, 1L, 2L, 1L)
   )
 })

From d204fb2da656084f330f6feeda17fb6525871e2d Mon Sep 17 00:00:00 2001
From: Jeffrey Girard <me@jmgirard.com>
Date: Sun, 18 Oct 2026 13:14:54 -0500
Subject: [PATCH 10/10] Move the NEWS entries to a new devel section after the
 0.18.0 release

---
 NEWS.md | 67 +++++++++++++++++++++++++++++++--------------------------
 1 file changed, 36 insertions(+), 31 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index f860b1c04..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
@@ -155,37 +189,8 @@
   `get_marginaltrends()` now use the Sidak adjustment and say so in a message.
   The Tukey adjustment is only appropriate for pairwise comparisons, but these
   slopes are each tested against 0. `emmeans::test()` makes the same
-  substitution. Slopes that are compared pairwise, for example with
-  `comparison = ~pairwise` or `post_process = ~pairwise`, keep the Tukey
-  adjustment.
-
-* `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.
-  * 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.
+  substitution. Slopes that are compared, for example with the `comparison`
+  or `post_process` argument, keep the Tukey adjustment.
 
 * `estimate_means()` with `backend = "emmeans"` now shows the response level
   names for models from `ordinal::clmm()`. Before, the `Response` column showed