# Anova-like summary for brm-model?

**URL:** <https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119>\
**Category:** brms\
**Tags:** brms\
**Created:** [March 20, 2025, 1:42pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119 "2025-03-20T13:42:19Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![striatum](https://avatars.discourse-cdn.com/v4/letter/s/ac8455/32.png) [@striatum](https://discourse.mc-stan.org/u/striatum)\
**Post date:** [March 20, 2025, 1:42pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/1 "2025-03-20T13:42:20Z")

</div>

The emmeans package function joint\_test can produce an anova-like summary for categorical predictors (factors), for example, the overall effect of FactorA, FactorB, and FactorA:FactorB. This can be done even for brm model (from brmsfit), which might be useful even though only to help “traditionalists” and/or for purely didactic purposes. However, it feels weird to have p-values rather than something “more Bayesian”, like HPD intervals etc. Thus, I was wondering how could you compute such a summary of the overall effects – e.g., a HPD interval for FactorA, not for its individual levels? Thanks!

---

<div class="post-metadata">

**Author:** ![Solomon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/solomon/32/2495_2.png) [@Solomon](https://discourse.mc-stan.org/u/Solomon)\
**Post date:** [March 20, 2025, 6:17pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/2 "2025-03-20T18:17:45Z")

</div>

Hey @striatum, would you mind showing a minimal reproducible example?

---

<div class="post-metadata">

**Author:** ![striatum](https://avatars.discourse-cdn.com/v4/letter/s/ac8455/32.png) [@striatum](https://discourse.mc-stan.org/u/striatum)\
**Post date:** [March 20, 2025, 7:09pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/3 "2025-03-20T19:09:31Z")

</div>

Thank you @Solomon for helping me! I am taking my old post example: [How to properly compare interacting levels](https://discourse.mc-stan.org/t/how-to-properly-compare-interacting-levels/20457/3). The model is slightly revised, but is nice and simple leaving aside random effects, covariates etc.

```no-highlight
# uniformly generated means
m = runif(4, 3, 18)

# random values, random means, same std. dev.
y = list()
for (i in 1:4) {
    y[[i]] = rnorm(30, m[i], 0.8)
}
y = unlist(y)

# factorial design
F1 = c(rep('A', 30), rep('B', 30))
F2 = c('I', 'J')
design_matrix = expand.grid(F1=F1, F2=F2)

# dataset
dat = cbind(design_matrix, y)

library(brms)
library(emmeans)

model1 <- brm(y ~
    F1 * F2,
    data = dat,
    chains=4, iter=4000, cores=4)

model1.emm = emmeans(model1, ~ F1 * F2)
joint_tests(model1.emm)
 # model term df1 df2 F.ratio Chisq p.value
 # F1 1 Inf 5288.825 5288.825 <.0001
 # F2 1 Inf 18.793 18.793 <.0001
 # F1:F2 1 Inf 56.745 56.745 <.0001

```

I appreciate the anova-style overall effect of F1, F2, and the interaction, but credible intervals might be more in Bayesian-spirit.

---

<div class="post-metadata">

**Author:** ![Solomon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/solomon/32/2495_2.png) [@Solomon](https://discourse.mc-stan.org/u/Solomon)\
**Post date:** [March 20, 2025, 7:40pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/4 "2025-03-20T19:40:17Z")

</div>

Hmm. I don’t think I’m of any help, here. But yes I agree, that kind of breakdown seems very non-Bayesian. I just wouldn’t do it.

---

<div class="post-metadata">

**Author:** ![striatum](https://avatars.discourse-cdn.com/v4/letter/s/ac8455/32.png) [@striatum](https://discourse.mc-stan.org/u/striatum)\
**Post date:** [March 20, 2025, 7:57pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/5 "2025-03-20T19:57:58Z")

</div>

Thanks @Solomon! So nothing meaningful that would give something like HPD for the “whole” F1 or F2?

---

<div class="post-metadata">

**Author:** ![Solomon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/solomon/32/2495_2.png) [@Solomon](https://discourse.mc-stan.org/u/Solomon)\
**Post date:** [March 21, 2025, 3:05pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/6 "2025-03-21T15:05:00Z")

</div>

To my mind, if you care about the coefficients, just describe them with point estimates and 95% intervals or perhaps even plot them. If you care about contrasts, then just compute the contrasts of interest.

---

<div class="post-metadata">

**Author:** ![striatum](https://avatars.discourse-cdn.com/v4/letter/s/ac8455/32.png) [@striatum](https://discourse.mc-stan.org/u/striatum)\
**Post date:** [March 24, 2025, 9:59am UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/7 "2025-03-24T09:59:14Z")

</div>

A follow-up, if I may, @Solomon (and with apologies if this is plain dumb): what I am asking about is the so-called **omnibus test** of the whole factor effect, not the differences/contrasts or the coefficients. Would you consider averages of (per)centile points to be meaningful quantities? If so, I could get an average value at 95% low/high, i.e., something like a credible interval for a given factor.

---

<div class="post-metadata">

**Author:** ![amynang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/amynang/32/19114_2.png) [@amynang](https://discourse.mc-stan.org/u/amynang)\
**Post date:** [March 24, 2025, 10:52am UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/8 "2025-03-24T10:52:45Z")

</div>

It’s been a while, but if I remember correctly, these tests are based on log likelihood ratios of different models; one that contains the focal term versus one that doesn’t. Can you clarify when you say that you would rather have an interval rather than a p-value, what would it be an interval of?

+1 for

> [@Solomon](#):
>
> I just wouldn’t do it.

---

<div class="post-metadata">

**Author:** ![striatum](https://avatars.discourse-cdn.com/v4/letter/s/ac8455/32.png) [@striatum](https://discourse.mc-stan.org/u/striatum)\
**Post date:** [March 24, 2025, 11:25am UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/9 "2025-03-24T11:25:53Z")

</div>

@Solomon, can you summarise briefly why you would not do that?

---

<div class="post-metadata">

**Author:** ![Solomon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/solomon/32/2495_2.png) [@Solomon](https://discourse.mc-stan.org/u/Solomon)\
**Post date:** [March 24, 2025, 1:57pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/10 "2025-03-24T13:57:51Z")

</div>

I’m not a fan of NHST. Though one can do NHST as a Bayesian, to my eye Bayesian software such as brms makes it particularly convenient to switch from hypothesis testing to predictions and effect sizes. For more on this, see Figure 1 from [Kruschke & Liddell (2018)](https://doi.org/10.3758/s13423-016-1221-4). In the context of that figure, I’m a bottom-row kind of guy.

---

<div class="post-metadata">

**Author:** ![FeliOtte](https://avatars.discourse-cdn.com/v4/letter/f/bb73d2/32.png) [@FeliOtte](https://discourse.mc-stan.org/u/FeliOtte)\
**Post date:** [April 2, 2025, 2:22pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/11 "2025-04-02T14:22:41Z")

</div>

Hi there!  
I was after the same thing for a while and ultimately decided to report conditional effects instead. They also give you an idea of the overall effect and you don’t have to leave the Bayesian sphere to do so. You can get actual numbers instead of a plot from conditional\_effects like this:

```no-highlight
condeff <- conditional_effects(mdl_name, 
categorical = TRUE, #only if outcome variable is categorical, adapt as appropriate
re_formula=NULL) #means that random effects are included in the calculation

```

I don’t know if that’s a possible alternative in your case, but maybe it’s something to consider. Good luck either way!

---

<div class="post-metadata">

**Author:** ![jsocolar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jsocolar/32/2486_2.png) [@jsocolar](https://discourse.mc-stan.org/u/jsocolar)\
**Post date:** [April 2, 2025, 6:32pm UTC](https://discourse.mc-stan.org/t/anova-like-summary-for-brm-model/39119/12 "2025-04-02T18:32:26Z")

</div>

In a Bayesian setting, we generally know a priori that the effects of the different factor levels are not all jointly zero (we typically use a prior that assigns only infinitesimal probability mass to this possibility). Questions like “does this factor belong in the model” are often viewed as problems of model selection, and are often judged based on a model’s predictive performance (e.g. as assessed by cross validation). If desired you can construct a fully Bayesian test of whether inclusion of a factor variable improves the model’s predictive performance based on leave-one-out cross validation (LOO-CV). `brms` includes functionality from the `loo` package to make this quite easy to do.

Another alternative, which is not available inside of `brms` and is rarely recommended on this forum is to marginalize over a binary indicator variable that indicates which of two models is the “correct” one. This _is_ a way to assign nonzero prior mass to the possibility that the effects of all levels in a factor are jointly zero. But that’s usually not actually a realistic statement of prior beliefs, which I think is why you won’t see it mentioned much here. I personally find it nicer and more sensible to view the problem as a question about predictive performance and model selection, while acknowledging that it’s vanishingly unlikely that the population level differences between the factor levels are literally zero.