2.1 Pre-processing

This analysis exercise uses thefev_dat datasetcontained in themmrm-package:

> data(fev_data, package = "mmrm")

It is an artificial (simulated) dataset of a clinical trialinvestigating the effect of an active treatment onFEV1(forced expired volume in one second), compared to placebo.FEV1 is a measure of how quickly the lungs can be emptiedand low levels may indicate chronic obstructive pulmonary disease(COPD).

The dataset is atibble with 800 rows and the followingnotable variables:

• USUBJID (subject ID)
• AVISIT (visit number, factor)
• VISITN (visit number, numeric)
• ARMCD (treatment,TRT orPBO)
• RACE (3-category race)
• SEX (female or male)
• FEV1_BL (FEV1 at baseline, %)
• FEV1 (FEV1 at study visits)
• WEIGHT (weighting variable)

The primary endpoint for the analysis is change from baseline inFEV1, which we derive below and denoteFEV1_CHG.

> fev_data <- fev_data |>+   mutate("FEV1_CHG" = FEV1 - FEV1_BL)

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

> fev_data <- brm_data(+   data = fev_data,+   outcome = "FEV1_CHG",+   group = "ARMCD",+   time = "AVISIT",+   patient = "USUBJID",+   baseline = "FEV1_BL",+   reference_group = "PBO",+   covariates = c("RACE", "SEX")+ ) |>+   brm_data_chronologize(order = "VISITN")

The following table shows the first rows of the dataset.

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

2.2 Descriptivestatistics

Table of baseline characteristics:

