Calculating contrasts of marginal effects with marginaleffects for brms model

Question

I have fitted a logistic model with brms and want to calculate the average marginal effects (AMEs).

library(brms)

model <- brm(formula = outcome ~ var1 + var2 + var3, family = bernoulli(), data = data)

var1 is categorical and has two levels (1, 0), var3 is categorical and has three levels ("a", "b", "c"). var3 is another variable whose value I don't want to set/average when calculating the marginal effects.

Now, I want to calculate the AME/slope of var1 for each level of var2. (I don't want to use average covariates (MEMs).)

library(marginaleffects)

slopes <- avg_slopes(
  model,
  variables = "var1",
  by = "var2"
) %>% posterior_draws()

This gives me:

 Term  Contrast           var2  Estimate   2.5 %    97.5 %
 var1  mean(1) - mean(0)  a     0.0361     -0.1098  0.1735
 var1  mean(1) - mean(0)  b     0.0618     -0.0454  0.1666
 var1  mean(1) - mean(0)  c    -0.0788     -0.1667  0.0177

Now, I need the pairwise contrasts between these contrasts. I.e., the differences in Estimate for a - b, a - c and b - c.

I could do it manually like so:

slopes_a <- slopes %>% filter(var2 == "a")
slopes_b <- slopes %>% filter(var2 == "b")
slopes_c <- slopes %>% filter(var2 == "c")

df <- data.frame(
  `a - b` = slopes_a$draw - slopes_b$draw,
  `a - c` = slopes_a$draw - slopes_c$draw,
  `b - c` = slopes_b$draw - slopes_c$draw
)

Is this a valid way and is there perhaps a better/nicer way, e.g., directly with the marginaleffects package?

Answer 1

Use the hypothesis argument. This is well documented in the man page at ?avg_slopes, and also in a detailed vignette: https://marginaleffects.com/vignettes/hypothesis.html

library(brms)
library(marginaleffects)

dat <- transform(mtcars, cyl = factor(cyl))
mod <- brm(vs ~ mpg * cyl, family = bernoulli(), data = dat)
avg_slopes(mod, 
    variables = "mpg", 
    by = "cyl")

     Term    Contrast cyl  Estimate   2.5 % 97.5 %
      mpg mean(dY/dX)   6 -1.80e-01 -0.2568 0.0216
      mpg mean(dY/dX)   4  2.71e-03 -0.0382 0.0682
      mpg mean(dY/dX)   8  3.66e-09 -0.0469 0.0666

avg_slopes(mod, 
    variables = "mpg", 
    by = "cyl", 
    hypothesis = "pairwise")

      Term Estimate   2.5 % 97.5 %
     6 - 4 -0.18166 -0.2780 0.0172
     6 - 8 -0.17877 -0.2763 0.0257
     4 - 8  0.00294 -0.0679 0.0827
```