Difference betweem brms.mmrm and mmrm estimates

[yarriofultramar]

I wanted to compare results of MMRM models that was run using mmrm package with results of brms.mmrm. I do not have much experience with BRMS so this might be something obvious but when I look at summaries between frequentist model

model<-mmrm(Change~ARM+Baseline+Site_pooled+Gender+ Visit +Visit*ARM + us(Visit|USUBJIDID), data=data, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear")

and brms.mmrm:

data_mmrm <- brm_data(
  data = filter(baza, FAS==1),
  outcome = "Change",
  group = "ARM",
  time = "Visit",
  patient = "USUBJID",
  baseline = "Baseline",
  reference_group = "Placebo",
  reference_time = NULL,
  covariates = c("Site_pooled", "Gender")
)

model <- brm_model(data, formula = formula, refresh = 0)

I get totally different estimates (using summary(model)).
Although the predicted values (response) are very close to observed data.

draws <- brm_marginal_draws(model = model)
marginal_summaries<-brm_marginal_summaries(draws, level = 0.90)
marginal_summaries_res<-marginal_summaries %>%
  select(-mcse) %>%
  pivot_wider(id_cols=c("marginal", "time", "group"), names_from = "statistic", values_from = "value") %>%
  select(
    marginal, time, group,  mean, median,  lower, upper, sd
  )

I am sure I made some mistake but I do not understand why this is a case? Should I set priors manually? If so how would I know what would be their parameters?
I tried mixed effect model:

model_lmer<-brm(Change ~ Baseline + ARM + ARM:Visit + Visit + Site_pooled + Gender + (Visit|USUBJID), data=data)

And here estimates were much closer to frequentist mmrm.

I appreciate your hard work and would be very obliged for some pointers.

[wlandau]

Would you post a reproducible example? It would help to be able to run the code and see what you see. Otherwise, it's hard to say what might be happening.

I get totally different estimates (using summary(model)).

Are the parameterizations the same between models?

[wlandau]

Looks like development brms (version 2.22.9) has the same issue. I will work on a reprex with brms only (without brms.mmrm) and then post an issue to https://github.com/paul-buerkner/brms/issues if I can prove the problem is there.

[yarriofultramar]

Thank you very much! It is first time for me that I see such an issue in statistical program. I tried brms with mixed model (VISIT|USUBID) and it gave much more sensible results.

[wlandau]

I have a reprex that uses only brms (without brms.mmrm). It turns out unstructured covariance breaks when the rows of the data are sorted a certain (sensible) way.

library(brms)
library(emmeans)
library(tidyverse)
options(brms.backend = "rstan")

data(fev_data, package = "mmrm")
data_unsorted <- fev_data |>
  filter(!is.na(FEV1))
data_sorted <- data_unsorted |>
  arrange(ARMCD, AVISIT, USUBJID)

marginals_data <- data_unsorted |>
  group_by(ARMCD, AVISIT) |>
  summarize(
    mean = mean(FEV1),
    lower = mean - 1.96 * sd(FEV1) / sqrt(n()),
    upper = mean + 1.96 * sd(FEV1) / sqrt(n()),
    .groups = "drop"
  ) |>
  mutate(AVISIT = as.factor(as.character(AVISIT)))

marginals_model <- function(data) {
  marginal_unsorted <- brm(
    data = data,
    formula = FEV1 ~ 0 + ARMCD + AVISIT + unstr(time = AVISIT, gr = USUBJID),
  ) |>
    emmeans(specs = c("ARMCD", "AVISIT")) |>
    as.data.frame() |>
    rename(mean = emmean, lower = lower.HPD, upper = upper.HPD)
}

marginals <- bind_rows(
  sorted = marginals_model(data_sorted),
  unsorted = marginals_model(data_unsorted),
  data = marginals_data,
  .id = "source"
)

ggplot(marginals) +
  geom_point(
    aes(x = AVISIT, y = mean, color = source),
    position = position_dodge(width = 0.5)
  ) +
  geom_errorbar(
    aes(x = AVISIT, ymin = lower, ymax = upper, color = source),
    position = position_dodge(width = 0.5)
  ) +
  facet_wrap(~ ARMCD)

sessionInfo()
#> R version 4.4.0 (2024-04-24)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS 15.3.1
#> 
#> Matrix products: default
#> BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/Indiana/Indianapolis
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] lubridate_1.9.3 forcats_1.0.0   stringr_1.5.1   dplyr_1.1.4     purrr_1.0.2     readr_2.1.5    
#>  [7] tidyr_1.3.1     tibble_3.2.1    ggplot2_3.5.1   tidyverse_2.0.0 emmeans_1.10.4  brms_2.22.9    
#> [13] Rcpp_1.0.13-1  
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.5         tensorA_0.36.2.1     QuickJSR_1.3.1       processx_3.8.4      
#>  [5] inline_0.3.20        lattice_0.22-6       callr_3.7.6          tzdb_0.4.0          
#>  [9] ps_1.8.1             vctrs_0.6.5          tools_4.4.0          generics_0.1.3      
#> [13] curl_6.0.1           stats4_4.4.0         parallel_4.4.0       sandwich_3.1-0      
#> [17] pkgconfig_2.0.3      Matrix_1.7-0         checkmate_2.3.2      distributional_0.5.0
#> [21] RcppParallel_5.1.9   lifecycle_1.0.4      farver_2.1.2         compiler_4.4.0      
#> [25] textshaping_0.4.0    Brobdingnag_1.2-9    munsell_0.5.1        codetools_0.2-20    
#> [29] bayesplot_1.11.1     pillar_1.10.1        MASS_7.3-61          StanHeaders_2.32.10 
#> [33] bridgesampling_1.1-2 abind_1.4-8          multcomp_1.4-26      nlme_3.1-165        
#> [37] posterior_1.6.0      rstan_2.32.6         tidyselect_1.2.1     mvtnorm_1.2-5       
#> [41] stringi_1.8.4        labeling_0.4.3       splines_4.4.0        grid_4.4.0          
#> [45] colorspace_2.1-1     cli_3.6.3            magrittr_2.0.3       loo_2.8.0           
#> [49] survival_3.7-0       pkgbuild_1.4.4       TH.data_1.1-2        withr_3.0.2         
#> [53] scales_1.3.0         backports_1.5.0      estimability_1.5.1   timechange_0.3.0    
#> [57] matrixStats_1.4.1    gridExtra_2.3        ragg_1.3.3           zoo_1.8-12          
#> [61] hms_1.1.3            coda_0.19-4.1        V8_6.0.0             rstantools_2.4.0    
#> [65] rlang_1.1.4          xtable_1.8-4         glue_1.8.0           jsonlite_1.8.9      
#> [69] rstudioapi_0.16.0    R6_2.5.1             systemfonts_1.1.0

[wlandau]

I just posted paul-buerkner/brms#1747. Hopefully this will be fixed on the brms side. Not sure what brms.mmrm could do because I'm not sure if sorting by USUBJID first would fix the bug in all cases. (Nor should it be necessary.)

[wlandau]

The bug in brms was fixed in the development version (pak::pkg_install("paul-buerkner/brms")).

[yarriofultramar]

Thank you for help!