> fev_data |>+   select(ARMCD, USUBJID, SEX, RACE, FEV1_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 of change from baseline in FEV1 over 52 weeks:

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

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

> fev_data |>+   group_by(ARMCD) |>+   ggplot(aes(x = AVISIT, y = FEV1_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 FEV1",+     fill = "Treatment"+   ) ++   scale_fill_manual(values = c("darkgoldenrod2", "coral2")) ++   theme_bw()

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

3.1 Bayesian model

The formula for the Bayesian model includes additive effects forbaseline, study visit, race, sex, and study-arm-by-visitinteraction.

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

We fit the model usingbrms.mmrm::brm_model(). To ensurea good basis of comparison with the frequentist model, we put anextremely diffuse prior on the intercept. The parameters already havediffuse flexible priors by default.

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

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

> summary(b_mmrm_fit)#>  Family: gaussian #>   Links: mu = identity; sigma = log #> Formula: FEV1_CHG ~ FEV1_BL + ARMCD:AVISIT + AVISIT + RACE + SEX + unstr(time = AVISIT, gr = USUBJID) #>          sigma ~ 0 + AVISIT#>    Data: data[!is.na(data[[attr(data, "brm_outcome")]]), ] (Number of observations: 537) #>   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(VIS1,VIS2)     0.36      0.08     0.18     0.52 1.00    48758    26086#> cortime(VIS1,VIS3)     0.14      0.10    -0.05     0.33 1.00    49018    26172#> cortime(VIS2,VIS3)     0.04      0.10    -0.16     0.23 1.00    49178    25472#> cortime(VIS1,VIS4)     0.17      0.11    -0.06     0.38 1.00    49528    25555#> cortime(VIS2,VIS4)     0.11      0.09    -0.07     0.28 1.00    49509    24007#> cortime(VIS3,VIS4)     0.01      0.10    -0.18     0.21 1.00    45294    24353#> #> Regression Coefficients:#>                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS#> Intercept                     24.34      1.41    21.60    27.09 1.00    43678#> FEV1_BL                       -0.84      0.03    -0.89    -0.78 1.00    56286#> AVISIT2                        4.80      0.82     3.20     6.40 1.00    31792#> AVISIT3                       10.37      0.83     8.73    12.01 1.00    29808#> AVISIT4                       15.20      1.33    12.60    17.83 1.00    35293#> RACEBlackorAfricanAmerican     1.41      0.59     0.27     2.55 1.00    46945#> RACEWhite                      5.46      0.63     4.23     6.69 1.00    47801#> SEXFemale                      0.35      0.51    -0.64     1.36 1.00    49193#> AVISITVIS1:ARMCDTRT            3.98      1.07     1.89     6.06 1.00    31646#> AVISITVIS2:ARMCDTRT            3.93      0.83     2.31     5.55 1.00    46665#> AVISITVIS3:ARMCDTRT            2.98      0.68     1.65     4.31 1.00    52457#> AVISITVIS4:ARMCDTRT            4.41      1.68     1.06     7.71 1.00    45307#> sigma_AVISITVIS1               1.83      0.06     1.71     1.95 1.00    50344#> sigma_AVISITVIS2               1.59      0.06     1.47     1.71 1.00    48248#> sigma_AVISITVIS3               1.33      0.06     1.21     1.46 1.00    48058#> sigma_AVISITVIS4               2.28      0.06     2.16     2.41 1.00    51078#>                            Tail_ESS#> Intercept                     25394#> FEV1_BL                       24494#> AVISIT2                       24396#> AVISIT3                       24188#> AVISIT4                       24810#> RACEBlackorAfricanAmerican    25405#> RACEWhite                     23816#> SEXFemale                     24919#> AVISITVIS1:ARMCDTRT           26255#> AVISITVIS2:ARMCDTRT           23809#> AVISITVIS3:ARMCDTRT           24705#> AVISITVIS4:ARMCDTRT           25026#> sigma_AVISITVIS1              26156#> sigma_AVISITVIS2              24526#> sigma_AVISITVIS3              24328#> sigma_AVISITVIS4              23975#> #> 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 Frequentistmodel

The formula for the frequentist model is the same, except for thedifferent syntax for specifying the covariance structure of the MMRM. Wefit the model below.

> f_mmrm_fit <- mmrm::mmrm(+   formula = FEV1_CHG ~ FEV1_BL + ARMCD:AVISIT + AVISIT + RACE + SEX ++     us(AVISIT | USUBJID),+   data = mutate(+     fev_data,+     AVISIT = factor(as.character(AVISIT), ordered = FALSE)+   )+ )

The parameter summaries of the frequentist model are below.

> summary(f_mmrm_fit)#> mmrm fit#> #> Formula:     #> FEV1_CHG ~ FEV1_BL + ARMCD:AVISIT + AVISIT + RACE + SEX + us(AVISIT |  #>     USUBJID)#> Data:        #> mutate(fev_data, AVISIT = factor(as.character(AVISIT), ordered = FALSE)) (used #> 537 observations from 197 subjects with maximum 4 timepoints)#> Covariance:  unstructured (10 variance parameters)#> Method:      Satterthwaite#> Vcov Method: Asymptotic#> Inference:   REML#> #> Model selection criteria:#>      AIC      BIC   logLik deviance #>   3381.4   3414.2  -1680.7   3361.4 #> #> Coefficients: #>                                Estimate Std. Error        df t value Pr(>|t|)#> (Intercept)                    24.35372    1.40754 257.97000  17.302  < 2e-16#> FEV1_BL                        -0.84022    0.02777 190.27000 -30.251  < 2e-16#> AVISITVIS2                      4.79036    0.79848 144.82000   5.999 1.51e-08#> AVISITVIS3                     10.36601    0.81318 157.08000  12.748  < 2e-16#> AVISITVIS4                     15.19231    1.30857 139.25000  11.610  < 2e-16#> RACEBlack or African American   1.41921    0.57874 169.56000   2.452 0.015211#> RACEWhite                       5.45679    0.61626 157.54000   8.855 1.65e-15#> SEXFemale                       0.33812    0.49273 166.43000   0.686 0.493529#> AVISITVIS1:ARMCDTRT             3.98329    1.04540 142.32000   3.810 0.000206#> AVISITVIS2:ARMCDTRT             3.93076    0.81351 142.26000   4.832 3.46e-06#> AVISITVIS3:ARMCDTRT             2.98372    0.66567 129.61000   4.482 1.61e-05#> AVISITVIS4:ARMCDTRT             4.40400    1.66049 132.88000   2.652 0.008970#>                                  #> (Intercept)                   ***#> FEV1_BL                       ***#> AVISITVIS2                    ***#> AVISITVIS3                    ***#> AVISITVIS4                    ***#> RACEBlack or African American *  #> RACEWhite                     ***#> SEXFemale                        #> AVISITVIS1:ARMCDTRT           ***#> AVISITVIS2:ARMCDTRT           ***#> AVISITVIS3:ARMCDTRT           ***#> AVISITVIS4:ARMCDTRT           ** #> ---#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#> #> Covariance estimate:#>         VIS1    VIS2    VIS3    VIS4#> VIS1 37.8301 11.3255  3.4796 10.6844#> VIS2 11.3255 23.5476  0.7760  5.5103#> VIS3  3.4796  0.7760 13.8037  0.5683#> VIS4 10.6844  5.5103  0.5683 92.9625

4.1 Extract estimatesfrom Bayesian model

We extract and standardize the Bayesian estimates.

> b_mmrm_draws <- b_mmrm_fit |>+   as_draws_df()> visit_levels <- sort(unique(as.character(fev_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("Intercept", "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 = "_") |>+       gsub(pattern = "avisitvis", replacement = "avisit")+   )

4.2 Extract estimatesfrom 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)) |>+   mutate(variable = gsub("avisitvis", "avisit", variable)) |>+   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() |>+     gsub(pattern = "avisitvis", replacement = "avisit"),+   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 eachmodel and the differences between estimates (Bayesian minusfrequentist). The second table shows the standard errors of theseestimates and differences between standard errors. In each table, the“Relative” column shows the relative difference (the difference dividedby the frequentist quantity).

Because of the different statistical paradigms and estimationprocedures, especially regarding the covariance parameters, it would notbe realistic to expect the Bayesian and frequentist approaches to yieldvirtually identical results. Nevertheless, the absolute and relativedifferences in the table below show strong agreement betweenbrms.mmrm andmmrm.

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