# Emmeans or hypothesis for getting comparisons with Mixed Model for Repeated Measures (MMRM) fitted with brms

**URL:** <https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831>\
**Category:** brms\
**Tags:** brms\
**Created:** [December 20, 2022, 9:50am UTC](https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831 "2022-12-20T09:50:37Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Bjoern](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bjoern/32/5429_2.png) [@Bjoern](https://discourse.mc-stan.org/u/Bjoern)\
**Post date:** [December 20, 2022, 9:50am UTC](https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831/1 "2022-12-20T09:50:37Z")

</div>

I am trying to understand whether I should use `hypothesis` (I tried with and without `robust=T`) from `brms` or `emmeans` + `pairs` or `contrast` from the `emmeans` package to get treatment comparisons at different visits from a Mixed Model for Repeated Measures (MMRM) fitted with `brms`.

The example data is a simulated randomized trial with 3 doses of a drug compared with a placebo, with the continuous outcome measured at baseline and at 4 post-baseline visits. We want treatment comparisons (primarily vs. placebo) at each post-baseline visit using a typical MMRM model (using the `unstr` covariance matrix just added in [this branch](https://github.com/paul-buerkner/brms/tree/unstr-corr) that just got merged to master).

I’m asking because I get some small differences (pretty small, but perhaps one approach is just obviously the right way to do this?), namely with `emmeans::contrast`:

```plaintext
AVISIT = visit1:
 contrast estimate lower.HPD upper.HPD
 TRT01P10 - TRT01P0 0.050 -0.2331 0.359
 TRT01P20 - TRT01P0 0.119 -0.1842 0.437
 TRT01P40 - TRT01P0 0.292 -0.0637 0.647

AVISIT = visit2:
 contrast estimate lower.HPD upper.HPD
 TRT01P10 - TRT01P0 0.494 0.1408 0.832
 TRT01P20 - TRT01P0 0.535 0.1937 0.878
 TRT01P40 - TRT01P0 0.699 0.3215 1.089

AVISIT = visit3:
 contrast estimate lower.HPD upper.HPD
 TRT01P10 - TRT01P0 0.756 0.3742 1.096
 TRT01P20 - TRT01P0 0.800 0.4352 1.130
 TRT01P40 - TRT01P0 0.924 0.5069 1.329

AVISIT = visit4:
 contrast estimate lower.HPD upper.HPD
 TRT01P10 - TRT01P0 0.641 0.2109 1.051
 TRT01P20 - TRT01P0 1.077 0.6580 1.483
 TRT01P40 - TRT01P0 0.846 0.3558 1.315

Point estimate displayed: median 
HPD interval probability: 0.95 

```

vs. with `brms::hypothesis` with `robust=T`

```plaintext
                               Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (TRT01P10) < 0 0.04997481 0.1498573 -0.23652611 0.3573566 0.584786054 0.36900     
2 (TRT01P10+`AVISITvisit2:TRT01P10`) < 0 0.49446603 0.1798226 0.13724712 0.8308284 0.002757583 0.00275     
3 (TRT01P10+`AVISITvisit3:TRT01P10`) < 0 0.75551675 0.1834758 0.38120334 1.1083434 0.000000000 0.00000     
4 (TRT01P10+`AVISITvisit4:TRT01P10`) < 0 0.64062921 0.2176190 0.22142103 1.0649631 0.002004008 0.00200     
5 (TRT01P20) < 0 0.11927446 0.1552015 -0.18766045 0.4340984 0.280409731 0.21900     
6 (TRT01P20+`AVISITvisit2:TRT01P20`) < 0 0.53504101 0.1784003 0.19045258 0.8752969 0.001753068 0.00175     
7 (TRT01P20+`AVISITvisit3:TRT01P20`) < 0 0.79969931 0.1736051 0.44967037 1.1500322 0.000000000 0.00000     
8 (TRT01P20+`AVISITvisit4:TRT01P20`) < 0 1.07718263 0.2135734 0.65742548 1.4814281 0.000000000 0.00000     
9 (TRT01P40) < 0 0.29201577 0.1818559 -0.06070752 0.6529470 0.057082452 0.05400     
10 (TRT01P40+`AVISITvisit2:TRT01P40`) < 0 0.69867399 0.2013159 0.31675064 1.0847861 0.000000000 0.00000     
11 (TRT01P40+`AVISITvisit3:TRT01P40`) < 0 0.92408798 0.2100295 0.51222037 1.3386063 0.000000000 0.00000     
12 (TRT01P40+`AVISITvisit4:TRT01P40`) < 0 0.84635702 0.2515857 0.37197939 1.3315527 0.000000000 0.00000 

```

vs. `brms::hypothesis` with `robust=F`:

```plaintext
                               Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (TRT01P10) < 0 0.04986414 0.1506882 -0.23652611 0.3573566 0.584786054 0.36900     
2 (TRT01P10+`AVISITvisit2:TRT01P10`) < 0 0.49150227 0.1800661 0.13724712 0.8308284 0.002757583 0.00275     
3 (TRT01P10+`AVISITvisit3:TRT01P10`) < 0 0.75078906 0.1858033 0.38120334 1.1083434 0.000000000 0.00000     
4 (TRT01P10+`AVISITvisit4:TRT01P10`) < 0 0.64224199 0.2188358 0.22142103 1.0649631 0.002004008 0.00200     
5 (TRT01P20) < 0 0.12040181 0.1579221 -0.18766045 0.4340984 0.280409731 0.21900     
6 (TRT01P20+`AVISITvisit2:TRT01P20`) < 0 0.53414896 0.1770733 0.19045258 0.8752969 0.001753068 0.00175     
7 (TRT01P20+`AVISITvisit3:TRT01P20`) < 0 0.79689191 0.1774612 0.44967037 1.1500322 0.000000000 0.00000     
8 (TRT01P20+`AVISITvisit4:TRT01P20`) < 0 1.07723592 0.2113137 0.65742548 1.4814281 0.000000000 0.00000     
9 (TRT01P40) < 0 0.29439702 0.1811998 -0.06070752 0.6529470 0.057082452 0.05400     
10 (TRT01P40+`AVISITvisit2:TRT01P40`) < 0 0.69977352 0.1968975 0.31675064 1.0847861 0.000000000 0.00000     
11 (TRT01P40+`AVISITvisit3:TRT01P40`) < 0 0.92335447 0.2135107 0.51222037 1.3386063 0.000000000 0.00000     
12 (TRT01P40+`AVISITvisit4:TRT01P40`) < 0 0.84717469 0.2472672 0.37197939 1.3315527 0.000000000 0.00000

```

Here’s my data simulation & model fitting code:

```plaintext
library(tidyverse)
library(brms)
library(posterior)
library(mvtnorm)
library(dqrng)
library(emmeans)

#### Simulate some data ----

set.seed(1234)
# Correlation matrix between visits (baseline + 4 post-baseline visits)
corr_matrix <- matrix(c(1, 0.6, 0.48, 0.4, 0.375,
                        0.6, 1, 0.6, 0.48, 0.4,
                        0.48, 0.6, 1, 0.6, 0.48,
                        0.4, 0.48, 0.6, 1, 0.6,
                        0.375, 0.4, 0.48, 0.6, 1),
                      nrow=5,ncol=5,byrow = T)

# Standard deviations by visit (baseline + 4 post-baseline visits)
sds <- sqrt(c(0.75, 0.8, 0.85, 0.95, 1.1))

cov_matrix <- t(sds * t(sds*corr_matrix))

apply_colnames <- function(x,y){
  colnames(x) <- y
  return(x)
}

set.seed(1234)

# Simulate from multivariate normal for control group 
# (before adding treatment effect later)
# We simulate 1000 patients and then apply inclusion criteria and keep the
# first 200 that meet them.
simulated_data <- rmvnorm(n=1000, mean=rep(0, 5), sigma = cov_matrix) %>%
  # turn into tibble
  apply_colnames(c("BASE", paste0("visit", 1:4))) %>%
  as_tibble() %>%
  # Apply inclusion criteria and keep first 200 patients
  filter(BASE>0) %>%
  filter(row_number()<=200) %>%
  # Assign subject ID, treatment group and create missing data
  mutate(USUBJID = row_number(),
         TRT01P = dqsample(x=c(0L, 10L, 20L, 40L), size=200, replace=T),
         # Simulate dropouts
         dropp2 = plogis(visit1-2),
         dropp3 = plogis(visit2-2),
         dropp4 = plogis(visit3-2),
         dropv = case_when(runif(n=n())<dropp2 ~ 2L,
                           runif(n=n())<dropp3 ~ 3L,
                           runif(n=n())<dropp4 ~ 4L,
                           TRUE ~ 5L),
         visit2 = ifelse(dropv<=2L, NA_real_, visit2),
         visit3 = ifelse(dropv<=3L, NA_real_, visit3),
         visit4 = ifelse(dropv<=3L, NA_real_, visit4)) %>%
  dplyr::select(-dropp2, -dropp3, -dropp4, -dropv) %>%
  # Turn data into long-format
  pivot_longer(cols=starts_with("visit"), names_to = "AVISIT", values_to="AVAL") %>%
  mutate(
    # Assign visit days
    ADY = case_when(AVISIT=="visit1" ~ 2L*7L,
                    AVISIT=="visit2" ~ 4L*7L,
                    AVISIT=="visit3" ~ 8L*7L,
                    AVISIT=="visit4" ~ 12L*7L),
    # Assume rising treatment effect over time (half there by week 3) with an 
    # Emax dose response (ED50 = 5 mg)
    AVAL = AVAL + 0.9 * (ADY/7)^3/((ADY/7)^3+3^3) * TRT01P/(TRT01P+5),
    # Change from baseline = value - baseline
    CHG = AVAL - BASE,
    TRT01P=factor(TRT01P)) %>%
  relocate(USUBJID, TRT01P, AVISIT, ADY, AVAL, CHG, BASE) %>%
  # Discard missing data
  filter(!is.na(AVAL))

#### Fit a model ----

brmfit1 <- brm(
  formula=bf( CHG ~ 1 + AVISIT + TRT01P + BASE + TRT01P:AVISIT + BASE:AVISIT,  
              autocor=~unstr(time=AVISIT, gr=USUBJID), sigma ~ AVISIT),
  data = simulated_data,
  control = list(adapt_delta=0.99),
  # Explicitly specifiyng quite weak priors (given the variability in the data).
  prior= prior(class="b", student_t(3, 0, 10)) + 
    prior(class="Intercept", student_t(3, 0, 25)))

#### Compare emmeans with hypothesis ----

emm1 <- emmeans(brmfit1, ~ TRT01P | AVISIT, weights="proportional")
#pairs(emm1) # Also works
contrast(emm1, adjust="none", method="trt.vs.ctrl") 

contrasts_of_interest <- c("TRT01P10", "TRT01P10 + `AVISITvisit2:TRT01P10`", 
                           "TRT01P10 + `AVISITvisit3:TRT01P10`",
                           "TRT01P10 + `AVISITvisit4:TRT01P10`",
                           "TRT01P20", "TRT01P20 + `AVISITvisit2:TRT01P20`", 
                           "TRT01P20 + `AVISITvisit3:TRT01P20`",
                           "TRT01P20 + `AVISITvisit4:TRT01P20`",
                           "TRT01P40", "TRT01P40 + `AVISITvisit2:TRT01P40`", 
                           "TRT01P40 + `AVISITvisit3:TRT01P40`",
                           "TRT01P40 + `AVISITvisit4:TRT01P40`") %>%
  paste0(" < 0")

get_hypothesis <- function(brmfit, contrast){
  hypothesis(brmfit, contrast, alpha = 0.025, robust=T)$hypothesis
}

map_dfr(contrasts_of_interest, function(x) get_hypothesis(brmfit1, x))

```

Please also provide the following information in addition to your question:

- Operating System: Windows 10
- brms Version: Latest version as of Friday of [GitHub - paul-buerkner/brms at unstr-corr](https://github.com/paul-buerkner/brms/tree/unstr-corr) (brms help says “Package brms version 2.18.5”)

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [December 20, 2022, 6:36pm UTC](https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831/2 "2022-12-20T18:36:22Z")

</div>

Can you explain which differences you are refering to exactly? Some differences in uncertainty intervals are expected because of emmeans uses HPDs by default while brms uses (quantile-based) credible intervals.

---

<div class="post-metadata">

**Author:** ![Bjoern](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bjoern/32/5429_2.png) [@Bjoern](https://discourse.mc-stan.org/u/Bjoern)\
**Post date:** [December 21, 2022, 7:52am UTC](https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831/3 "2022-12-21T07:52:43Z")

</div>

Ah, thanks, I think that’s the explanation. The point estimates all seem the same with `robust=T` and all the differences were small differences in the credible intervals that could very plausibly be explained by HPD vs. quantile.

To confirm/make it the same if desired, as far as I could see, there’s no easy way to get HPD intervals out of hypothesis (when I looked it seems like quantiles are hard-coded as the method [here](https://github.com/paul-buerkner/brms/blob/master/R/hypothesis.R)) and I couldn’t see a quantile options for `emmeans`, either (but didn’t dig into their codebase).

I’ll admit for many basic questions the `emmeans` syntax seems very convenient, but it’s kind of funny if that then comes with an automatic unavoidable choice of HPD.

Perhaps I shouldn’t exaggerate the differences between HPD and quantile-based in pracctice, and most models I work with have nice unimodal posteriors that are kind of approximately normal on the right scale. In that case there’ll not be much of a difference (and then we don’t need to worry too much about invariance under transformations either, which is a potential thing with non-linear link functions). Thoughts on that are certainly welcome…

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [December 21, 2022, 2:08pm UTC](https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831/4 "2022-12-21T14:08:23Z")

</div>

The emmeans interface is much more powerful than that of hypothesis. With respect to uncertainty intervals I prefer CIs over HPDs because the former can be more robustly estimated via draws. Not sure if we can use CIs in emmeans but if not, it would definitely be a good addition to the package.

---

<div class="post-metadata">

**Author:** ![Mat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mat/32/20100_2.png) [@Mat](https://discourse.mc-stan.org/u/Mat)\
**Post date:** [July 30, 2025, 7:06am UTC](https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831/5 "2025-07-30T07:06:29Z")

</div>

FYI, this feature has now been implemented in `emmeans:`

> <https://github.com/rvlenth/emmeans/issues/538>
>
> Hi,
> 
> Currently, \`emmeans\` returns HPD intervals for Bayesian models, but there’s… no option to use quantile-based credible intervals instead (as done by \`brms::hypothesis()\`). It would be great to add a \`ci\_type\` argument, e.g., \`ci\_type = "hpd"\` (default) or \`ci\_type = "quantile"\`, to let users choose.
> 
> Paul Bürkner, the creator of \`brms\`, also suggested this would be a useful addition:
> \[https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831/4\](https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831/4)
> 
> Thanks!

---

<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:** [July 30, 2025, 8:20am UTC](https://discourse.mc-stan.org/t/emmeans-or-hypothesis-for-getting-comparisons-with-mixed-model-for-repeated-measures-mmrm-fitted-with-brms/29831/6 "2025-07-30T08:20:19Z")

</div>

FWIW, this was already possible, in a slightly more roundabout way, with a tidybayes-emmeans combo:

```r
fit %>% emmeans(...) %>% 
  contrast(...) %>%
  gather_emmeans_draws() %>% 
  mean_qi(.width = .95)

```