- 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!