BCVA data comparison between Bayesian and frequentist MMRMs

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.

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)

Data

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

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.

Fitting MMRMs

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    64232    23380
#> cortime(VIS01,VIS03)     0.31      0.03     0.25     0.37 1.00    63974    26624
#> cortime(VIS02,VIS03)     0.05      0.03    -0.02     0.11 1.00    70110    24587
#> cortime(VIS01,VIS04)     0.21      0.03     0.15     0.27 1.00    42555    26079
#> cortime(VIS02,VIS04)     0.13      0.03     0.07     0.20 1.00    50575    25891
#> cortime(VIS03,VIS04)    -0.01      0.03    -0.07     0.05 1.00    50372    25192
#> cortime(VIS01,VIS05)     0.17      0.03     0.11     0.23 1.00    44730    27906
#> cortime(VIS02,VIS05)     0.12      0.03     0.05     0.18 1.00    50419    27572
#> cortime(VIS03,VIS05)    -0.01      0.03    -0.07     0.06 1.00    53185    27533
#> cortime(VIS04,VIS05)     0.38      0.03     0.32     0.43 1.00    50371    27017
#> cortime(VIS01,VIS06)     0.26      0.03     0.20     0.32 1.00    41105    27790
#> cortime(VIS02,VIS06)     0.20      0.03     0.14     0.27 1.00    46110    26785
#> cortime(VIS03,VIS06)     0.04      0.03    -0.02     0.11 1.00    50559    27215
#> cortime(VIS04,VIS06)     0.40      0.03     0.35     0.46 1.00    51210    27168
#> cortime(VIS05,VIS06)     0.39      0.03     0.33     0.45 1.00    54492    24829
#> cortime(VIS01,VIS07)     0.07      0.04    -0.00     0.13 1.00    64797    24890
#> cortime(VIS02,VIS07)     0.09      0.03     0.02     0.15 1.00    66962    24910
#> cortime(VIS03,VIS07)    -0.00      0.04    -0.07     0.07 1.00    61101    24110
#> cortime(VIS04,VIS07)     0.15      0.03     0.08     0.21 1.00    62311    22787
#> cortime(VIS05,VIS07)     0.19      0.03     0.13     0.26 1.00    62582    24390
#> cortime(VIS06,VIS07)     0.21      0.04     0.14     0.28 1.00    68580    24400
#> cortime(VIS01,VIS08)     0.05      0.04    -0.02     0.12 1.00    66066    23819
#> cortime(VIS02,VIS08)     0.10      0.04     0.03     0.17 1.00    69331    24606
#> cortime(VIS03,VIS08)    -0.03      0.04    -0.10     0.04 1.00    65051    23999
#> cortime(VIS04,VIS08)     0.17      0.03     0.10     0.24 1.00    70470    24091
#> cortime(VIS05,VIS08)     0.17      0.04     0.10     0.24 1.00    72172    24032
#> cortime(VIS06,VIS08)     0.16      0.04     0.09     0.23 1.00    68767    24574
#> cortime(VIS07,VIS08)     0.05      0.04    -0.02     0.13 1.00    68724    22432
#> cortime(VIS01,VIS09)     0.03      0.04    -0.04     0.10 1.00    70365    23435
#> cortime(VIS02,VIS09)    -0.01      0.04    -0.08     0.07 1.00    67029    24108
#> cortime(VIS03,VIS09)    -0.04      0.04    -0.12     0.03 1.00    65481    22947
#> cortime(VIS04,VIS09)     0.12      0.04     0.04     0.19 1.00    70778    24290
#> cortime(VIS05,VIS09)     0.09      0.04     0.01     0.16 1.00    67795    23738
#> cortime(VIS06,VIS09)     0.17      0.04     0.09     0.24 1.00    68382    23236
#> cortime(VIS07,VIS09)     0.02      0.04    -0.06     0.09 1.00    60320    23116
#> cortime(VIS08,VIS09)     0.06      0.04    -0.02     0.14 1.00    66326    24022
#> cortime(VIS01,VIS10)     0.02      0.04    -0.06     0.10 1.00    56720    25693
#> cortime(VIS02,VIS10)     0.13      0.04     0.05     0.20 1.00    57584    25219
#> cortime(VIS03,VIS10)     0.02      0.04    -0.06     0.10 1.00    60710    25228
#> cortime(VIS04,VIS10)     0.31      0.04     0.24     0.38 1.00    56727    25218
#> cortime(VIS05,VIS10)     0.24      0.04     0.16     0.31 1.00    63407    25064
#> cortime(VIS06,VIS10)     0.30      0.04     0.22     0.37 1.00    59518    24562
#> cortime(VIS07,VIS10)     0.06      0.04    -0.03     0.15 1.00    68160    23806
#> cortime(VIS08,VIS10)     0.09      0.04     0.01     0.18 1.00    71536    24639
#> cortime(VIS09,VIS10)     0.08      0.05    -0.01     0.17 1.00    72879    23828
#> 
#> 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    52968    23424
#> BCVA_BL                 -0.00      0.00    -0.01     0.00 1.00    55733    22421
#> AVISITVIS02              0.28      0.07     0.14     0.42 1.00    29523    26038
#> AVISITVIS03              0.46      0.07     0.33     0.59 1.00    43908    25273
#> AVISITVIS04              0.86      0.08     0.71     1.01 1.00    27128    26680
#> AVISITVIS05              0.97      0.09     0.80     1.14 1.00    27041    26249
#> AVISITVIS06              1.34      0.09     1.17     1.51 1.00    28423    26075
#> AVISITVIS07              1.42      0.11     1.21     1.63 1.00    33236    26105
#> AVISITVIS08              1.71      0.12     1.49     1.94 1.00    33031    26727
#> AVISITVIS09              2.00      0.13     1.75     2.25 1.00    37527    27192
#> AVISITVIS10              2.10      0.14     1.82     2.38 1.00    31905    25606
#> RACEBlack                1.04      0.06     0.93     1.15 1.00    50375    27044
#> RACEWhite                2.01      0.05     1.90     2.11 1.00    48898    25367
#> AVISITVIS01:ARMCDTRT     0.54      0.06     0.42     0.66 1.00    33402    25937
#> AVISITVIS02:ARMCDTRT     0.73      0.08     0.57     0.88 1.00    45687    26156
#> AVISITVIS03:ARMCDTRT     1.01      0.09     0.83     1.19 1.00    44466    25516
#> AVISITVIS04:ARMCDTRT     1.10      0.10     0.91     1.30 1.00    35139    25412
#> AVISITVIS05:ARMCDTRT     1.38      0.11     1.16     1.61 1.00    34453    27482
#> AVISITVIS06:ARMCDTRT     1.63      0.12     1.39     1.86 1.00    34605    26525
#> AVISITVIS07:ARMCDTRT     2.02      0.14     1.74     2.29 1.00    42883    25319
#> AVISITVIS08:ARMCDTRT     2.35      0.15     2.05     2.64 1.00    42545    25257
#> AVISITVIS09:ARMCDTRT     2.66      0.17     2.33     2.98 1.00    42305    25539
#> AVISITVIS10:ARMCDTRT     3.07      0.18     2.72     3.43 1.00    39524    25368
#> sigma_AVISITVIS01       -0.01      0.02    -0.05     0.03 1.00    62714    25549
#> sigma_AVISITVIS02        0.23      0.02     0.18     0.27 1.00    70306    23424
#> sigma_AVISITVIS03        0.36      0.02     0.31     0.40 1.00    65928    25658
#> sigma_AVISITVIS04        0.44      0.02     0.40     0.49 1.00    54413    26689
#> sigma_AVISITVIS05        0.57      0.02     0.52     0.61 1.00    62557    26764
#> sigma_AVISITVIS06        0.58      0.02     0.53     0.63 1.00    53370    26108
#> sigma_AVISITVIS07        0.69      0.03     0.64     0.74 1.00    71710    24191
#> sigma_AVISITVIS08        0.74      0.03     0.69     0.79 1.00    74050    24224
#> sigma_AVISITVIS09        0.80      0.03     0.75     0.85 1.00    75622    23406
#> sigma_AVISITVIS10        0.84      0.03     0.79     0.90 1.00    63708    25100
#> 
#> 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

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.085  < 2e-16 ***
#> BCVA_BL              -9.934e-04  2.156e-03  9.905e+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.636e-02  9.796e+02  11.222  < 2e-16 ***
#> AVISITVIS05           9.638e-01  8.634e-02  9.630e+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.944  < 2e-16 ***
#> AVISITVIS09           1.996e+00  1.283e-01  7.784e+02  15.549  < 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.891  < 2e-16 ***
#> RACEWhite             2.005e+00  5.198e-02  9.768e+02  38.573  < 2e-16 ***
#> AVISITVIS01:ARMCDTRT  5.391e-01  6.282e-02  9.859e+02   8.582  < 2e-16 ***
#> AVISITVIS02:ARMCDTRT  7.248e-01  7.984e-02  9.804e+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.714  < 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.172  < 2e-16 ***
#> AVISITVIS10:ARMCDTRT  3.072e+00  1.815e-01  6.620e+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.9713  0.0630  0.4371  0.3315  0.3056 0.4688  0.1325  0.1020  0.0611 0.0587
#> VIS02 0.0630  1.5618  0.0871  0.2685  0.2635 0.4637  0.2181  0.2777 -0.0153 0.3761
#> VIS03 0.4371  0.0871  2.0221 -0.0216 -0.0190 0.1102 -0.0048 -0.0993 -0.1322 0.0719
#> VIS04 0.3315  0.2685 -0.0216  2.4114  1.0476 1.1410  0.4626  0.5660  0.4085 1.1479
#> VIS05 0.3056  0.2635 -0.0190  1.0476  3.0916 1.2594  0.6910  0.6308  0.3594 0.9999
#> VIS06 0.4688  0.4637  0.1102  1.1410  1.2594 3.1853  0.7541  0.6094  0.6823 1.2559
#> VIS07 0.1325  0.2181 -0.0048  0.4626  0.6910 0.7541  3.9273  0.2306  0.0723 0.3020
#> VIS08 0.1020  0.2777 -0.0993  0.5660  0.6308 0.6094  0.2306  4.3272  0.2683 0.4660
#> VIS09 0.0611 -0.0153 -0.1322  0.4085  0.3594 0.6823  0.0723  0.2683  4.8635 0.4140
#> VIS10 0.0587  0.3761  0.0719  1.1479  0.9999 1.2559  0.3020  0.4660  0.4140 5.3520

