I’m trying to obtain marginal means from a model fitted with `brms`

. I know I can do so easily using the `emmeans`

package. But I can also do this by getting the draws manually, and then calculate the averages. The problem is that I’m getting similar means, but rather different variances in the posterior distributions of the marginal means.

I can demonstrate with a simple example, actually using `rstanarm`

for convenience. Here I do the two calculations, and plot them side by side. The `emmeans`

package clearly has wider posteriors, but I don’t understand why.

```
library(rstanarm)
library(emmeans)
library(tidybayes)
library(dplyr)
library(modelr)
mtcars2 <- mutate(mtcars, cyl = factor(cyl), vs = factor(vs))
m1 <- stan_lmer(qsec ~ vs + (vs | cyl), mtcars2)
summary(m1)
emm1 <- emmeans(m1, ~vs) %>%
gather_emmeans_draws()
emm2 <- data_grid(mtcars2, vs, cyl) %>%
add_fitted_draws(m1) %>%
group_by(.draw, vs) %>%
summarise(.value = mean(.value))
emms <- bind_rows(emmeans = emm1, manual = emm2, .id = 'method')
ggplot(emms, aes(.value, vs)) +
geom_halfeyeh() +
facet_grid(method ~ .)
```

- Operating System: Windows 10
- brms Version: 2.6.0