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 Tail_ESS
#> cortime(VIS01,VIS02)     0.05      0.03    -0.01     0.11 1.00    64943    24580
#> cortime(VIS01,VIS03)     0.31      0.03     0.25     0.37 1.00    64795    24810
#> cortime(VIS02,VIS03)     0.05      0.03    -0.01     0.11 1.00    75095    23845
#> cortime(VIS01,VIS04)     0.21      0.03     0.15     0.27 1.00    48997    25597
#> cortime(VIS02,VIS04)     0.13      0.03     0.07     0.20 1.00    53822    27463
#> cortime(VIS03,VIS04)    -0.01      0.03    -0.07     0.05 1.00    49973    25432
#> cortime(VIS01,VIS05)     0.17      0.03     0.11     0.23 1.00    47899    25766
#> cortime(VIS02,VIS05)     0.12      0.03     0.05     0.18 1.00    55002    27003
#> cortime(VIS03,VIS05)    -0.01      0.03    -0.07     0.06 1.00    50643    26319
#> cortime(VIS04,VIS05)     0.38      0.03     0.32     0.43 1.00    55423    27026
#> cortime(VIS01,VIS06)     0.26      0.03     0.20     0.32 1.00    44737    28050
#> cortime(VIS02,VIS06)     0.20      0.03     0.14     0.27 1.00    54726    27792
#> cortime(VIS03,VIS06)     0.04      0.03    -0.02     0.11 1.00    52319    27395
#> cortime(VIS04,VIS06)     0.40      0.03     0.35     0.46 1.00    51509    26799
#> cortime(VIS05,VIS06)     0.39      0.03     0.34     0.45 1.00    55687    26874
#> cortime(VIS01,VIS07)     0.07      0.04    -0.00     0.13 1.00    67222    23189
#> cortime(VIS02,VIS07)     0.09      0.03     0.02     0.15 1.00    68013    23631
#> cortime(VIS03,VIS07)    -0.00      0.03    -0.07     0.07 1.00    67332    24340
#> cortime(VIS04,VIS07)     0.15      0.03     0.08     0.21 1.00    66288    24838
#> cortime(VIS05,VIS07)     0.19      0.03     0.13     0.26 1.00    70461    24208
#> cortime(VIS06,VIS07)     0.21      0.04     0.14     0.28 1.00    69109    24028
#> cortime(VIS01,VIS08)     0.05      0.04    -0.02     0.12 1.00    72277    22000
#> cortime(VIS02,VIS08)     0.10      0.04     0.03     0.17 1.00    73280    24125
#> cortime(VIS03,VIS08)    -0.03      0.04    -0.10     0.04 1.00    68920    24860
#> cortime(VIS04,VIS08)     0.17      0.03     0.10     0.24 1.00    67547    23658
#> cortime(VIS05,VIS08)     0.17      0.04     0.10     0.24 1.00    71694    24584
#> cortime(VIS06,VIS08)     0.16      0.04     0.09     0.23 1.00    69518    24126
#> cortime(VIS07,VIS08)     0.05      0.04    -0.02     0.13 1.00    62240    22292
#> cortime(VIS01,VIS09)     0.03      0.04    -0.04     0.10 1.00    68423    23721
#> cortime(VIS02,VIS09)    -0.01      0.04    -0.08     0.07 1.00    74240    21638
#> cortime(VIS03,VIS09)    -0.04      0.04    -0.12     0.03 1.00    70100    22754
#> cortime(VIS04,VIS09)     0.12      0.04     0.04     0.19 1.00    69789    22965
#> cortime(VIS05,VIS09)     0.09      0.04     0.02     0.16 1.00    72416    24164
#> cortime(VIS06,VIS09)     0.17      0.04     0.09     0.24 1.00    73834    24991
#> cortime(VIS07,VIS09)     0.02      0.04    -0.06     0.09 1.00    70843    23470
#> cortime(VIS08,VIS09)     0.06      0.04    -0.02     0.14 1.00    68629    22953
#> cortime(VIS01,VIS10)     0.02      0.04    -0.06     0.10 1.00    61683    25514
#> cortime(VIS02,VIS10)     0.13      0.04     0.05     0.20 1.00    62160    26376
#> cortime(VIS03,VIS10)     0.02      0.04    -0.06     0.10 1.00    62572    25325
#> cortime(VIS04,VIS10)     0.31      0.04     0.24     0.38 1.00    62792    23595
#> cortime(VIS05,VIS10)     0.24      0.04     0.16     0.31 1.00    64921    24888
#> cortime(VIS06,VIS10)     0.30      0.04     0.22     0.37 1.00    64007    23874
#> cortime(VIS07,VIS10)     0.06      0.05    -0.03     0.15 1.00    68579    23223
#> cortime(VIS08,VIS10)     0.09      0.04     0.01     0.18 1.00    67211    24274
#> cortime(VIS09,VIS10)     0.08      0.05    -0.01     0.17 1.00    68075    23303
#> 
#> Population-Level Effects: 
#>                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept                4.29      0.17     3.95     4.62 1.00    54732    23560
#> BCVA_BL                 -0.00      0.00    -0.01     0.00 1.00    57530    22363
#> AVISITVIS02              0.28      0.07     0.14     0.42 1.00    31423    27345
#> AVISITVIS03              0.46      0.07     0.33     0.59 1.00    43974    27529
#> AVISITVIS04              0.86      0.08     0.71     1.01 1.00    28624    26174
#> AVISITVIS05              0.96      0.09     0.79     1.14 1.00    28737    25221
#> AVISITVIS06              1.33      0.09     1.17     1.50 1.00    28571    26120
#> AVISITVIS07              1.42      0.11     1.20     1.63 1.00    35204    26051
#> AVISITVIS08              1.71      0.12     1.48     1.94 1.00    34176    25768
#> AVISITVIS09              2.00      0.13     1.74     2.25 1.00    35133    25230
#> AVISITVIS10              2.10      0.14     1.82     2.38 1.00    32940    27567
#> RACEBlack                1.04      0.05     0.93     1.14 1.00    50156    25546
#> RACEWhite                2.01      0.05     1.90     2.11 1.00    53389    27446
#> AVISITVIS01:ARMCDTRT     0.54      0.06     0.42     0.67 1.00    34020    27641
#> AVISITVIS02:ARMCDTRT     0.72      0.08     0.57     0.88 1.00    50604    24914
#> AVISITVIS03:ARMCDTRT     1.01      0.09     0.83     1.19 1.00    47195    27254
#> AVISITVIS04:ARMCDTRT     1.10      0.10     0.91     1.30 1.00    37121    27502
#> AVISITVIS05:ARMCDTRT     1.38      0.12     1.15     1.61 1.00    37968    27516
#> AVISITVIS06:ARMCDTRT     1.63      0.12     1.40     1.86 1.00    34971    26561
#> AVISITVIS07:ARMCDTRT     2.02      0.14     1.75     2.29 1.00    45423    25443
#> AVISITVIS08:ARMCDTRT     2.35      0.15     2.05     2.64 1.00    42623    24864
#> AVISITVIS09:ARMCDTRT     2.66      0.17     2.33     2.99 1.00    42236    26456
#> AVISITVIS10:ARMCDTRT     3.07      0.18     2.72     3.43 1.00    38944    24776
#> sigma_AVISITVIS01       -0.01      0.02    -0.05     0.03 1.00    66786    24815
#> sigma_AVISITVIS02        0.23      0.02     0.18     0.27 1.00    73930    23895
#> sigma_AVISITVIS03        0.36      0.02     0.31     0.40 1.00    71987    25366
#> sigma_AVISITVIS04        0.44      0.02     0.40     0.49 1.00    59027    25866
#> sigma_AVISITVIS05        0.57      0.02     0.52     0.61 1.00    61245    24199
#> sigma_AVISITVIS06        0.58      0.02     0.53     0.63 1.00    56690    26883
#> sigma_AVISITVIS07        0.69      0.03     0.64     0.74 1.00    77367    22628
#> sigma_AVISITVIS08        0.74      0.03     0.69     0.79 1.00    71907    22165
#> sigma_AVISITVIS09        0.80      0.03     0.75     0.85 1.00    70162    24546
#> sigma_AVISITVIS10        0.84      0.03     0.79     0.90 1.00    64936    25928
#> 
#> 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  VIS10
#> VIS01 0.9712  0.0630  0.4371  0.3314  0.3055 0.4686  0.1324  0.1019  0.0610 0.0585
#> VIS02 0.0630  1.5618  0.0871  0.2685  0.2635 0.4636  0.2180  0.2776 -0.0153 0.3762
#> VIS03 0.4371  0.0871  2.0221 -0.0216 -0.0189 0.1102 -0.0048 -0.0993 -0.1322 0.0719
#> VIS04 0.3314  0.2685 -0.0216  2.4114  1.0476 1.1409  0.4625  0.5660  0.4086 1.1481
#> VIS05 0.3055  0.2635 -0.0189  1.0476  3.0915 1.2592  0.6909  0.6307  0.3593 0.9999
#> VIS06 0.4686  0.4636  0.1102  1.1409  1.2592 3.1852  0.7539  0.6093  0.6821 1.2559
#> VIS07 0.1324  0.2180 -0.0048  0.4625  0.6909 0.7539  3.9273  0.2306  0.0723 0.3017
#> VIS08 0.1019  0.2776 -0.0993  0.5660  0.6307 0.6093  0.2306  4.3272  0.2682 0.4658
#> VIS09 0.0610 -0.0153 -0.1322  0.4086  0.3593 0.6821  0.0723  0.2682  4.8635 0.4138
#> VIS10 0.0585  0.3762  0.0719  1.1481  0.9999 1.2559  0.3017  0.4658  0.4138 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	4.2879	4.2881	−0.0001	0.0000
bcva_bl	−0.0010	−0.0010	0.0000	0.0083
avisitvis02	0.2816	0.2810	0.0006	0.0021
avisitvis03	0.4579	0.4573	0.0006	0.0013
avisitvis04	0.8570	0.8570	0.0000	0.0001
avisitvis05	0.9646	0.9638	0.0008	0.0008
avisitvis06	1.3346	1.3339	0.0007	0.0005
avisitvis07	1.4175	1.4167	0.0008	0.0005
avisitvis08	1.7110	1.7107	0.0002	0.0001
avisitvis09	1.9959	1.9956	0.0003	0.0002
avisitvis10	2.1015	2.1005	0.0010	0.0005
raceblack	1.0386	1.0382	0.0003	0.0003
racewhite	2.0056	2.0051	0.0006	0.0003
avisitvis01:armcdtrt	0.5393	0.5391	0.0002	0.0003
avisitvis02:armcdtrt	0.7243	0.7248	−0.0005	−0.0006
avisitvis03:armcdtrt	1.0108	1.0115	−0.0007	−0.0007
avisitvis04:armcdtrt	1.1047	1.1042	0.0005	0.0005
avisitvis05:armcdtrt	1.3833	1.3834	−0.0001	−0.0001
avisitvis06:armcdtrt	1.6293	1.6301	−0.0008	−0.0005
avisitvis07:armcdtrt	2.0156	2.0160	−0.0003	−0.0002
avisitvis08:armcdtrt	2.3476	2.3469	0.0007	0.0003
avisitvis09:armcdtrt	2.6579	2.6585	−0.0005	−0.0002
avisitvis10:armcdtrt	3.0717	3.0722	−0.0005	−0.0002
sigma_avisitvis01	0.9895	0.9855	0.0040	0.0040
sigma_avisitvis02	1.2555	1.2497	0.0058	0.0047
sigma_avisitvis03	1.4288	1.4220	0.0068	0.0048
sigma_avisitvis04	1.5563	1.5529	0.0035	0.0022
sigma_avisitvis05	1.7631	1.7583	0.0048	0.0027
sigma_avisitvis06	1.7882	1.7847	0.0035	0.0019
sigma_avisitvis07	1.9928	1.9817	0.0111	0.0056
sigma_avisitvis08	2.0926	2.0802	0.0124	0.0059
sigma_avisitvis09	2.2204	2.2053	0.0150	0.0068
sigma_avisitvis10	2.3276	2.3134	0.0142	0.0061
correlation_vis01_vis02	0.0490	0.0511	−0.0021	−0.0420
correlation_vis01_vis03	0.3086	0.3119	−0.0033	−0.0106
correlation_vis02_vis03	0.0483	0.0490	−0.0008	−0.0154
correlation_vis01_vis04	0.2125	0.2165	−0.0040	−0.0184
correlation_vis02_vis04	0.1349	0.1383	−0.0034	−0.0248
correlation_vis03_vis04	−0.0107	−0.0098	−0.0009	0.0910
correlation_vis01_vis05	0.1723	0.1763	−0.0041	−0.0230
correlation_vis02_vis05	0.1165	0.1199	−0.0034	−0.0283
correlation_vis03_vis05	−0.0082	−0.0076	−0.0006	0.0833
correlation_vis04_vis05	0.3764	0.3837	−0.0073	−0.0190
correlation_vis01_vis06	0.2617	0.2665	−0.0047	−0.0177
correlation_vis02_vis06	0.2038	0.2079	−0.0041	−0.0198
correlation_vis03_vis06	0.0423	0.0434	−0.0011	−0.0260
correlation_vis04_vis06	0.4041	0.4117	−0.0076	−0.0184
correlation_vis05_vis06	0.3938	0.4013	−0.0075	−0.0186
correlation_vis01_vis07	0.0656	0.0678	−0.0022	−0.0326
correlation_vis02_vis07	0.0859	0.0880	−0.0021	−0.0240
correlation_vis03_vis07	−0.0021	−0.0017	−0.0003	0.2020
correlation_vis04_vis07	0.1459	0.1503	−0.0044	−0.0292
correlation_vis05_vis07	0.1936	0.1983	−0.0046	−0.0234
correlation_vis06_vis07	0.2080	0.2132	−0.0052	−0.0242
correlation_vis01_vis08	0.0476	0.0497	−0.0021	−0.0413
correlation_vis02_vis08	0.1041	0.1068	−0.0027	−0.0249
correlation_vis03_vis08	−0.0332	−0.0336	0.0003	−0.0096
correlation_vis04_vis08	0.1707	0.1752	−0.0045	−0.0255
correlation_vis05_vis08	0.1682	0.1724	−0.0043	−0.0248
correlation_vis06_vis08	0.1595	0.1641	−0.0046	−0.0282
correlation_vis07_vis08	0.0538	0.0559	−0.0022	−0.0385
correlation_vis01_vis09	0.0269	0.0281	−0.0012	−0.0426
correlation_vis02_vis09	−0.0065	−0.0056	−0.0010	0.1763
correlation_vis03_vis09	−0.0418	−0.0421	0.0004	−0.0084
correlation_vis04_vis09	0.1160	0.1193	−0.0033	−0.0280
correlation_vis05_vis09	0.0894	0.0927	−0.0033	−0.0351
correlation_vis06_vis09	0.1692	0.1733	−0.0041	−0.0235
correlation_vis07_vis09	0.0153	0.0165	−0.0012	−0.0748
correlation_vis08_vis09	0.0565	0.0585	−0.0020	−0.0337
correlation_vis01_vis10	0.0230	0.0257	−0.0027	−0.1053
correlation_vis02_vis10	0.1264	0.1301	−0.0037	−0.0287
correlation_vis03_vis10	0.0215	0.0219	−0.0003	−0.0142
correlation_vis04_vis10	0.3115	0.3196	−0.0081	−0.0254
correlation_vis05_vis10	0.2384	0.2458	−0.0075	−0.0303
correlation_vis06_vis10	0.2957	0.3042	−0.0085	−0.0279
correlation_vis07_vis10	0.0627	0.0658	−0.0031	−0.0478
correlation_vis08_vis10	0.0931	0.0968	−0.0037	−0.0380
correlation_vis09_vis10	0.0780	0.0811	−0.0031	−0.0380