Comparison

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

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 = "_")
+   )

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

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.2871	4.2881	−0.0009	−0.0002
bcva_bl	−0.0010	−0.0010	0.0000	−0.0005
avisitvis02	0.2814	0.2810	0.0004	0.0015
avisitvis03	0.4579	0.4573	0.0006	0.0013
avisitvis04	0.8578	0.8570	0.0008	0.0010
avisitvis05	0.9656	0.9638	0.0018	0.0019
avisitvis06	1.3353	1.3339	0.0014	0.0010
avisitvis07	1.4179	1.4167	0.0012	0.0008
avisitvis08	1.7125	1.7107	0.0017	0.0010
avisitvis09	1.9965	1.9956	0.0009	0.0005
avisitvis10	2.1011	2.1005	0.0006	0.0003
raceblack	1.0386	1.0382	0.0003	0.0003
racewhite	2.0056	2.0051	0.0005	0.0003
avisitvis01:armcdtrt	0.5403	0.5391	0.0012	0.0023
avisitvis02:armcdtrt	0.7254	0.7248	0.0006	0.0009
avisitvis03:armcdtrt	1.0122	1.0115	0.0007	0.0007
avisitvis04:armcdtrt	1.1041	1.1042	0.0000	0.0000
avisitvis05:armcdtrt	1.3828	1.3834	−0.0006	−0.0004
avisitvis06:armcdtrt	1.6290	1.6301	−0.0011	−0.0007
avisitvis07:armcdtrt	2.0165	2.0160	0.0006	0.0003
avisitvis08:armcdtrt	2.3460	2.3469	−0.0010	−0.0004
avisitvis09:armcdtrt	2.6583	2.6585	−0.0002	−0.0001
avisitvis10:armcdtrt	3.0728	3.0723	0.0005	0.0002
sigma_avisitvis01	0.9895	0.9855	0.0039	0.0040
sigma_avisitvis02	1.2557	1.2497	0.0060	0.0048
sigma_avisitvis03	1.4289	1.4220	0.0068	0.0048
sigma_avisitvis04	1.5569	1.5529	0.0040	0.0026
sigma_avisitvis05	1.7632	1.7583	0.0049	0.0028
sigma_avisitvis06	1.7883	1.7847	0.0036	0.0020
sigma_avisitvis07	1.9930	1.9817	0.0112	0.0057
sigma_avisitvis08	2.0924	2.0802	0.0122	0.0058
sigma_avisitvis09	2.2203	2.2053	0.0150	0.0068
sigma_avisitvis10	2.3283	2.3134	0.0149	0.0064
correlation_vis01_vis02	0.0493	0.0512	−0.0019	−0.0367
correlation_vis01_vis03	0.3086	0.3119	−0.0033	−0.0107
correlation_vis02_vis03	0.0482	0.0490	−0.0008	−0.0169
correlation_vis01_vis04	0.2129	0.2166	−0.0037	−0.0172
correlation_vis02_vis04	0.1348	0.1383	−0.0035	−0.0256
correlation_vis03_vis04	−0.0107	−0.0098	−0.0009	0.0901
correlation_vis01_vis05	0.1724	0.1764	−0.0040	−0.0226
correlation_vis02_vis05	0.1165	0.1199	−0.0035	−0.0288
correlation_vis03_vis05	−0.0084	−0.0076	−0.0008	0.1023
correlation_vis04_vis05	0.3767	0.3837	−0.0069	−0.0180
correlation_vis01_vis06	0.2621	0.2665	−0.0044	−0.0166
correlation_vis02_vis06	0.2037	0.2079	−0.0042	−0.0200
correlation_vis03_vis06	0.0423	0.0434	−0.0011	−0.0262
correlation_vis04_vis06	0.4043	0.4117	−0.0074	−0.0180
correlation_vis05_vis06	0.3939	0.4013	−0.0074	−0.0186
correlation_vis01_vis07	0.0655	0.0679	−0.0024	−0.0348
correlation_vis02_vis07	0.0856	0.0881	−0.0025	−0.0283
correlation_vis03_vis07	−0.0021	−0.0017	−0.0004	0.2066
correlation_vis04_vis07	0.1461	0.1503	−0.0042	−0.0278
correlation_vis05_vis07	0.1934	0.1983	−0.0049	−0.0246
correlation_vis06_vis07	0.2080	0.2132	−0.0052	−0.0244
correlation_vis01_vis08	0.0481	0.0497	−0.0017	−0.0336
correlation_vis02_vis08	0.1042	0.1068	−0.0026	−0.0245
correlation_vis03_vis08	−0.0333	−0.0336	0.0003	−0.0088
correlation_vis04_vis08	0.1711	0.1752	−0.0041	−0.0233
correlation_vis05_vis08	0.1682	0.1725	−0.0042	−0.0246
correlation_vis06_vis08	0.1598	0.1641	−0.0044	−0.0267
correlation_vis07_vis08	0.0535	0.0559	−0.0024	−0.0429
correlation_vis01_vis09	0.0271	0.0281	−0.0010	−0.0367
correlation_vis02_vis09	−0.0062	−0.0056	−0.0007	0.1176
correlation_vis03_vis09	−0.0413	−0.0421	0.0009	−0.0212
correlation_vis04_vis09	0.1158	0.1193	−0.0035	−0.0296
correlation_vis05_vis09	0.0896	0.0927	−0.0031	−0.0331
correlation_vis06_vis09	0.1692	0.1733	−0.0042	−0.0240
correlation_vis07_vis09	0.0153	0.0165	−0.0013	−0.0778
correlation_vis08_vis09	0.0567	0.0585	−0.0018	−0.0303
correlation_vis01_vis10	0.0231	0.0257	−0.0027	−0.1040
correlation_vis02_vis10	0.1266	0.1301	−0.0035	−0.0267
correlation_vis03_vis10	0.0215	0.0219	−0.0003	−0.0149
correlation_vis04_vis10	0.3116	0.3195	−0.0079	−0.0249
correlation_vis05_vis10	0.2384	0.2458	−0.0074	−0.0303
correlation_vis06_vis10	0.2959	0.3042	−0.0083	−0.0272
correlation_vis07_vis10	0.0629	0.0659	−0.0030	−0.0451
correlation_vis08_vis10	0.0932	0.0968	−0.0036	−0.0374
correlation_vis09_vis10	0.0780	0.0811	−0.0032	−0.0390

