# Counterfactual prediction in a mediation model without all predictors
**URL:** https://discourse.mc-stan.org/t/counterfactual-prediction-in-a-mediation-model-without-all-predictors/36029
**Category:** brms
**Tags:** posterior-predictive, brms
**Created:** [July 28, 2024, 4:13am UTC](https://discourse.mc-stan.org/t/counterfactual-prediction-in-a-mediation-model-without-all-predictors/36029 "2024-07-28T04:13:55Z")
**Posts on this page:** 2
**Page:** 1
### Author:  [@Uri](https://discourse.mc-stan.org/u/Uri)
#### Post date: [July 28, 2024, 4:13am UTC](https://discourse.mc-stan.org/t/counterfactual-prediction-in-a-mediation-model-without-all-predictors/36029/1 "2024-07-28T04:13:55Z")
- 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](https://xcelab.net/rm/)” 2nd edition book and [S. Kurz’s translation](https://bookdown.org/content/4857/) 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:
```no-highlight
# 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).
```no-highlight
# 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:
```no-highlight
predict(bh5.2, newdata=tibble(m=seq(from=-3, to=3, length.out=30), resp="d") # resp for clarity
```
Thanks!
---
### Author:  [@joels](https://discourse.mc-stan.org/u/joels)
#### Post date: [July 29, 2024, 9:56pm UTC](https://discourse.mc-stan.org/t/counterfactual-prediction-in-a-mediation-model-without-all-predictors/36029/2 "2024-07-29T21:56:44Z")
The `bayestestR` package has a `mediation()` function (and a [mediation Vignette](https://easystats.github.io/bayestestR/articles/mediation.html)) that summarizes Bayesian mediation models, including models fit with `brms`. The [function code](https://github.com/easystats/bayestestR/blob/main/R/mediation.R) might also be helpful in case you need to do some custom programming to get what you need.