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 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).
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
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 | 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"
+ )
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 |