# A parameterization-agnostic way to get posterior samples of marginal means in {brms}

I am developing an R package on top of `brms`, and the fixed effect parameterization depends on arbitrary inputs from the user. My goal is to produce custom posterior summaries based on the marginal means, without having to make assumptions about the choice of underlying parameters. Is there a way to get posterior samples of the marginal means? The `emmeans` package produces posterior means and credible intervals of the marginal means, but for my purposes, it would be great to have direct access to the underlying posterior samples.

``````library(brms)
library(tidyr)
data <- expand_grid(
unit = seq_len(100),
factor1 = c("a", "b"),
factor2 = c("y", "z")
)
data\$outcome <- rnorm(n = nrow(data))
data
#> # A tibble: 400 × 4
#>     unit factor1 factor2 outcome
#>    <int> <chr>   <chr>     <dbl>
#>  1     1 a       y       -0.0438
#>  2     1 a       z       -0.506
#>  3     1 b       y        1.03
#>  4     1 b       z       -0.320
#>  5     2 a       y        0.935
#>  6     2 a       z       -1.11
#>  7     2 b       y       -1.74
#>  8     2 b       z        1.28
#>  9     3 a       y        0.426
#> 10     3 a       z       -0.0901
#> # ℹ 390 more rows
#> # ℹ Use `print(n = ...)` to see more rows

formula <- brmsformula(formula = outcome ~ factor1 * factor2)
model <- brm(formula = formula, data = data)
model # The parameterization is arbitrary.
#> ...
#> Population-Level Effects:
#>                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept            -0.09      0.09    -0.27     0.10 1.00     2776     2976
#> factor1b              0.09      0.13    -0.17     0.34 1.00     2502     2156
#> factor2z             -0.08      0.13    -0.33     0.17 1.00     2483     2678
#> factor1b:factor2z     0.05      0.18    -0.30     0.42 1.00     2241     1972
#> ...

# Would be great to access the underlying posterior samples of the
# marginal summaries below.
emmeans(object = model, specs = ~factor1:factor2)
#> factor1 factor2   emmean lower.HPD upper.HPD
#> a       y       -0.09058    -0.270    0.0924
#> b       y        0.00758    -0.166    0.1856
#> a       z       -0.16373    -0.341    0.0183
#> b       z       -0.02090    -0.202    0.1594
``````

Hmm… after reading Sophisticated models in emmeans, the answer turns out to be surprisingly simple: `as.mcmc()` just works on `emmeans()` objects produced from `brms` models!

``````means <- emmeans(object = model, specs = ~factor1:factor2)
#> Markov Chain Monte Carlo (MCMC) output:
#> Start = 1
#> End = 7
#> Thinning interval = 1
#>   factor1 a, factor2 y factor1 b, factor2 y factor1 a, factor2 z factor1 b, factor2 z
#> 1         -0.080343590          -0.03301226            0.1116829         -0.007073059
#> 2          0.123850908           0.10712421            0.1137746          0.074622100
#> 3         -0.004248606          -0.11752113            0.1320306         -0.162218602
#> 4          0.080808490           0.03389035            0.1664620          0.081042827
#> 5         -0.063756760           0.14679267            0.1415305          0.086451079
#> 6          0.067092005          -0.05948340            0.1002536         -0.177094880
#> 7         -0.006326991          -0.14891293            0.1451574         -0.107361345

dim(as.mcmc(means)[])
#>  1000    4