Pairwise comparisons of slopes

[andieich]

Is it possible to do pairwise comparisons of slopes in an ANCOVA?

Consider this model

library(palmerpenguins)
library(modelbased)
m <- lm(flipper_length_mm ~ body_mass_g * species, data = penguins)

I want to recreate the following, which is easy to do in emmeans:

library(emmeans)
emtrends(m, ~ species, var = "body_mass_g") %>% 
    pairs()

It gives

 contrast           estimate      SE  df t.ratio p.value
 Adelie - Chinstrap -0.00529 0.00196 327  -2.704  0.0197
 Adelie - Gentoo    -0.00273 0.00138 327  -1.978  0.1192
 Chinstrap - Gentoo  0.00256 0.00196 327   1.305  0.3933

Running what I thought was the equivalent in modelbased:

estimate_contrasts(
  m,
  contrast = "species",
  trend = "body_mass_g",
  p_adjust = "tukey"
)

Gives

Marginal Contrasts Analysis

Level1    | Level2    | Difference |   SE |         95% CI | t(327) |      p
----------------------------------------------------------------------------
Chinstrap | Adelie    |       8.05 | 1.23 | [ 5.64, 10.46] |   6.57 | < .001
Gentoo    | Adelie    |      15.55 | 1.20 | [13.20, 17.90] |  13.01 | < .001
Gentoo    | Chinstrap |       7.50 | 1.44 | [ 4.67, 10.33] |   5.21 | < .001

Variable predicted: flipper_length_mm
Predictors contrasted: species
Predictors averaged: body_mass_g (4.2e+03)
p-value adjustment method: Tukey

It seems to be comparisons for flipper_length_mm along the mean of body_mass_g.

Is there a way to do pairwise comparisons of slopes? I would like to use the modelbased package for an introductory stats class and I like its functions and output.

[strengejacke]

This is how you can do it for contrasting slopes. The idea is: you have by to stratify, contrast to define for which variable you would like to calculate the (pairwise) comparisons or contrasts, and trend(or slope in a future version) to define the variable for which you would like to calculate the marginal effects / slope.

library(palmerpenguins)
#> 
#> Attaching package: 'palmerpenguins'
#> The following objects are masked from 'package:datasets':
#> 
#>     penguins, penguins_raw
library(modelbased)
library(emmeans)
#> Welcome to emmeans.
#> Caution: You lose important information if you filter this package's results.
#> See '? untidy'

m <- lm(flipper_length_mm ~ body_mass_g * species, data = penguins)

emtrends(m, ~ species, var = "body_mass_g") |> pairs()
#>  contrast           estimate      SE  df t.ratio p.value
#>  Adelie - Chinstrap -0.00523 0.00195 336  -2.683  0.0209
#>  Adelie - Gentoo    -0.00236 0.00135 336  -1.746  0.1896
#>  Chinstrap - Gentoo  0.00287 0.00195 336   1.468  0.3078
#> 
#> P value adjustment: tukey method for comparing a family of 3 estimates

estimate_contrasts(
  m,
  contrast = "body_mass_g",
  by = "species",
  p_adjust = "tukey"
)
#> Tukey adjustment requires at least 2 groups. P-values were not adjusted.
#> Marginal Contrasts Analysis
#> 
#> Level1    | Level2    | Difference |       SE |        95% CI | t(336) |     p
#> ------------------------------------------------------------------------------
#> Chinstrap | Adelie    |   5.23e-03 | 1.95e-03 | [ 0.00, 0.01] |   2.68 | 0.008
#> Gentoo    | Adelie    |   2.36e-03 | 1.35e-03 | [ 0.00, 0.01] |   1.75 | 0.082
#> Gentoo    | Chinstrap |  -2.87e-03 | 1.96e-03 | [-0.01, 0.00] |  -1.47 | 0.144
#> 
#> Variable predicted: flipper_length_mm
#> Predictors contrasted: body_mass_g
#> Predictors averaged: body_mass_g (4.2e+03)
#> p-value adjustment method: Tukey

^{Created on 2026-02-17 with reprex v2.1.1}

However, I realized that the Tukey-adjustment is not yet working for slopes. Will try to fix is.

[andieich]

Great, thanks a lot for your quick explanations and the fix!

[strengejacke]

Should work now when #600 is merged:

data(penguins)
m <- lm(flipper_len ~ body_mass * species, data = penguins)

pairs(emmeans::emtrends(m, ~species, var = "body_mass"))
#>  contrast           estimate      SE  df t.ratio p.value
#>  Adelie - Chinstrap -0.00523 0.00195 336  -2.683  0.0209
#>  Adelie - Gentoo    -0.00236 0.00135 336  -1.746  0.1896
#>  Chinstrap - Gentoo  0.00287 0.00195 336   1.468  0.3078
#> 
#> P value adjustment: tukey method for comparing a family of 3 estimates

modelbased::estimate_contrasts(
  m,
  contrast = "body_mass",
  by = "species",
  p_adjust = "tukey"
)
#> Marginal Contrasts Analysis
#> 
#> Level1    | Level2    | Difference |       SE |        95% CI | t(336) |     p
#> ------------------------------------------------------------------------------
#> Chinstrap | Adelie    |   5.23e-03 | 1.95e-03 | [ 0.00, 0.01] |   2.68 | 0.021
#> Gentoo    | Adelie    |   2.36e-03 | 1.35e-03 | [ 0.00, 0.01] |   1.75 | 0.190
#> Gentoo    | Chinstrap |  -2.87e-03 | 1.96e-03 | [-0.01, 0.00] |  -1.47 | 0.309
#> 
#> Variable predicted: flipper_len
#> Predictors contrasted: body_mass
#> Predictors averaged: body_mass (4.2e+03)
#> p-value adjustment method: Tukey

^{Created on 2026-02-17 with reprex v2.1.1}

[strengejacke]

My comment wasn't complete, so here another try:

• in all estimate_*() functions: you have by to stratify
• in estimate_contrast(): contrast to define for which variable you would like to calculate the (pairwise) comparisons or contrasts
• in estimate_slopes(): trend (or slope in a future version, see Renaming argument trend into slope for estimate_slopes()? #596) to define the variable for which you would like to calculate the marginal effects / slope.

[strengejacke]

btw, see also vignettes under "Case Studies" -> "Contrasts": https://easystats.github.io/modelbased/articles/overview_of_vignettes.html