Standard deviations don't add up; variances do: If $X_1$ and $X_2$ are uncorrelated with std. deviation $s_1$ and $s_2$ respectively, then $X_1 + X_2$ has variance $(s_1^2 + s_2^2)$ and std. deviation $(s_1^2 + s_2^2)^{\text{½}}$.
So this is wrong:
rnorm(n, sd = sd_actor_Intercept + sd_block_Intercept)
Here's how to sample the combined effect of a new actor and a new block:
rnorm(n, sd = sqrt(sd_actor_Intercept^2 + sd_block_Intercept^2))
However, there's an easier way: you can use the brms::posterior_epred function to draw from the expected value of the posterior predictive distribution.
For illustration, I calculate the average marginal effects both ways using the prosocial chimpanzees example from Statistical Rethinking.
library("posterior") # To summarise posterior samples with `summarise_draws`.
library("tidyverse")
library("brms")
data(chimpanzees, package = "rethinking")
b12.5 <-
brm(
# In this model, actors and blocks are independent and therefore uncorrelated.
pulled_left ~ prosoc_left * condition + (1 | actor) + (1 | block),
family = bernoulli,
data = chimpanzees,
backend = "cmdstanr"
)
Method 1: Compute average marginal effects "by hand":
as_draws_df(b12.5) |>
# First we add the "fixed effects" ie. the population effects,
# for the four combinations of condition and prosocial option.
mutate(
`0/0` = b_Intercept,
`0/1` = b_Intercept + b_condition,
`1/0` = b_Intercept + b_prosoc_left,
`1/1` = b_Intercept + b_prosoc_left + `b_prosoc_left:condition`
) |>
# Then we draw random effects for a new actor and a new block.
mutate(
re_actor = rnorm(n(), sd = sd_actor__Intercept),
re_block = rnorm(n(), sd = sd_block__Intercept),
# Or equivalently:
re = rnorm(n(), sd = sqrt(sd_actor__Intercept^2 + sd_block__Intercept^2))
) |>
# Add the fixed and random effects, transform to the probability scale.
# NOTE: For convenience, we re-use the random effect draws, `re`,
# for each combination of the fixed effects.
mutate(
across(`0/0`:`1/1`, \(fe) inv_logit_scaled(fe + re))
) |>
# Summarise the posterior samples.
select(`0/0`:`1/1`) |>
summarise_draws()
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0/0 0.568 0.620 0.335 0.447 0.0273 0.992 1.00 2975. 3519.
#> 2 0/1 0.512 0.515 0.339 0.489 0.0177 0.988 1.00 3053. 3115.
#> 3 1/0 0.650 0.756 0.320 0.318 0.0483 0.996 1.00 3048. 3414.
#> 4 1/1 0.684 0.802 0.312 0.266 0.0602 0.997 1.00 2936. 3656.
Method 2: Compute average marginal effects using the {brms} machinery.
# All four combination of the fixed effects.
new_d <- crossing(prosoc_left = 0:1, condition = 0:1)
new_d
#> A tibble: 4 × 2
#> prosoc_left condition
#> <int> <int>
#> 1 0 0
#> 2 0 1
#> 3 1 0
#> 4 1 1
b12.5 |>
posterior_epred(
newdata = new_d,
# We want to draw new levels for the grouping factors, actor and block.
allow_new_levels = TRUE,
# And we want to generate these draws by sampling from the (multivariate)
# normal distribution of the actor and block variables.
sample_new_levels = "gaussian"
) |>
summarise_draws()
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 ...1 0.575 0.631 0.335 0.440 0.0209 0.994 1.00 3326. 3310.
#> 2 ...2 0.519 0.530 0.339 0.486 0.0138 0.990 1.00 3393. 3232.
#> 3 ...3 0.656 0.762 0.320 0.311 0.0392 0.996 1.00 3373. 3240.
#> 4 ...4 0.639 0.741 0.324 0.337 0.0336 0.996 1.00 3359. 3193.
I'm not happy that the variable names are ...1, ...2, ... but I'm sure there's an elegant way to fix this.
See also this blog post by Andrew Heiss: A guide to correctly calculating posterior predictions and average marginal effects with multilievel Bayesian models.
