Comparison of Bayesian and frequentist MMRM fits

About

This vignette provides an example comparison of a Bayesian MMRM fit, obtained by brms.mmrm::brm_model(), and a frequentist MMRM fit, obtained by mmrm::mmrm(). An overview of parameter estimates and differences by type of MMRM is given in the summary (Table 4) at the end.

Prerequisites

Load required packages:

> packages <- c(
+   "dplyr",
+   "tidyr",
+   "ggplot2",
+   "gt",
+   "gtsummary",
+   "purrr",
+   "parallel",
+   "brms.mmrm",
+   "mmrm",
+   "emmeans",
+   "posterior"
+ )
> invisible(lapply(packages, library, character.only = TRUE))

Set seed:

> set.seed(123L)

Data

Pre-processing

The fev_dat dataset is used that is contained in the mmrm-package:

> data(fev_data, package = "mmrm")

It is an artificial (simulated) dataset of a clinical trial investigating the effect of an active treatment on FEV1 (forced expired volume in one second), compared to placebo. FEV1 is a measure of how quickly the lungs can be emptied and low levels may indicate chronic obstructive pulmonary disease (COPD).

The dataset is a tibble with 800 rows and 7 variables:

USUBJID (subject ID),
AVISIT (visit number),
ARMCD (treatment, TRT or PBO),
RACE (3-category race),
SEX (sex),
FEV1_BL (FEV1 at baseline, %),
FEV1 (FEV1 at study visits),
WEIGHT (weighting variable).

In this analysis, change from baseline in FEV1 is used as an endpoint and a corresponding new variable (FEV1_CHG) is created.

Create and preprocess dataset for MMRM analysis:

> fev_data <- fev_data |>
+   mutate("FEV1_CHG" = FEV1 - FEV1_BL)
> fev_data <- brm_data(
+   data = fev_data,
+   outcome = "FEV1_CHG",
+   role = "change",
+   group = "ARMCD",
+   time = "AVISIT",
+   patient = "USUBJID",
+   baseline = "FEV1_BL",
+   reference_group = "PBO",
+   covariates = c("RACE", "SEX")
+ ) |>
+   mutate(ARMCD = factor(ARMCD), AVISIT = factor(AVISIT))

First rows of dataset:

> head(fev_data) |>
+   gt() |>
+   tab_caption(caption = md("Table 1. First rows of dataset."))

Table 1. First rows of dataset.
FEV1_CHG	FEV1_BL	ARMCD	AVISIT	USUBJID	RACE	SEX
NA	45.02477	PBO	VIS1	PT2	Asian	Male
-13.569552	45.02477	PBO	VIS2	PT2	Asian	Male
-8.145878	45.02477	PBO	VIS3	PT2	Asian	Male
3.783324	45.02477	PBO	VIS4	PT2	Asian	Male
NA	43.50070	PBO	VIS1	PT3	Black or African American	Female
-7.513705	43.50070	PBO	VIS2	PT3	Black or African American	Female

Descriptive statistics

Table of baseline characteristics:

> fev_data |>
+   select(ARMCD, USUBJID, SEX, RACE, FEV1_BL) |>
+   distinct() |>
+   select(-USUBJID) |>
+   tbl_summary(
+     by = c(ARMCD),
+     statistic = list(
+       all_continuous() ~ "{mean} ({sd})",
+       all_categorical() ~ "{n} / {N} ({p}%)"
+     )
+   ) |>
+   modify_caption("Table 2. Baseline characteristics.")

Table 2. Baseline characteristics.
Characteristic	PBO, N = 105¹	TRT, N = 95¹
SEX
Male	50 / 105 (48%)	44 / 95 (46%)
Female	55 / 105 (52%)	51 / 95 (54%)
RACE
Asian	38 / 105 (36%)	32 / 95 (34%)
Black or African American	46 / 105 (44%)	29 / 95 (31%)
White	21 / 105 (20%)	34 / 95 (36%)
FEV1_BL	40 (9)	40 (9)
¹ n / N (%); Mean (SD)

Table of change from baseline in FEV1 over 52 weeks:

> fev_data |>
+   pull(AVISIT) |>
+   unique() |>
+   sort() |>
+   purrr::map(
+     .f = ~ fev_data |>
+       filter(AVISIT %in% .x) |>
+       tbl_summary(
+         by = ARMCD,
+         include = FEV1_CHG,
+         type = FEV1_CHG ~ "continuous2",
+         statistic = FEV1_CHG ~ c("{mean} ({sd})",
+                                  "{median} ({p25}, {p75})",
+                                  "{min}, {max}"),
+         label = list(FEV1_CHG = paste("Visit ", .x))
+       )
+   ) |>
+   tbl_stack(quiet = TRUE) |>
+   modify_caption("Table 3. Change from baseline.")

