About

This vignette uses the bcva_data dataset from the mmrm package to compare 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 (Tables 4 and 5) at the end.

1 Prerequisites

This comparison workflow requires the following packages.

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

We set a seed for the random number generator to ensure statistical reproducibility.

> set.seed(123L)

2 Data

2.1 Pre-processing

This analysis exercise uses the bcva_data dataset contained in the mmrm package:

> data(bcva_data, package = "mmrm")

According to https://openpharma.github.io/mmrm/latest-tag/articles/mmrm_review_methods.html:

The BCVA dataset contains data from a randomized longitudinal ophthalmology trial evaluating the change in baseline corrected visual acuity (BCVA) over the course of 10 visits. BCVA corresponds to the number of letters read from a visual acuity chart.

The dataset is a tibble with 800 rows and 7 variables. The primary endpoint for the analysis is BCVA_CHG.

USUBJID (subject ID)
AVISITN (visit number, numeric)
AVISIT (visit number, factor)
ARMCD (treatment, TRT or CTL)
RACE (3-category race)
BCVA_BL (BCVA at baseline)
BCVA_CHG (BCVA change from baseline)

The rest of the pre-processing steps create factors for the study arm and visit and apply the usual checking and standardization steps of brms.mmrm::brm_data().

> bcva_data <- brm_data(
+   data = bcva_data,
+   outcome = "BCVA_CHG",
+   role = "change",
+   group = "ARMCD",
+   time = "AVISIT",
+   patient = "USUBJID",
+   baseline = "BCVA_BL",
+   reference_group = "CTL",
+   covariates = "RACE"
+ ) |>
+   mutate(ARMCD = factor(ARMCD), AVISIT = factor(AVISIT))

The following table shows the first rows of the dataset.

> head(bcva_data) |>
+   gt() |>
+   tab_caption(caption = md("Table 1. First rows of the pre-processed `bcva_data` dataset."))

Table 1. First rows of the pre-processed `bcva_data` dataset.
BCVA_CHG	BCVA_BL	ARMCD	AVISIT	USUBJID	RACE
5.058546	71.70881	CTL	VIS01	3	Asian
4.018582	71.70881	CTL	VIS02	3	Asian
3.572535	71.70881	CTL	VIS03	3	Asian
4.822669	71.70881	CTL	VIS04	3	Asian
7.348768	71.70881	CTL	VIS05	3	Asian
6.076615	71.70881	CTL	VIS06	3	Asian

2.2 Descriptive statistics

Table of baseline characteristics:

