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

Interfaces
1

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
2

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)
head(as.mcmc(means)[[1]])
#> 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)[[1]])
#> [1] 1000    4

colnames(head(as.mcmc(means)[[1]]))
#> "factor1 a, factor2 y" "factor1 b, factor2 y" "factor1 a, factor2 z" "factor1 b, factor2 z"
1 Like