Calculating average marginal effects in models with random slopes

Interfaces brms
7

@jackbailey Thanks, I took a look and like your approach of doing the calculations and weighting by unique predictor combos to save time. brmsmargins matches your code in cases where they both work.
Although the original uncertainty around whether these are accurate remains for me.

dummy <- mtcars
dummy$Time <- mtcars$am
m3 <- brm(mpg ~ Time + (1 + Time | cyl), data = dummy, cores = 4)
summary(m3)

which gives:

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: mpg ~ Time + (1 + Time | cyl) 
   Data: dummy (Number of observations: 32) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~cyl (Number of levels: 3) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)           5.23      2.84     1.91    12.13 1.00     1783     1967
sd(Time)                3.12      2.40     0.13     9.07 1.00     1387     1677
cor(Intercept,Time)     0.30      0.51    -0.80     0.98 1.00     2404     1988

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    18.86      2.77    13.13    24.55 1.00     1353     1568
Time          2.79      2.45    -2.04     7.88 1.00     1852     1738

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.12      0.43     2.42     4.11 1.00     3105     2706

and the AME for Time

time_ame1 <- bayes_dydx.default(m3, variable = "Time", stepsize = 1)
time_ame1$est %>% quantile(probs = c(.5, .025, .975))

which gives:

       50%       2.5%      97.5% 
 2.5753360 -0.1291799  5.3944124

or

time_ame2 <- brmsmargins(m3,
   add = data.frame(Time = c(0, 1)),
   CI = .95, CIType = "ETI",
   contrasts = matrix(c(-1, 1)))

time_ame2$ContrastSummary

which gives:

          M      Mdn         LL       UL PercentROPE PercentMID   CI CIType
1: 2.600315 2.575336 -0.1291799 5.394412          NA         NA 0.95    ETI
   ROPE  MID      Label
1: <NA> <NA> Contrast_1

To me a notable difference is the population estimate straight from brms is:
for Time: 2.79 [-2.04, 7.88]
compared to the AME for Time: 2.58 [-0.13, 5.39]
with the latter having considerably less uncertainty

On a separate note, if you’re interested in various marginal estimates after brms models, you may find some use in the additional capacity of brmsmargins

dummy <- mtcars
dummy$Time <- mtcars$am
dummy$Tx <- mtcars$vs
m2 <- brm(mpg ~ Time + Time:Tx + (1 + Time + Time:Tx | cyl), data = dummy, cores = 4)
summary(m2)

time_ame2 <- brmsmargins(m2,
            at = expand.grid(Time = 0:1, Tx = 0:1),
            CI = .95, CIType = "ETI",
            contrasts = cbind(
              "Time 0 v 1: Ctrl" = c(-1, +1, 0, 0),
              "Time 0 v 1: Tx"   = c(0, 0, -1, +1),
              "Tx v Ctrl: Time 0" = c(-1, 0, +1, 0),
              "Tx v Ctrl: Time 1" = c(0, -1, 0, +1)))

time_ame2$ContrastSummary

which gives:

          M      Mdn        LL        UL PercentROPE PercentMID   CI CIType
1: 1.494765 1.470411 -1.666221  4.573508          NA         NA 0.95    ETI
2: 5.146171 5.184606 -7.693666 17.631245          NA         NA 0.95    ETI
3: 0.000000 0.000000  0.000000  0.000000          NA         NA 0.95    ETI
4: 3.651406 3.751300 -9.708911 15.813571          NA         NA 0.95    ETI
   ROPE  MID             Label
1: <NA> <NA>  Time 0 v 1: Ctrl
2: <NA> <NA>    Time 0 v 1: Tx
3: <NA> <NA> Tx v Ctrl: Time 0
4: <NA> <NA> Tx v Ctrl: Time 1

these are counterfactuals so what if the population was all given the control treatment at baseline, etc. Also obviously it requires some more manual specification of desired contrasts to be applied. There also is a subset argument, which I use in longitudinal studies as I do not necessarily want to upweight the sample distribution of values in people who completed more time points in a longitudinal study. Anyway see
here if interested.

1 Like