Counterfactual prediction in a mediation model without all predictors

Interfaces brms
,
1
  • Operating System: Debian 12
  • brms Version: 2.21.0

I’m new to Bayesian statistics and modeling. I’m using R. McElreath’s “statistical rethinking” 2nd edition book and S. Kurz’s translation into tidyverse and brms. I tried to find an answer but I don’t have the vocabulary in the area to search efficiently.

In chapter 5 practice question H2 we have a model of M → A → D and we are asked to find the effect of changing M.

My model looks like this:

# Loading data
data(WaffleDivorce, package="rethinking")

wd <- WaffleDivorce
wd %<>% mutate(
  d = rethinking::standardize(Divorce),
  m = rethinking::standardize(Marriage),
  a = rethinking::standardize(MedianAgeMarriage)
)

rm(WaffleDivorce)
# Running model
a_model <- bf(a ~ 1 + m) 
d_model <- bf(d ~ 1 + a)
bh5.2 <- brm(data = wd,
             family = gaussian,
             a_model + d_model + set_rescor(FALSE),
             prior = c(
               prior(normal(0, 0.2), class=Intercept, resp=a),
               prior(normal(0, 0.5), class=b, resp=a),
               prior(exponential(1), class=sigma, resp=a),
               prior(normal(0, 0.2), class=Intercept, resp=d),
               prior(normal(0, 0.5), class=b, resp=d),
               prior(exponential(1), class=sigma, resp=d)
             ),
             iter=2000, warmup=1000, cores=4, chains=4,
             seed=5,
             file="fits/bh05.02")

To tackle this issue I want to use a function like predict but it requires all the parameters to be supplied (M & A) and in the book McElreath explains that we first need to calculate A from M (impossible using predict).

I can use a multi-step approach and predict A from M using resp parameter and than predict D from A or use the posterior draws (which is the method in my code).

# Getting counterfactual est
d5.2 <- as_draws_df(bh5.2)
set.seed(5)
h5.2.post <- d5.2 %>%
  expand_grid(m = seq(from=-3, to=3, length.out=30)) %>%
  mutate(a_sim = rnorm(n(), mean=b_a_Intercept + b_a_m * m, sd = sigma_a),
         d_sim = rnorm(n(), mean=b_d_Intercept + b_d_a * a_sim, sd = sigma_d)
         ) %>%
  pivot_longer(cols=ends_with("sim")) %>%
  group_by(m, name) %>%
  summarise(mean=mean(value),
            ll = quantile(value, prob=0.025),
            ul = quantile(value, prob=1-0.025))

Is there a way in brms to get such estimate in a single step?
Something more like this:

predict(bh5.2, newdata=tibble(m=seq(from=-3, to=3, length.out=30), resp="d") # resp for clarity

Thanks!

2

The bayestestR package has a mediation() function (and a mediation Vignette) that summarizes Bayesian mediation models, including models fit with brms. The function code might also be helpful in case you need to do some custom programming to get what you need.