Function to get means of a 2x3 ANOVA model

Hi all,

please imagine a simple 3x2 design. E.g. intervention (no treatment, normal dose, double dose) where we expect the effect of it to differ between sexes. If I fit a simple model (e.g. in brms or rstan), then I can calculate for each subgroup (no treament & female, no treatment &male, normal dose & male etc…) the expected value \mu from the posterior distributions. Then, I can use those to calculate differences between groups as I wish.

My question is: is there already a function that gives you the expected value mu for each subgroup? I know how to calculate them manually (see a reproducible example below for what I mean using mtcars dataset).


# import data
d <- mtcars 
d$cyl <- as.factor(d$cyl)
d$vs <- as.factor(d$vs)
# specify model 
f <- bf(mpg ~ cyl * vs)
m <- brm(
  data = d,
  formula = f
# calculate mu for the different subgroups 
posterior_means = as_draws_df(m) %>%
    cyl4vs0 = b_Intercept,
    cyl4vs1 = b_Intercept + b_vs1 ,
    cyl6vs0 = b_Intercept + b_cyl6 ,
    cyl6vs1 = b_Intercept + b_cyl6 + `b_cyl6:vs1`, 
    cyl8vs0 = b_Intercept + b_cyl8 ,
    cyl8vs1 = b_Intercept + b_cyl8 + `b_cyl8:vs1`) %>% 

# now we can calcuate what the difference in mu is between the subgroups if we 
# wanted to know, e.g. posterior difference of cyl4vs0 and cyl6vs1
test <- posterior_means$cyl4vs0 - posterior_means$cyl6vs1

I think the solution is to use the function posterior_linpred and feed it a new dataframe with the required conditions. Please correct me if I am wrong. But I compared it to the manual calculation and they are identical.

1 Like

It works for your Gaussian model with identity link, but typically you want to use posterior_epred for this. Check out the brms manual and look at the difference between posterior_linpred and posterior_epred. The former returns draws of the linear predictor before applying any link function, while the latter draws of the expected value of the posterior predictive distribution.