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

Interfaces brms
1

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 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:

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

                               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:

                               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:

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:

1 Like
2

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.

1 Like
3

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) 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…

4

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.

1 Like