> 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.1709	0.1709	−0.0001	−0.0005
bcva_bl	0.0022	0.0022	0.0000	−0.0007
avisitvis02	0.0709	0.0707	0.0003	0.0040
avisitvis03	0.0675	0.0672	0.0004	0.0054
avisitvis04	0.0763	0.0764	0.0000	−0.0006
avisitvis05	0.0867	0.0863	0.0004	0.0047
avisitvis06	0.0869	0.0865	0.0004	0.0050
avisitvis07	0.1078	0.1071	0.0007	0.0067
avisitvis08	0.1157	0.1145	0.0012	0.0105
avisitvis09	0.1302	0.1283	0.0018	0.0143
avisitvis10	0.1411	0.1400	0.0011	0.0075
raceblack	0.0551	0.0550	0.0002	0.0030
racewhite	0.0520	0.0520	0.0001	0.0011
avisitvis01:armcdtrt	0.0635	0.0628	0.0007	0.0117
avisitvis02:armcdtrt	0.0798	0.0798	−0.0001	−0.0007
avisitvis03:armcdtrt	0.0922	0.0916	0.0005	0.0058
avisitvis04:armcdtrt	0.1003	0.1004	−0.0001	−0.0009
avisitvis05:armcdtrt	0.1150	0.1147	0.0003	0.0028
avisitvis06:armcdtrt	0.1189	0.1189	0.0001	0.0007
avisitvis07:armcdtrt	0.1387	0.1382	0.0006	0.0042
avisitvis08:armcdtrt	0.1487	0.1474	0.0013	0.0086
avisitvis09:armcdtrt	0.1651	0.1644	0.0008	0.0046
avisitvis10:armcdtrt	0.1829	0.1815	0.0014	0.0076

27/03/2024