> 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.1710	0.1709	0.0001	0.0005
bcva_bl	0.0022	0.0022	0.0000	0.0017
avisitvis02	0.0713	0.0707	0.0006	0.0091
avisitvis03	0.0671	0.0672	−0.0001	−0.0014
avisitvis04	0.0760	0.0764	−0.0003	−0.0046
avisitvis05	0.0869	0.0863	0.0006	0.0066
avisitvis06	0.0866	0.0865	0.0001	0.0017
avisitvis07	0.1084	0.1071	0.0014	0.0129
avisitvis08	0.1159	0.1145	0.0014	0.0123
avisitvis09	0.1305	0.1283	0.0022	0.0168
avisitvis10	0.1409	0.1400	0.0009	0.0062
raceblack	0.0546	0.0550	−0.0003	−0.0058
racewhite	0.0521	0.0520	0.0001	0.0020
avisitvis01:armcdtrt	0.0634	0.0628	0.0006	0.0093
avisitvis02:armcdtrt	0.0806	0.0798	0.0007	0.0094
avisitvis03:armcdtrt	0.0922	0.0916	0.0005	0.0060
avisitvis04:armcdtrt	0.0996	0.1004	−0.0008	−0.0079
avisitvis05:armcdtrt	0.1159	0.1147	0.0013	0.0110
avisitvis06:armcdtrt	0.1191	0.1189	0.0002	0.0017
avisitvis07:armcdtrt	0.1390	0.1382	0.0008	0.0061
avisitvis08:armcdtrt	0.1500	0.1474	0.0027	0.0180
avisitvis09:armcdtrt	0.1677	0.1644	0.0033	0.0200
avisitvis10:armcdtrt	0.1831	0.1815	0.0016	0.0091

BCVA data comparison between Bayesian and frequentist MMRMs

09/02/2024