Getting adjusted average treatment effect when factors specified as random effects

I’m trying to calculate a covariate-adjusted average treatment effect (ATE) for an experiment. I think I know how to do so using brms when variables are treated as fixed effects, but I’m unsure if I’m doing it correctly when I want to use partial pooling and specify the factors as random effects. Here’s some example code:

First install the necessary packages and make some fake data:

# List required packages
list.of.packages <- c("brms", "tidyverse", "randomizr", "tidybayes")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]

# Install packages if missing
if(length(new.packages)) install.packages(new.packages)

# Set seed so results reproducible
set.seed(26)

# Generate data
dat <- tibble(
  y = rnorm(n = 100, mean = 0, sd = 1),
  treatment = randomizr::simple_ra(N = 100, 
                           prob_each = c(.22, .22, .22, .22, .12), 
                           conditions = c("treat1", "treat2","treat3","treat4","control")),
  gender = randomizr::simple_ra(N = 100, 
                                prob_each = c(.48, .48, .04), 
                                conditions = c("male", "female","other")))

Now when I typically compute an adjusted ATE (i.e., what is the effect of the treatment once I control for pre-existing differences on the categorical covariate), I typically use the marginal effects call from brms like this (the plot for treatment alone I think shows the covariate adjusted ATE):

# Fixed effects
fixed <- brm(y ~ treatment + gender + treatment*gender, data = dat)

# Marginal effects
marginal_effects(fixed)

Now I’d like to get the ATE while adjusting for gender and using partial pooling. Here’s my stab at it:

# Unadjusted
random <- brm(y ~ (1 | treatment) + (1 | gender) + (1| treatment:gender), data = dat)

# Adjusted means
tidybayes::get_variables(random)
random %>%
  tidybayes::spread_draws(b_Intercept, r_treatment[condition, ]) %>%
  tidybayes::median_qi(adjusted_ate = b_Intercept + r_treatment)

So my question is whether this is the best (or correct way) to get a covariate adjusted ATE when specifying factors as random (instead of fixed) effects? I ask because I can’t get the ATE or covariate adjusted ATE using the marginal_effects call when a factor is specified as a random effect.

2 Likes

Hey Andrew!

This seems correct. You can get the same result using

coef(random, robust = TRUE)$treatment

Cheers!
Max

1 Like

I am not familiar with the term ATE, but one thing that I find weird is that you seem to ignore the interaction (1| treatment:gender) term in your estimate of the effect. This might be OK for some use cases, but I am not 100% sure it is always appropriate. Using the model you’ve given, I can look at the mean of the interaction coefficients:

inter_coef <- coef(random, summary = FALSE)$`treatment:gender`
hist(rowMeans(inter_coef[,,1]))

And I see that for many samples the average is quite far from zero, so the interaction effects may take a bit of the weight from the main effect… (i.e. they all may be positive for a given treatment so taking only the main effect of treatment would underestimate the total effect you model). This might cancel out on average, but this is AFAIK not guaranteed to happen.

I usually find it easier to make inferences in such cases via the predictions of the model, as this makes it explicit which predictors do I use and in what way.

I’ve written about those issues recently for other people at Inferences marginal of random effects in multi-level models and Multivariate model interpreting residual correlation and group correlation so check if that helps (but feel free to ask for clarifications if that’s not clear).

Best of luck with your model!