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.
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)
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."))
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 |
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.")
Characteristic | CTL, N = 4941 | TRT, N = 5061 |
---|---|---|
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) |
1 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.")
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()
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).
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
This section compares the Bayesian posterior parameter estimates from
brms.mmrm
to the frequentist parameter estimates of the
mmrm
package.
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 = "_")
+ )
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())
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"
+ )
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"
+ )
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 |