Pairwise comparisons of slopes (#600)

* Pairwise comparisons of slopes
Fixes #599

* formatting, use "datasets"

Author: strengejacke

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.13.1.8
+Version: 0.13.1.9
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",

R/estimate_slopes.R

@@ -108,7 +108,7 @@
 #' @examplesIf all(insight::check_if_installed(c("marginaleffects", "emmeans"), quietly = TRUE)) && getRversion() >= "4.5.0"
 #' \dontrun{
 #' # marginal effects with different `estimate` options
-#' data(penguins)
+#' data(penguins, package = "datasets")
 #' penguins$long_bill <- factor(datawizard::categorize(penguins$bill_len), labels = c("short", "long"))
 #' m <- glm(long_bill ~ sex + species + island * bill_dep, data = penguins, family = "binomial")
 #'
@@ -121,28 +121,27 @@
 #' marginaleffects::avg_slopes(m, variables = "bill_dep", by = "island")
 #' }
 #' @export
-estimate_slopes <- function(model,
-                            trend = NULL,
-                            by = NULL,
-                            predict = NULL,
-                            ci = 0.95,
-                            estimate = NULL,
-                            transform = NULL,
-                            p_adjust = "none",
-                            keep_iterations = FALSE,
-                            backend = NULL,
-                            verbose = TRUE,
-                            ...) {
+estimate_slopes <- function(
+  model,
+  trend = NULL,
+  by = NULL,
+  predict = NULL,
+  ci = 0.95,
+  estimate = NULL,
+  transform = NULL,
+  p_adjust = "none",
+  keep_iterations = FALSE,
+  backend = NULL,
+  verbose = TRUE,
+  ...
+) {
   # Process argument ---------------------------------------------------------
   if (is.null(backend)) {
     backend <- getOption("modelbased_backend", "marginaleffects")
   }
 
   # validate input
-  estimate <- insight::validate_argument(
-    estimate,
-    c("typical", "specific", "average")
-  )
+  estimate <- insight::validate_argument(estimate, c("typical", "specific", "average"))
 
   if (backend == "emmeans") {
     # Emmeans ----------------------------------------------------------------

R/p_adjust.R

@@ -8,9 +8,17 @@
 
   # extract information
   datagrid <- attributes(params)$datagrid
-  focal <- .safe(insight::trim_ws(gsub("=.*", "\\1", attributes(params)$contrast)))
+
+  # when contrasting slopes, we need "by" as focal term
+  if (inherits(params, "estimate_slopes")) {
+    focal <- .safe(insight::trim_ws(gsub("=.*", "\\1", attributes(params)$by)))
+  } else {
+    focal <- .safe(insight::trim_ws(gsub("=.*", "\\1", attributes(params)$contrast)))
+  }
+
   # extract degrees of freedom
   dof <- .safe(params$df[1])
+
   if (is.null(dof)) {
     dof <- insight::get_df(model, type = "wald", verbose = FALSE)
   }

man/estimate_slopes.Rd

@@ -703,7 +703,7 @@ test_that("estimate_contrasts - p.adjust", {
   skip_if(getRversion() < "4.5.0")
   skip_if_not_installed("emmeans")
 
-  data(penguins)
+  data(penguins, package = "datasets")
   m_spec <- lm(body_mass ~ species, data = penguins)
   out1 <- estimate_contrasts(m_spec, p_adjust = "tukey")
   out2 <- estimate_contrasts(m_spec, p_adjust = "tukey", backend = "emmeans")
@@ -1709,7 +1709,7 @@ test_that("estimate_contrast, informative error when `by` and `contrast` are the
 
 test_that("estimate_contrast, works with aov (when no statistic is extracted)", {
   skip_if(getRversion() < "4.5.0")
-  data(penguins)
+  data(penguins, package = "datasets")
   fit <- aov(formula = body_mass ~ species, data = penguins)
 
   out1 <- marginaleffects::avg_predictions(fit, by = "species", hypothesis = ~pairwise)
@@ -1736,7 +1736,7 @@ test_that("estimate_contrast, works with aov (when no statistic is extracted)",
 test_that("estimate_contrast, slopes with different estimate options", {
   skip_if(getRversion() < "4.5.0")
   skip_if_not_installed("datawizard")
-  data(penguins)
+  data(penguins, package = "datasets")
   penguins$long_bill <- factor(
     datawizard::categorize(penguins$bill_len),
     labels = c("short", "long")
@@ -1766,3 +1766,26 @@ test_that("estimate_contrast, comparison-options as strings", {
   out <- estimate_contrasts(mod2, contrast = "cyl_helmert", comparison = "poly")
   expect_equal(out$Difference, c(-8.17673, 0.92996), tolerance = 1e-4)
 })
+
+
+test_that("estimate_contrast, p-adjust tukey works for contrasting slopes", {
+  skip_if(getRversion() < "4.5.0")
+  skip_if_not_installed("emmeans")
+
+  data(penguins, package = "datasets")
+  m <- lm(flipper_len ~ body_mass * species, data = penguins)
+
+  out1 <- as.data.frame(pairs(emmeans::emtrends(m, ~species, var = "body_mass")))
+  out2 <- estimate_contrasts(
+    m,
+    contrast = "body_mass",
+    by = "species",
+    p_adjust = "tukey"
+  )
+
+  # Note: p-values from emmeans::emtrends() + pairs() and estimate_contrasts()
+  # can differ slightly due to different underlying calculation/adjustment
+  # methods, especially with Tukey p-adjustment. A tolerance of 1e-2 is used
+  # here to avoid fragile tests while still ensuring close agreement.
+  expect_equal(out1$p.value, out2$p, tolerance = 1e-2)
+})