> bcva_data |>
+   select(ARMCD, USUBJID, RACE, BCVA_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	CTL, N = 494¹	TRT, N = 506¹
RACE
Asian	151 / 494 (31%)	146 / 506 (29%)
Black	149 / 494 (30%)	168 / 506 (33%)
White	194 / 494 (39%)	192 / 506 (38%)
BCVA_BL	75 (10)	75 (10)
¹ n / N (%); Mean (SD)

Table of change from baseline in BCVA over 52 weeks:

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

Table 3. Change from baseline.
Characteristic	CTL, N = 494	TRT, N = 506
Visit VIS01
Mean (SD)	5.32 (1.23)	5.86 (1.33)
Median (IQR)	5.34 (4.51, 6.17)	5.86 (4.98, 6.81)
Range	1.83, 9.02	2.28, 10.30
Unknown	12	5
Visit VIS02
Mean (SD)	5.59 (1.49)	6.33 (1.45)
Median (IQR)	5.53 (4.64, 6.47)	6.36 (5.34, 7.31)
Range	0.29, 10.15	2.35, 10.75
Unknown	13	7
Visit VIS03
Mean (SD)	5.79 (1.61)	6.79 (1.71)
Median (IQR)	5.73 (4.64, 6.89)	6.82 (5.66, 7.93)
Range	1.53, 11.46	1.13, 11.49
Unknown	23	17
Visit VIS04
Mean (SD)	6.18 (1.73)	7.29 (1.82)
Median (IQR)	6.14 (5.05, 7.40)	7.24 (6.06, 8.54)
Range	0.45, 11.49	2.07, 11.47
Unknown	36	18
Visit VIS05
Mean (SD)	6.28 (1.97)	7.68 (1.94)
Median (IQR)	6.18 (4.96, 7.71)	7.75 (6.48, 8.93)
Range	0.87, 11.53	2.24, 14.10
Unknown	40	35
Visit VIS06
Mean (SD)	6.69 (1.97)	8.31 (1.98)
Median (IQR)	6.64 (5.26, 8.14)	8.29 (6.92, 9.73)
Range	1.35, 12.95	1.93, 14.38
Unknown	84	48
Visit VIS07
Mean (SD)	6.78 (2.09)	8.78 (2.11)
Median (IQR)	6.91 (5.46, 8.29)	8.67 (7.44, 10.25)
Range	-1.54, 11.99	3.21, 14.36
Unknown	106	78
Visit VIS08
Mean (SD)	7.08 (2.25)	9.40 (2.26)
Median (IQR)	7.08 (5.57, 8.67)	9.35 (7.96, 10.86)
Range	0.97, 13.71	2.28, 15.95
Unknown	123	86
Visit VIS09
Mean (SD)	7.39 (2.33)	10.01 (2.50)
Median (IQR)	7.47 (5.77, 9.05)	10.01 (8.19, 11.73)
Range	0.04, 14.61	4.22, 18.09
Unknown	167	114
Visit VIS10
Mean (SD)	7.49 (2.58)	10.59 (2.36)
Median (IQR)	7.40 (5.73, 9.01)	10.71 (9.03, 12.24)
Range	-0.08, 15.75	3.24, 16.40
Unknown	213	170

The following figure shows the primary endpoint over the four study visits in the data.

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

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

3 Fitting MMRMs

3.1 Bayesian model

The formula for the Bayesian model includes additive effects for baseline, study visit, race, and study-arm-by-visit interaction.

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

We fit the model using brms.mmrm::brm_model(). The computation takes several minutes because of the size of the dataset. To ensure a good basis of comparison with the frequentist model, we put an extremely diffuse prior on the intercept. The parameters already have diffuse flexible priors by default.

> b_mmrm_fit <- brm_model(
+   data = filter(bcva_data, !is.na(BCVA_CHG)),
+   formula = b_mmrm_formula,
+   prior = brms::prior(class = "Intercept", prior = "student_t(3, 0, 1000)"),
+   iter = 10000,
+   warmup = 2000,
+   chains = 4,
+   cores = 4,
+   refresh = 0
+ )

Here is a posterior summary of model parameters, including fixed effects and pairwise correlation among visits within patients.

> summary(b_mmrm_fit)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = log 
#> Formula: BCVA_CHG ~ BCVA_BL + ARMCD:AVISIT + AVISIT + RACE + unstr(time = AVISIT, gr = USUBJID) 
#>          sigma ~ 0 + AVISIT
#>    Data: data[!is.na(data[[attr(data, "brm_outcome")]]), ] (Number of observations: 8605) 
#>   Draws: 4 chains, each with iter = 10000; warmup = 2000; thin = 1;
#>          total post-warmup draws = 32000
#> 
#> Correlation Structures:
#>                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> cortime(VIS01,VIS02)     0.80      0.01     0.78     0.82 1.00    25814
#> cortime(VIS01,VIS03)     0.80      0.01     0.77     0.82 1.00    27894
#> cortime(VIS02,VIS03)     0.71      0.02     0.67     0.74 1.00    25786
#> cortime(VIS01,VIS04)     0.70      0.01     0.67     0.73 1.00    21340
#> cortime(VIS02,VIS04)     0.68      0.02     0.64     0.71 1.00    21716
#> cortime(VIS03,VIS04)     0.60      0.02     0.56     0.64 1.00    18655
#> cortime(VIS01,VIS05)     0.61      0.02     0.57     0.64 1.00    22549
#> cortime(VIS02,VIS05)     0.59      0.02     0.55     0.63 1.00    21432
#> cortime(VIS03,VIS05)     0.53      0.03     0.48     0.58 1.00    19140
#> cortime(VIS04,VIS05)     0.71      0.02     0.68     0.74 1.00    25313
#> cortime(VIS01,VIS06)     0.51      0.02     0.47     0.56 1.00    22245
#> cortime(VIS02,VIS06)     0.53      0.02     0.48     0.57 1.00    21179
#> cortime(VIS03,VIS06)     0.48      0.03     0.42     0.53 1.00    18362
#> cortime(VIS04,VIS06)     0.68      0.02     0.64     0.71 1.00    24781
#> cortime(VIS05,VIS06)     0.68      0.02     0.64     0.71 1.00    27072
#> cortime(VIS01,VIS07)     0.26      0.04     0.19     0.33 1.00    29308
#> cortime(VIS02,VIS07)     0.30      0.04     0.23     0.37 1.00    27723
#> cortime(VIS03,VIS07)     0.30      0.04     0.22     0.37 1.00    25913
#> cortime(VIS04,VIS07)     0.41      0.03     0.34     0.48 1.00    28633
#> cortime(VIS05,VIS07)     0.46      0.03     0.40     0.52 1.00    32096
#> cortime(VIS06,VIS07)     0.51      0.03     0.45     0.56 1.00    30156
#> cortime(VIS01,VIS08)     0.09      0.04     0.02     0.16 1.00    36119
#> cortime(VIS02,VIS08)     0.16      0.04     0.08     0.23 1.00    34156
#> cortime(VIS03,VIS08)     0.15      0.04     0.08     0.23 1.00    32502
#> cortime(VIS04,VIS08)     0.31      0.04     0.23     0.38 1.00    34914
#> cortime(VIS05,VIS08)     0.35      0.03     0.28     0.41 1.00    37071
#> cortime(VIS06,VIS08)     0.40      0.03     0.33     0.46 1.00    34382
#> cortime(VIS07,VIS08)     0.36      0.03     0.29     0.42 1.00    38088
#> cortime(VIS01,VIS09)    -0.05      0.03    -0.12     0.01 1.00    46334
#> cortime(VIS02,VIS09)    -0.01      0.04    -0.07     0.06 1.00    42920
#> cortime(VIS03,VIS09)     0.05      0.04    -0.03     0.12 1.00    42833
#> cortime(VIS04,VIS09)     0.16      0.04     0.09     0.23 1.00    44165
#> cortime(VIS05,VIS09)     0.21      0.03     0.14     0.27 1.00    46855
#> cortime(VIS06,VIS09)     0.33      0.03     0.27     0.39 1.00    47695
#> cortime(VIS07,VIS09)     0.29      0.03     0.23     0.36 1.00    49402
#> cortime(VIS08,VIS09)     0.39      0.03     0.33     0.44 1.00    50252
#> cortime(VIS01,VIS10)    -0.10      0.03    -0.16    -0.03 1.00    43885
#> cortime(VIS02,VIS10)    -0.00      0.03    -0.06     0.06 1.00    42961
#> cortime(VIS03,VIS10)     0.04      0.03    -0.03     0.10 1.00    40552
#> cortime(VIS04,VIS10)     0.22      0.03     0.16     0.28 1.00    44812
#> cortime(VIS05,VIS10)     0.27      0.03     0.21     0.33 1.00    49906
#> cortime(VIS06,VIS10)     0.38      0.03     0.32     0.44 1.00    49501
#> cortime(VIS07,VIS10)     0.33      0.03     0.27     0.39 1.00    53935
#> cortime(VIS08,VIS10)     0.41      0.03     0.36     0.47 1.00    53862
#> cortime(VIS09,VIS10)     0.44      0.03     0.39     0.49 1.00    54958
#>                      Tail_ESS
#> cortime(VIS01,VIS02)    24405
#> cortime(VIS01,VIS03)    24039
#> cortime(VIS02,VIS03)    23777
#> cortime(VIS01,VIS04)    24870
#> cortime(VIS02,VIS04)    25165
#> cortime(VIS03,VIS04)    23933
#> cortime(VIS01,VIS05)    24830
#> cortime(VIS02,VIS05)    23805
#> cortime(VIS03,VIS05)    23943
#> cortime(VIS04,VIS05)    25526
#> cortime(VIS01,VIS06)    24409
#> cortime(VIS02,VIS06)    24107
#> cortime(VIS03,VIS06)    24098
#> cortime(VIS04,VIS06)    25983
#> cortime(VIS05,VIS06)    24029
#> cortime(VIS01,VIS07)    24355
#> cortime(VIS02,VIS07)    23774
#> cortime(VIS03,VIS07)    22936
#> cortime(VIS04,VIS07)    23514
#> cortime(VIS05,VIS07)    24212
#> cortime(VIS06,VIS07)    23912
#> cortime(VIS01,VIS08)    25501
#> cortime(VIS02,VIS08)    24810
#> cortime(VIS03,VIS08)    22271
#> cortime(VIS04,VIS08)    22823
#> cortime(VIS05,VIS08)    23153
#> cortime(VIS06,VIS08)    25061
#> cortime(VIS07,VIS08)    24467
#> cortime(VIS01,VIS09)    25154
#> cortime(VIS02,VIS09)    24271
#> cortime(VIS03,VIS09)    24686
#> cortime(VIS04,VIS09)    24026
#> cortime(VIS05,VIS09)    23467
#> cortime(VIS06,VIS09)    25686
#> cortime(VIS07,VIS09)    24723
#> cortime(VIS08,VIS09)    24237
#> cortime(VIS01,VIS10)    24788
#> cortime(VIS02,VIS10)    25198
#> cortime(VIS03,VIS10)    24858
#> cortime(VIS04,VIS10)    23634
#> cortime(VIS05,VIS10)    23341
#> cortime(VIS06,VIS10)    23791
#> cortime(VIS07,VIS10)    24641
#> cortime(VIS08,VIS10)    22307
#> cortime(VIS09,VIS10)    24701
#> 
#> Regression Coefficients:
#>                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept                7.27      0.27     6.74     7.80 1.00    53640
#> BCVA_BL                  0.00      0.00    -0.01     0.01 1.00    84617
#> AVISITVIS02              0.05      0.09    -0.13     0.23 1.00    18165
#> AVISITVIS03              0.12      0.10    -0.07     0.31 1.00    18271
#> AVISITVIS04              0.12      0.10    -0.07     0.32 1.00    18320
#> AVISITVIS05              0.11      0.10    -0.08     0.30 1.00    17714
#> AVISITVIS06              0.19      0.10    -0.01     0.39 1.00    17322
#> AVISITVIS07              0.02      0.11    -0.19     0.23 1.00    18801
#> AVISITVIS08              0.14      0.10    -0.06     0.35 1.00    18460
#> AVISITVIS09              0.25      0.11     0.03     0.47 1.00    21095
#> AVISITVIS10              0.02      0.12    -0.20     0.25 1.00    21735
#> RACEBlack                0.01      0.09    -0.16     0.18 1.00    47505
#> RACEWhite               -0.14      0.08    -0.29     0.02 1.00    47750
#> AVISITVIS01:ARMCDTRT     0.05      0.13    -0.21     0.31 1.00    11737
#> AVISITVIS02:ARMCDTRT    -0.01      0.13    -0.26     0.25 1.00    14570
#> AVISITVIS03:ARMCDTRT    -0.01      0.14    -0.28     0.26 1.00    14593
#> AVISITVIS04:ARMCDTRT    -0.06      0.14    -0.34     0.21 1.00    14824
#> AVISITVIS05:ARMCDTRT    -0.12      0.13    -0.38     0.14 1.00    14287
#> AVISITVIS06:ARMCDTRT    -0.27      0.13    -0.53    -0.00 1.00    14099
#> AVISITVIS07:ARMCDTRT     0.14      0.14    -0.15     0.41 1.00    16235
#> AVISITVIS08:ARMCDTRT    -0.04      0.14    -0.32     0.24 1.00    14965
#> AVISITVIS09:ARMCDTRT     0.01      0.16    -0.30     0.32 1.00    16977
#> AVISITVIS10:ARMCDTRT     0.12      0.16    -0.20     0.43 1.00    17595
#> sigma_AVISITVIS01        0.92      0.02     0.88     0.96 1.00    23844
#> sigma_AVISITVIS02        0.92      0.02     0.87     0.96 1.00    26609
#> sigma_AVISITVIS03        0.96      0.02     0.92     1.01 1.00    27236
#> sigma_AVISITVIS04        0.96      0.02     0.92     1.00 1.00    27047
#> sigma_AVISITVIS05        0.91      0.02     0.87     0.95 1.00    27612
#> sigma_AVISITVIS06        0.91      0.02     0.86     0.95 1.00    26448
#> sigma_AVISITVIS07        0.94      0.02     0.90     0.99 1.00    28187
#> sigma_AVISITVIS08        0.92      0.02     0.88     0.97 1.00    27560
#> sigma_AVISITVIS09        0.98      0.02     0.93     1.02 1.00    29896
#> sigma_AVISITVIS10        0.96      0.03     0.91     1.01 1.00    31353
#>                      Tail_ESS
#> Intercept               22503
#> BCVA_BL                 22111
#> AVISITVIS02             21902
#> AVISITVIS03             22699
#> AVISITVIS04             23308
#> AVISITVIS05             21610
#> AVISITVIS06             21869
#> AVISITVIS07             22286
#> AVISITVIS08             23739
#> AVISITVIS09             23571
#> AVISITVIS10             23915
#> RACEBlack               25965
#> RACEWhite               27488
#> AVISITVIS01:ARMCDTRT    19550
#> AVISITVIS02:ARMCDTRT    20053
#> AVISITVIS03:ARMCDTRT    20491
#> AVISITVIS04:ARMCDTRT    21692
#> AVISITVIS05:ARMCDTRT    21524
#> AVISITVIS06:ARMCDTRT    21221
#> AVISITVIS07:ARMCDTRT    21230
#> AVISITVIS08:ARMCDTRT    20811
#> AVISITVIS09:ARMCDTRT    21823
#> AVISITVIS10:ARMCDTRT    22135
#> sigma_AVISITVIS01       23146
#> sigma_AVISITVIS02       23958
#> sigma_AVISITVIS03       23743
#> sigma_AVISITVIS04       23701
#> sigma_AVISITVIS05       22645
#> sigma_AVISITVIS06       23962
#> sigma_AVISITVIS07       23653
#> sigma_AVISITVIS08       23029
#> sigma_AVISITVIS09       22924
#> sigma_AVISITVIS10       23092
#> 
#> 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).

3.2 Frequentist model

The formula for the frequentist model is the same, except for the different syntax for specifying the covariance structure of the MMRM. We fit the model below.

> f_mmrm_fit <- mmrm::mmrm(
+   formula = BCVA_CHG ~ BCVA_BL + ARMCD:AVISIT + AVISIT + RACE +
+     us(AVISIT | USUBJID),
+   data = bcva_data
+ )

The parameter summaries of the frequentist model are below.

> summary(f_mmrm_fit)
#> mmrm fit
#> 
#> Formula:     BCVA_CHG ~ BCVA_BL + ARMCD:AVISIT + AVISIT + RACE + us(AVISIT |  
#>     USUBJID)
#> Data:        bcva_data (used 8605 observations from 1000 subjects with maximum 
#> 10 timepoints)
#> Covariance:  unstructured (55 variance parameters)
#> Method:      Satterthwaite
#> Vcov Method: Asymptotic
#> Inference:   REML
#> 
#> Model selection criteria:
#>      AIC      BIC   logLik deviance 
#>  32181.0  32451.0 -16035.5  32071.0 
#> 
#> Coefficients: 
#>                        Estimate Std. Error         df t value Pr(>|t|)    
#> (Intercept)           4.288e+00  1.709e-01  1.065e+03  25.086  < 2e-16 ***
#> BCVA_BL              -9.933e-04  2.156e-03  9.906e+02  -0.461    0.645    
#> AVISITVIS02           2.810e-01  7.067e-02  9.995e+02   3.976 7.51e-05 ***
#> AVISITVIS03           4.573e-01  6.716e-02  9.747e+02   6.809 1.71e-11 ***
#> AVISITVIS04           8.570e-01  7.637e-02  9.795e+02  11.221  < 2e-16 ***
#> AVISITVIS05           9.638e-01  8.634e-02  9.629e+02  11.163  < 2e-16 ***
#> AVISITVIS06           1.334e+00  8.650e-02  9.450e+02  15.421  < 2e-16 ***
#> AVISITVIS07           1.417e+00  1.071e-01  8.698e+02  13.233  < 2e-16 ***
#> AVISITVIS08           1.711e+00  1.145e-01  8.467e+02  14.943  < 2e-16 ***
#> AVISITVIS09           1.996e+00  1.283e-01  7.784e+02  15.550  < 2e-16 ***
#> AVISITVIS10           2.101e+00  1.400e-01  7.025e+02  15.003  < 2e-16 ***
#> RACEBlack             1.038e+00  5.496e-02  1.011e+03  18.892  < 2e-16 ***
#> RACEWhite             2.005e+00  5.198e-02  9.769e+02  38.574  < 2e-16 ***
#> AVISITVIS01:ARMCDTRT  5.391e-01  6.281e-02  9.859e+02   8.583  < 2e-16 ***
#> AVISITVIS02:ARMCDTRT  7.248e-01  7.984e-02  9.803e+02   9.078  < 2e-16 ***
#> AVISITVIS03:ARMCDTRT  1.012e+00  9.163e-02  9.638e+02  11.039  < 2e-16 ***
#> AVISITVIS04:ARMCDTRT  1.104e+00  1.004e-01  9.653e+02  11.003  < 2e-16 ***
#> AVISITVIS05:ARMCDTRT  1.383e+00  1.147e-01  9.505e+02  12.065  < 2e-16 ***
#> AVISITVIS06:ARMCDTRT  1.630e+00  1.189e-01  9.157e+02  13.715  < 2e-16 ***
#> AVISITVIS07:ARMCDTRT  2.016e+00  1.382e-01  8.262e+02  14.592  < 2e-16 ***
#> AVISITVIS08:ARMCDTRT  2.347e+00  1.474e-01  8.041e+02  15.924  < 2e-16 ***
#> AVISITVIS09:ARMCDTRT  2.658e+00  1.644e-01  7.277e+02  16.173  < 2e-16 ***
#> AVISITVIS10:ARMCDTRT  3.072e+00  1.815e-01  6.621e+02  16.929  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Covariance estimate:
#>        VIS01   VIS02   VIS03   VIS04   VIS05  VIS06   VIS07   VIS08   VIS09
#> VIS01 0.9712  0.0630  0.4371  0.3314  0.3055 0.4686  0.1324  0.1019  0.0610
#> VIS02 0.0630  1.5618  0.0871  0.2685  0.2635 0.4636  0.2180  0.2776 -0.0153
#> VIS03 0.4371  0.0871  2.0221 -0.0216 -0.0189 0.1102 -0.0048 -0.0993 -0.1322
#> VIS04 0.3314  0.2685 -0.0216  2.4114  1.0476 1.1409  0.4625  0.5660  0.4086
#> VIS05 0.3055  0.2635 -0.0189  1.0476  3.0915 1.2592  0.6909  0.6307  0.3593
#> VIS06 0.4686  0.4636  0.1102  1.1409  1.2592 3.1852  0.7539  0.6093  0.6821
#> VIS07 0.1324  0.2180 -0.0048  0.4625  0.6909 0.7539  3.9273  0.2306  0.0723
#> VIS08 0.1019  0.2776 -0.0993  0.5660  0.6307 0.6093  0.2306  4.3272  0.2682
#> VIS09 0.0610 -0.0153 -0.1322  0.4086  0.3593 0.6821  0.0723  0.2682  4.8635
#> VIS10 0.0585  0.3762  0.0719  1.1481  0.9999 1.2559  0.3017  0.4658  0.4138
#>        VIS10
#> VIS01 0.0585
#> VIS02 0.3762
#> VIS03 0.0719
#> VIS04 1.1481
#> VIS05 0.9999
#> VIS06 1.2559
#> VIS07 0.3017
#> VIS08 0.4658
#> VIS09 0.4138
#> VIS10 5.3520

4 Comparison

This section compares the Bayesian posterior parameter estimates from brms.mmrm to the frequentist parameter estimates of the mmrm package.

4.1 Extract estimates from Bayesian model

We extract and standardize the Bayesian estimates.

> b_mmrm_draws <- b_mmrm_fit |>
+   as_draws_df()
> visit_levels <- sort(unique(as.character(bcva_data$AVISIT)))
> for (level in visit_levels) {
+   name <- paste0("b_sigma_AVISIT", level)
+   b_mmrm_draws[[name]] <- exp(b_mmrm_draws[[name]])
+ }
> b_mmrm_summary <- b_mmrm_draws |>
+   summarize_draws() |>
+   select(variable, mean, sd) |>
+   filter(!(variable %in% c("lprior", "lp__"))) |>
+   rename(bayes_estimate = mean, bayes_se = sd) |>
+   mutate(
+     variable = variable |>
+       tolower() |>
+       gsub(pattern = "b_", replacement = "") |>
+       gsub(pattern = "b_sigma_AVISIT", replacement = "sigma_") |>
+       gsub(pattern = "cortime", replacement = "correlation") |>
+       gsub(pattern = "__", replacement = "_")
+   )

4.2 Extract estimates from frequentist model

We extract and standardize the frequentist estimates.

> f_mmrm_fixed <- summary(f_mmrm_fit)$coefficients |>
+   as_tibble(rownames = "variable") |>
+   mutate(variable = tolower(variable)) |>
+   mutate(variable = gsub("(", "", variable, fixed = TRUE)) |>
+   mutate(variable = gsub(")", "", variable, fixed = TRUE)) |>
+   rename(freq_estimate = Estimate, freq_se = `Std. Error`) |>
+   select(variable, freq_estimate, freq_se)

> f_mmrm_variance <- tibble(
+   variable = paste0("sigma_AVISIT", visit_levels) |> tolower(),
+   freq_estimate = sqrt(diag(f_mmrm_fit$cov))
+ )

> f_diagonal_factor <- diag(1 / sqrt(diag(f_mmrm_fit$cov)))
> f_corr_matrix <- f_diagonal_factor %*% f_mmrm_fit$cov %*% f_diagonal_factor
> colnames(f_corr_matrix) <- visit_levels

> f_mmrm_correlation <- f_corr_matrix |>
+   as.data.frame() |>
+   as_tibble() |>
+   mutate(x1 = visit_levels) |>
+   pivot_longer(
+     cols = -any_of("x1"),
+     names_to = "x2",
+     values_to = "freq_estimate"
+   ) |>
+   filter(
+     as.numeric(gsub("[^0-9]", "", x1)) < as.numeric(gsub("[^0-9]", "", x2))
+   ) |>
+   mutate(variable = sprintf("correlation_%s_%s", x1, x2)) |>
+   select(variable, freq_estimate)

> f_mmrm_summary <- bind_rows(
+   f_mmrm_fixed,
+   f_mmrm_variance,
+   f_mmrm_correlation
+ ) |>
+   mutate(variable = gsub("\\s+", "", variable) |> tolower())

4.3 Summary

The first table below summarizes the parameter estimates from each model and the differences between estimates (Bayesian minus frequentist). The second table shows the standard errors of these estimates and differences between standard errors. In each table, the “Relative” column shows the relative difference (the difference divided by the frequentist quantity).

Because of the different statistical paradigms and estimation procedures, especially regarding the covariance parameters, it would not be realistic to expect the Bayesian and frequentist approaches to yield virtually identical results. Nevertheless, the absolute and relative differences in the table below show strong agreement between brms.mmrm and mmrm.

> b_f_comparison <- full_join(
+   x = b_mmrm_summary,
+   y = f_mmrm_summary,
+   by = "variable"
+ ) |>
+   mutate(
+     diff_estimate = bayes_estimate - freq_estimate,
+     diff_relative_estimate = diff_estimate / freq_estimate,
+     diff_se = bayes_se - freq_se,
+     diff_relative_se = diff_se / freq_se
+   ) |>
+   select(variable, ends_with("estimate"), ends_with("se"))

> table_estimates <- b_f_comparison |>
+   select(variable, ends_with("estimate"))
> gt(table_estimates) |>
+   fmt_number(decimals = 4) |>
+   tab_caption(
+     caption = md(
+       paste(
+         "Table 4. Comparison of parameter estimates between",
+         "Bayesian and frequentist MMRMs."
+       )
+     )
+   ) |>
+   cols_label(
+     variable = "Variable",
+     bayes_estimate = "Bayesian",
+     freq_estimate = "Frequentist",
+     diff_estimate = "Difference",
+     diff_relative_estimate = "Relative"
+   )

Table 4. Comparison of parameter estimates between Bayesian and frequentist MMRMs.
Variable	Bayesian	Frequentist	Difference	Relative
intercept	7.2672	4.2881	2.9792	0.6948
bcva_bl	0.0013	−0.0010	0.0023	−2.2679
avisitvis02	0.0518	0.2810	−0.2292	−0.8155
avisitvis03	0.1197	0.4573	−0.3376	−0.7383
avisitvis04	0.1205	0.8570	−0.7364	−0.8593
avisitvis05	0.1105	0.9638	−0.8533	−0.8854
avisitvis06	0.1898	1.3339	−1.1442	−0.8577
avisitvis07	0.0157	1.4167	−1.4010	−0.9889
avisitvis08	0.1418	1.7107	−1.5690	−0.9171
avisitvis09	0.2534	1.9956	−1.7422	−0.8730
avisitvis10	0.0199	2.1005	−2.0806	−0.9905
raceblack	0.0063	1.0382	−1.0319	−0.9939
racewhite	−0.1353	2.0051	−2.1404	−1.0675
avisitvis01:armcdtrt	0.0504	0.5391	−0.4887	−0.9065
avisitvis02:armcdtrt	−0.0055	0.7248	−0.7303	−1.0076
avisitvis03:armcdtrt	−0.0110	1.0115	−1.0225	−1.0109
avisitvis04:armcdtrt	−0.0607	1.1042	−1.1649	−1.0550
avisitvis05:armcdtrt	−0.1227	1.3834	−1.5061	−1.0887
avisitvis06:armcdtrt	−0.2682	1.6301	−1.8983	−1.1646
avisitvis07:armcdtrt	0.1358	2.0160	−1.8802	−0.9327
avisitvis08:armcdtrt	−0.0393	2.3469	−2.3862	−1.0167
avisitvis09:armcdtrt	0.0090	2.6585	−2.6494	−0.9966
avisitvis10:armcdtrt	0.1170	3.0722	−2.9552	−0.9619
sigma_avisitvis01	2.5121	0.9855	1.5266	1.5490
sigma_avisitvis02	2.4975	1.2497	1.2478	0.9984
sigma_avisitvis03	2.6200	1.4220	1.1980	0.8424
sigma_avisitvis04	2.6149	1.5529	1.0620	0.6839
sigma_avisitvis05	2.4814	1.7583	0.7231	0.4113
sigma_avisitvis06	2.4810	1.7847	0.6963	0.3901
sigma_avisitvis07	2.5621	1.9817	0.5804	0.2929
sigma_avisitvis08	2.5214	2.0802	0.4412	0.2121
sigma_avisitvis09	2.6557	2.2053	0.4504	0.2042
sigma_avisitvis10	2.6160	2.3134	0.3026	0.1308
correlation_vis01_vis02	0.8009	0.0511	0.7498	14.6643
correlation_vis01_vis03	0.7952	0.3119	0.4833	1.5494
correlation_vis02_vis03	0.7067	0.0490	0.6576	13.4167
correlation_vis01_vis04	0.7017	0.2165	0.4851	2.2405
correlation_vis02_vis04	0.6771	0.1383	0.5388	3.8945
correlation_vis03_vis04	0.6030	−0.0098	0.6128	−62.6953
correlation_vis01_vis05	0.6062	0.1763	0.4299	2.4384
correlation_vis02_vis05	0.5925	0.1199	0.4726	3.9407
correlation_vis03_vis05	0.5323	−0.0076	0.5399	−71.2641
correlation_vis04_vis05	0.7124	0.3837	0.3288	0.8569
correlation_vis01_vis06	0.5140	0.2665	0.2475	0.9290
correlation_vis02_vis06	0.5252	0.2079	0.3173	1.5264
correlation_vis03_vis06	0.4765	0.0434	0.4330	9.9729
correlation_vis04_vis06	0.6757	0.4117	0.2640	0.6413
correlation_vis05_vis06	0.6757	0.4013	0.2745	0.6840
correlation_vis01_vis07	0.2628	0.0678	0.1950	2.8751
correlation_vis02_vis07	0.3037	0.0880	0.2156	2.4490
correlation_vis03_vis07	0.2995	−0.0017	0.3012	−176.4416
correlation_vis04_vis07	0.4146	0.1503	0.2643	1.7587
correlation_vis05_vis07	0.4593	0.1983	0.2610	1.3163
correlation_vis06_vis07	0.5080	0.2132	0.2949	1.3833
correlation_vis01_vis08	0.0934	0.0497	0.0437	0.8784
correlation_vis02_vis08	0.1580	0.1068	0.0512	0.4791
correlation_vis03_vis08	0.1545	−0.0336	0.1881	−5.6048
correlation_vis04_vis08	0.3076	0.1752	0.1324	0.7558
correlation_vis05_vis08	0.3451	0.1724	0.1726	1.0011
correlation_vis06_vis08	0.3953	0.1641	0.2312	1.4088
correlation_vis07_vis08	0.3600	0.0559	0.3040	5.4355
correlation_vis01_vis09	−0.0549	0.0281	−0.0830	−2.9540
correlation_vis02_vis09	−0.0054	−0.0056	0.0002	−0.0276
correlation_vis03_vis09	0.0457	−0.0421	0.0879	−2.0854
correlation_vis04_vis09	0.1609	0.1193	0.0416	0.3486
correlation_vis05_vis09	0.2084	0.0927	0.1158	1.2494
correlation_vis06_vis09	0.3311	0.1733	0.1578	0.9106
correlation_vis07_vis09	0.2949	0.0165	0.2784	16.8264
correlation_vis08_vis09	0.3897	0.0585	0.3312	5.6660
correlation_vis01_vis10	−0.0950	0.0257	−0.1207	−4.7031
correlation_vis02_vis10	−0.0024	0.1301	−0.1325	−1.0181
correlation_vis03_vis10	0.0377	0.0219	0.0159	0.7261
correlation_vis04_vis10	0.2238	0.3196	−0.0958	−0.2998
correlation_vis05_vis10	0.2690	0.2458	0.0231	0.0942
correlation_vis06_vis10	0.3799	0.3042	0.0757	0.2490
correlation_vis07_vis10	0.3315	0.0658	0.2657	4.0375
correlation_vis08_vis10	0.4132	0.0968	0.3164	3.2689
correlation_vis09_vis10	0.4449	0.0811	0.3638	4.4851
intercept	7.3964	4.2881	3.1083	0.7249

> table_se <- b_f_comparison |>
+   select(variable, ends_with("se")) |>
+   filter(!is.na(freq_se))
> gt(table_se) |>
+   fmt_number(decimals = 4) |>
+   tab_caption(
+     caption = md(
+       paste(
+         "Table 5. Comparison of parameter standard errors between",
+         "Bayesian and frequentist MMRMs."
+       )
+     )
+   ) |>
+   cols_label(
+     variable = "Variable",
+     bayes_se = "Bayesian",
+     freq_se = "Frequentist",
+     diff_se = "Difference",
+     diff_relative_se = "Relative"
+   )

Table 5. Comparison of parameter standard errors between Bayesian and frequentist MMRMs.
Variable	Bayesian	Frequentist	Difference	Relative
intercept	0.2708	0.1709	0.0998	0.5840
bcva_bl	0.0033	0.0022	0.0012	0.5468
avisitvis02	0.0921	0.0707	0.0214	0.3033
avisitvis03	0.0955	0.0672	0.0284	0.4226
avisitvis04	0.1007	0.0764	0.0243	0.3188
avisitvis05	0.0976	0.0863	0.0113	0.1307
avisitvis06	0.0996	0.0865	0.0131	0.1512
avisitvis07	0.1059	0.1071	−0.0012	−0.0110
avisitvis08	0.1039	0.1145	−0.0106	−0.0925
avisitvis09	0.1108	0.1283	−0.0175	−0.1365
avisitvis10	0.1153	0.1400	−0.0247	−0.1764
raceblack	0.0859	0.0550	0.0310	0.5632
racewhite	0.0817	0.0520	0.0298	0.5724
avisitvis01:armcdtrt	0.1320	0.0628	0.0692	1.1013
avisitvis02:armcdtrt	0.1318	0.0798	0.0520	0.6511
avisitvis03:armcdtrt	0.1384	0.0916	0.0468	0.5108
avisitvis04:armcdtrt	0.1411	0.1004	0.0407	0.4058
avisitvis05:armcdtrt	0.1336	0.1147	0.0189	0.1648
avisitvis06:armcdtrt	0.1341	0.1189	0.0152	0.1280
avisitvis07:armcdtrt	0.1437	0.1382	0.0055	0.0402
avisitvis08:armcdtrt	0.1427	0.1474	−0.0046	−0.0315
avisitvis09:armcdtrt	0.1560	0.1644	−0.0084	−0.0511
avisitvis10:armcdtrt	0.1609	0.1815	−0.0206	−0.1135
intercept	0.0533	0.1709	−0.1176	−0.6883

BCVA data comparison between Bayesian and frequentist MMRMs

31/05/2024