Table 3. Change from baseline.
Characteristic	PBO, N = 105	TRT, N = 95
Visit VIS1
Mean (SD)	-8 (9)	-2 (10)
Median (IQR)	-9 (-16, 0)	-4 (-9, 7)
Range	-26, 12	-24, 20
Unknown	37	29
Visit VIS2
Mean (SD)	-3 (8)	2 (9)
Median (IQR)	-3 (-10, 1)	2 (-3, 8)
Range	-20, 15	-22, 23
Unknown	36	24
Visit VIS3
Mean (SD)	2 (8)	5 (9)
Median (IQR)	2 (-2, 8)	6 (0, 11)
Range	-15, 20	-19, 30
Unknown	34	37
Visit VIS4
Mean (SD)	8 (12)	13 (13)
Median (IQR)	6 (1, 18)	12 (5, 22)
Range	-20, 39	-14, 47
Unknown	38	28

> fev_data |>
+   group_by(ARMCD) |>
+   ggplot(aes(x = AVISIT, y = FEV1_CHG, fill = factor(ARMCD))) +
+   geom_hline(yintercept = 0, col = "grey", size = 1.2) +
+   geom_boxplot(na.rm = TRUE) +
+   labs(x = "Visit",
+        y = "Change from baseline in FEV1",
+        fill = "Treatment") +
+   scale_fill_manual(values = c("darkgoldenrod2", "coral2")) +
+   theme_bw()

Figure 1. Change from baseline in FEV1 over 4 visit time points.

Fitting MMRMs

Bayesian model

Model formula:

> b_mmrm_formula <- brm_formula(
+   data = fev_data,
+   intercept = TRUE,
+   baseline = TRUE,
+   group = FALSE,
+   time = TRUE,
+   baseline_time = FALSE,
+   group_time = TRUE,
+   correlation = "unstructured"
+ )
> print(b_mmrm_formula)
FEV1_CHG ~ FEV1_BL + ARMCD:AVISIT + AVISIT + RACE + SEX + unstr(time = AVISIT, gr = USUBJID) 
sigma ~ 0 + AVISIT

Fit model (using multiple cores, if available, to avoid long runtime):

> n_cores <- ifelse(
+   test = detectCores() > 1,
+   yes = min(2, detectCores()),
+   no = 1
+ )

> b_mmrm_fit <- brm_model(
+   data = filter(fev_data, !is.na(FEV1_CHG)),
+   formula = b_mmrm_formula,
+   iter = 1500,
+   chains = max(1, n_cores),
+   cores = n_cores,
+   refresh = 0
+ )

Summary of results:

> summary(b_mmrm_fit)
 Family: gaussian 
  Links: mu = identity; sigma = log 
Formula: FEV1_CHG ~ FEV1_BL + ARMCD:AVISIT + AVISIT + RACE + SEX + unstr(time = AVISIT, gr = USUBJID) 
         sigma ~ 0 + AVISIT
   Data: data[!is.na(data[[attr(data, "brm_outcome")]]), ] (Number of observations: 537) 
  Draws: 2 chains, each with iter = 1500; warmup = 750; thin = 1;
         total post-warmup draws = 1500

Correlation Structures:
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
cortime(VIS1,VIS2)     0.36      0.08     0.19     0.52 1.00     2550     1229
cortime(VIS1,VIS3)     0.14      0.10    -0.06     0.31 1.00     3032     1305
cortime(VIS2,VIS3)     0.04      0.10    -0.16     0.23 1.00     1928     1012
cortime(VIS1,VIS4)     0.17      0.11    -0.06     0.39 1.00     2638     1262
cortime(VIS2,VIS4)     0.11      0.09    -0.06     0.28 1.00     2404     1204
cortime(VIS3,VIS4)     0.01      0.10    -0.17     0.21 1.01     2579      855

Population-Level Effects: 
                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                     24.36      1.49    21.47    27.06 1.00     2163     1045
FEV1_BL                       -0.84      0.03    -0.90    -0.78 1.00     2979     1213
AVISITVIS2                     4.80      0.86     3.13     6.49 1.00     1604     1119
AVISITVIS3                    10.37      0.88     8.61    12.12 1.00     1338      877
AVISITVIS4                    15.16      1.30    12.62    17.79 1.00     1755     1273
RACEBlackorAfricanAmerican     1.40      0.58     0.29     2.54 1.00     1983     1032
RACEWhite                      5.44      0.61     4.31     6.64 1.00     2341     1238
SEXFemale                      0.35      0.51    -0.63     1.36 1.00     2932      965
AVISITVIS1:ARMCDTRT            3.98      1.10     1.78     6.10 1.00     1586      960
AVISITVIS2:ARMCDTRT            3.93      0.86     2.22     5.53 1.00     2091      838
AVISITVIS3:ARMCDTRT            2.98      0.67     1.66     4.33 1.00     3024     1304
AVISITVIS4:ARMCDTRT            4.46      1.66     1.22     7.72 1.01     2110     1307
sigma_AVISITVIS1               1.83      0.06     1.72     1.96 1.00     2312     1088
sigma_AVISITVIS2               1.59      0.06     1.48     1.71 1.00     2553     1228
sigma_AVISITVIS3               1.33      0.06     1.20     1.45 1.00     2458     1188
sigma_AVISITVIS4               2.28      0.06     2.17     2.40 1.00     2565     1006

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Frequentist model

Fit model:

> f_mmrm_fit <- mmrm::mmrm(
+   formula = FEV1_CHG ~ FEV1_BL + AVISIT + ARMCD:AVISIT + RACE + SEX +
+     us(AVISIT | USUBJID),
+   data = fev_data
+ )