I fitted MMRMs to the orthodont dataset using both mmrm::mmrm() and brms.mmrm::brm_model(). Apologies, the reprex has become a bit lengthy. Of note, the endpoint here is a response, not a chg-variable. I noted a few things that are summarized below:

Prerequisites

packages <- c("tidyr",  "dplyr", "stringr", "tictoc", "nlme")
invisible(lapply(packages, library, character.only = T))
library(brms.mmrm) # version brms.mmrm_0.0.2.9001
library(mmrm) #  mmrm_0.3.3
set.seed(123L)

Data

# data
data(Orthodont)
# age variable: brms.mmrm needs a character or factor 'time' variable
Orthodont$age_fct <- as.character(Orthodont$age)
Orthodont$age_fct[which(Orthodont$age_fct == "8")] <- "08" # changed as otherwise age = "10" becomes the reference level
Orthodont$age_fct <- factor(Orthodont$age_fct)
levels(Orthodont$age_fct)
#> [1] "08" "10" "12" "14"
# sex variable: reference level changed since `level_control` in brm_data() did not seem be effective
levels(Orthodont$Sex)
#> [1] "Male"   "Female"
Orthodont$Sex <- relevel(Orthodont$Sex, ref = "Female")
levels(Orthodont$Sex)
#> [1] "Female" "Male"
lapply(Orthodont, class) |> unlist()
#>  distance       age  Subject1  Subject2       Sex   age_fct 
#> "numeric" "numeric" "ordered"  "factor"  "factor"  "factor"
head(Orthodont)
#> Grouped Data: distance ~ age | Subject
#>   distance age Subject  Sex age_fct
#> 1     26.0   8     M01 Male      08
#> 2     25.0  10     M01 Male      10
#> 3     29.0  12     M01 Male      12
#> 4     31.0  14     M01 Male      14
#> 5     21.5   8     M02 Male      08
#> 6     22.5  10     M02 Male      10

mmrm

mmrm_fit <- mmrm(
  formula = distance  ~ age_fct + Sex * age_fct + us(age_fct | Subject),
  data = Orthodont,
  control = mmrm_control(method = "Satterthwaite")
)
summary(mmrm_fit)
#> mmrm fit
#> 
#> Formula:     distance ~ age_fct + Sex * age_fct + us(age_fct | Subject)
#> Data:        Orthodont (used 108 observations from 27 subjects with maximum 4 
#> timepoints)
#> Covariance:  unstructured (10 variance parameters)
#> Method:      Satterthwaite
#> Vcov Method: Asymptotic
#> Inference:   REML
#> 
#> Model selection criteria:
#>      AIC      BIC   logLik deviance 
#>      434      447     -207      414 
#> 
#> Coefficients: 
#>                   Estimate Std. Error      df t value Pr(>|t|)    
#> (Intercept)        21.1818     0.7017 25.0000  30.189  < 2e-16 ***
#> age_fct10           1.0455     0.6155 25.0000   1.699 0.101793    
#> age_fct12           1.9091     0.6068 25.0000   3.146 0.004241 ** 
#> age_fct14           2.9091     0.6729 25.0000   4.323 0.000215 ***
#> SexMale             1.6932     0.9115 25.0000   1.858 0.075038 .  
#> age_fct10:SexMale  -0.1080     0.7995 25.0000  -0.135 0.893671    
#> age_fct12:SexMale   0.9347     0.7883 25.0000   1.186 0.246904    
#> age_fct14:SexMale   1.6847     0.8741 25.0000   1.927 0.065385 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Covariance estimate:
#>        08     10     12     14
#> 08 5.4155 2.7168 3.9102 2.7102
#> 10 2.7168 4.1848 2.9272 3.3172
#> 12 3.9102 2.9272 6.4557 4.1307
#> 14 2.7102 3.3172 4.1307 4.9857

brms.mmrm

# Pre-process dataset
brms.mmrm_data <- brm_data(
  data = Orthodont, 
  outcome = "distance",
  role = "response",
  group = "Sex",
  time = "age_fct",
  patient = "Subject",
  baseline = NULL,
  level_control = "Female", # an effect for 'sexMale' is shown for both "Male" and "Female"
  level_baseline = "08" # an error occurs if set to NULL
)
# Model formula
brms.mmrm_formula <- brm_formula(
  data = brms.mmrm_data,
  intercept = TRUE,
  effect_baseline = FALSE,
  effect_group = TRUE,
  effect_time = TRUE,
  interaction_baseline = TRUE,
  interaction_group = TRUE,
  correlation = "unstructured"
)
print(brms.mmrm_formula)
#> distance ~ age_fct + Sex + Sex:age_fct + unstr(time = age_fct, gr = Subject) 
#> sigma ~ 0 + age_fct
# Model fit
options(mc.cores = 3)
tic()
brms.mmrm_fit <- brm_model(
  data = brms.mmrm_data,
  formula = brms.mmrm_formula,
  chains = 3,
  iter = 2000,
  refresh = 0,
  cores = 3
  )
toc()
#> 155.63 sec elapsed

summary(brms.mmrm_fit)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = log 
#> Formula: distance ~ age_fct + Sex + Sex:age_fct + unstr(time = age_fct, gr = Subject) 
#>          sigma ~ 0 + age_fct
#>    Data: data (Number of observations: 108) 
#>   Draws: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 3000
#> 
#> Correlation Structures:
#>                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> cortime(X08,X10)     0.46      0.14     0.15     0.71 1.00     1766     1808
#> cortime(X08,X12)     0.55      0.13     0.25     0.77 1.00     1865     2199
#> cortime(X10,X12)     0.43      0.15     0.08     0.69 1.00     2007     2106
#> cortime(X08,X14)     0.39      0.16     0.05     0.66 1.00     1675     1901
#> cortime(X10,X14)     0.62      0.12     0.34     0.81 1.00     1907     1897
#> cortime(X12,X14)     0.62      0.12     0.35     0.81 1.00     2109     2260
#> 
#> Population-Level Effects: 
#>                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept             21.17      0.71    19.76    22.62 1.00     1516     1838
#> age_fctX10             1.04      0.69    -0.34     2.44 1.00     1273     1540
#> age_fctX12             1.90      0.69     0.50     3.24 1.00     1420     1871
#> age_fctX14             2.89      0.75     1.36     4.40 1.01     1214     1615
#> SexMale                1.71      0.93    -0.13     3.54 1.00     1546     1955
#> age_fctX10:SexMale    -0.10      0.90    -1.86     1.72 1.00     1227     1675
#> age_fctX12:SexMale     0.94      0.90    -0.78     2.77 1.00     1419     1944
#> age_fctX14:SexMale     1.70      0.98    -0.23     3.60 1.00     1185     1602
#> sigma_age_fctX08       0.85      0.14     0.60     1.13 1.00     2121     2342
#> sigma_age_fctX10       0.71      0.13     0.46     0.99 1.00     2002     2045
#> sigma_age_fctX12       0.92      0.14     0.68     1.20 1.00     1949     2055
#> sigma_age_fctX14       0.78      0.13     0.55     1.07 1.00     1843     2039

Some observations and possible issues

First, the estimates are quite close, also the standard errors! The covariances look different, but maybe I missed something here.
We might want to add to the documentation how the ordering and the reference of the time variable (character or factor variable) is determined and that appropriate values need to be chosen to obtain the desired parameterization.
Changing the reference level of the groupvariable through level_control in brm_data() did not seem have an effect.
When role = "response" in brm_data(), I understood from the documentation that no level_baseline = ... is required, which makes sense to me. However, I received an error when trying level_baseline = NULL. (With effect_baseline = FALSE and effect_time = TRUE in brm_formula() the parameterization looked good.)

Perhaps we can briefly discuss this at the next meeting.

I tried out both packages briefly last week on the FEV-1 dataset built in to mmrm. Similar findings. I can share my code here at some point.

One comment about this:

The covariances look different, but maybe I missed something here.

This threw me at first also. The brms model for the residual variance involves a linear model for log(sd). The coefficients reported for sigma_age_fctXXX are summaries of the posterior for the log(residual sd) at these levels of the covariates. If you exponentiate and then square the medians, you'll probably find they are relatively close to the diagonal elements of the "covariance estimate" output from mmrm.

Thanks for the explanation, @andrew-bean! Then it's close too. The covariances generally appear to be a bit smaller with the brms.mmrm fit.

Variances and covariances from mmrm fit


mmrm_fit["cov"]
#> $cov
#>          08       10       12       14
#> 08 5.415458 2.716820 3.910231 2.710230
#> 10 2.716820 4.184773 2.927161 3.317160
#> 12 3.910231 2.927161 6.455743 4.130742
#> 14 2.710230 3.317160 4.130742 4.985741

Variances from brms.mmrm fit

summary(brms.mmrm_fit)[["fixed"]] |> 
  as_tibble(rownames = NA) |> 
  rownames_to_column("parameter") |>
  filter(str_detect(parameter, "sigma_")) |>
  select(parameter, Estimate) |>
  mutate(Variance = exp(Estimate)^2)
#> # A tibble: 4 × 3
#>   parameter        Estimate Variance
#>   <chr>               <dbl>    <dbl>
#> 1 sigma_age_fctX08    0.848     5.45
#> 2 sigma_age_fctX10    0.716     4.18
#> 3 sigma_age_fctX12    0.912     6.20
#> 4 sigma_age_fctX14    0.784     4.80

Covariances from brms.mmrm fit


summary(brms.mmrm_fit)[["cor_pars"]] |> 
  as_tibble(rownames = NA) |> 
  rownames_to_column("cortime") |>
  select(cortime, Estimate) |>
  mutate(Covariance = exp(Estimate)^2)
#> # A tibble: 6 × 3
#>   cortime          Estimate Covariance
#>   <chr>               <dbl>      <dbl>
#> 1 cortime(X08,X10)    0.449       2.46
#> 2 cortime(X08,X12)    0.548       2.99
#> 3 cortime(X10,X12)    0.423       2.33
#> 4 cortime(X08,X14)    0.383       2.15
#> 5 cortime(X10,X14)    0.624       3.48
#> 6 cortime(X12,X14)    0.621       3.46

Thanks for your work on this, @chstock (and for your comments, @andrew-bean)!

A first thought: interesting choice of dataset. I think this is one where we really need the full flexibility of a fully parameterized covariance matrix. Toeplitz doesn't fit because cortime(X12,X14) = 0.621 seems to be meaningfully different from the other elements on its off-diagonal band (cortime(X10,X12) = 0.432 and cortime(X8,X10) = 0.449), and AR(1) would struggle even more because the next furthest off-diagonal band has even higher positive correlations (cortime(X10,X14) = 0.624 and cortime(X08,X12) = 0.548).

When role = "response" in brm_data(), I understood from the documentation that no level_baseline = ... is required, which makes sense to me. However, I received an error when trying level_baseline = NULL. (With effect_baseline = FALSE and effect_time = TRUE in brm_formula() the parameterization looked good.)

Sorry about the confusion. The error behavior is correct and the documentation is incorrect. When the outcome variable is the raw response, we need to identify a baseline level so we can produce marginal estimates of change from baseline. Should be fixed in 25bfd4e9a97099c3a463359151837ac36bbea6f1.

Changing the reference level of the groupvariable through level_control in brm_data() did not seem have an effect.

Yes, level_control and level_baseline only apply to downstream post-processing functions like brm_marginal_draws(). These arguments do not control the fixed effect parameterization that brms chooses based on the data and formula.

We might want to add to the documentation how the ordering and the reference of the time variable (character or factor variable) is determined and that appropriate values need to be chosen to obtain the desired parameterization.

Looking at https://github.com/paul-buerkner/brms/blob/ca12a94a39318c10cdebe96425dc578f9ad926ab/R/data-predictor.R#L977-L996, the fixed effect parameterization seems to come from a call to model.matrix() which is called from within brms::make_standata(). The X field of the output has the fixed effects model matrix.

For our MMRM models with continuous outcomes, the model.matrix() call has two arguments: the dataset, and a special formula decorated with attributes which do not seem to affect the processing. brms does not seem to expose other arguments, but users can set the contrasts option. From ?options:

contrasts:

the default contrasts used in model fitting such as with aov or lm. A character vector of length two, the first giving the function to be used with unordered factors and the second the function to be used with ordered factors. By default the elements are named c("unordered", "ordered"), but the names are unused.

1652e71e073f118f6d66177fe10a6b47b1a52ccb adds explanations in the Rd files and usage vignette.

> 1. First, the estimates are quite close, also the standard errors! The covariances look different, but maybe I missed something here.

I put your results together side-by side and added the exp(sigma) correction in a reprex (expand "details" belown to see the all the code).

```r suppressPackageStartupMessages({ library(tidyverse) library(stringr) library(nlme) library(mmrm) library(brms) library(brms.mmrm) library(posterior) }) set.seed(123L) data(Orthodont) Orthodont$age_fct <- as.character(Orthodont$age) Orthodont$age_fct[which(Orthodont$age_fct == "8")] <- "08" Orthodont$age_fct <- factor(Orthodont$age_fct) Orthodont$Sex <- relevel(Orthodont$Sex, ref = "Female") Orthodont$Subject <- as.character(Orthodont$Subject) str(Orthodont) mmrm_fit <- mmrm( formula = distance ~ age_fct + Sex * age_fct + us(age_fct | Subject), data = Orthodont, control = mmrm_control(method = "Satterthwaite") ) summary(mmrm_fit) mmrm_fixed <- summary(mmrm_fit)$coefficients %>% as.data.frame() %>% mutate(variable = tolower(rownames(.))) %>% mutate(variable = gsub("(", "", variable, fixed = TRUE)) %>% mutate(variable = gsub(")", "", variable, fixed = TRUE)) %>% rename(mmrm_estimate = Estimate, mmrm_se = `Std. Error`) %>% as_tibble() %>% select(variable, mmrm_estimate, mmrm_se) mmrm_variance <- tibble( variable = c("sigma_08", "sigma_10", "sigma_12", "sigma_14"), mmrm_estimate = sqrt(diag(mmrm_fit$cov)) ) diagnoal_factor <- diag(1 / sqrt(diag(mmrm_fit$cov))) correlation_matrix <- diagnoal_factor %*% mmrm_fit$cov %*% diagnoal_factor levels <- sort(unique(as.character(Orthodont$age_fct))) colnames(correlation_matrix) <- levels mmrm_correlation <- correlation_matrix %>% as.data.frame() %>% as_tibble() %>% mutate(x1 = levels) %>% pivot_longer( cols = -any_of("x1"), names_to = "x2", values_to = "mmrm_estimate" ) %>% filter(as.numeric(x1) < as.numeric(x2)) %>% mutate(variable = sprintf("correlation_%s_%s", x1, x2)) %>% select(variable, mmrm_estimate) mmrm_summary <- bind_rows(mmrm_fixed, mmrm_variance, mmrm_correlation) brms.mmrm_data <- brm_data( data = Orthodont, outcome = "distance", role = "response", group = "Sex", time = "age_fct", patient = "Subject", baseline = NULL, level_control = "Female", level_baseline = "08" ) brms.mmrm_formula <- brm_formula( data = brms.mmrm_data, intercept = TRUE, effect_baseline = FALSE, effect_group = TRUE, effect_time = TRUE, interaction_baseline = TRUE, interaction_group = TRUE, correlation = "unstructured" ) options(mc.cores = 3) brms.mmrm_fit <- brm_model( data = brms.mmrm_data, formula = brms.mmrm_formula, chains = 3, iter = 2000, refresh = 0, cores = 3 ) brms_draws <- brms.mmrm_fit %>% as_draws_df() for (level in levels) { name <- paste0("b_sigma_age_fctX", level) brms_draws[[name]] <- exp(brms_draws[[name]]) } brms_summary <- brms_draws %>% summarize_draws() %>% select(variable, mean, sd) %>% filter(!(variable %in% c("lprior", "lp__"))) %>% rename(brms_estimate = mean, brms_se = sd) %>% mutate( variable = variable %>% tolower() %>% gsub(pattern = "b_", replacement = "") %>% gsub(pattern = "sigma_age_fct", replacement = "sigma_") %>% gsub(pattern = "cortime", replacement = "correlation") %>% gsub(pattern = "_x", replacement = "_") %>% gsub(pattern = "fctx", replacement = "fct") %>% gsub(pattern = "__", replacement = "_") ) results <- left_join(x = brms_summary, y = mmrm_summary, by = "variable") %>% mutate( diff_estimate = brms_estimate - mmrm_estimate, diff_se = brms_se - mmrm_se ) %>% select(variable, ends_with("estimate"), ends_with("se")) print(results) ```

I got these results:

# A tibble: 18 × 7
   variable          brms_estimate mmrm_estimate diff_estimate brms_se mmrm_se diff_se
   <chr>                     <num>         <num>         <num>   <num>   <num>   <num>
 1 intercept               21.2           21.2         0.00157   0.733   0.702  0.0316
 2 age_fct10                1.03           1.05       -0.0142    0.725   0.615  0.110 
 3 age_fct12                1.89           1.91       -0.0201    0.722   0.607  0.115 
 4 age_fct14                2.88           2.91       -0.0251    0.786   0.673  0.113 
 5 sexmale                  1.68           1.69       -0.0119    0.951   0.911  0.0395
 6 age_fct10:sexmale       -0.0874        -0.108       0.0206    0.928   0.799  0.129 
 7 age_fct12:sexmale        0.944          0.935       0.00960   0.956   0.788  0.167 
 8 age_fct14:sexmale        1.70           1.68        0.0163    1.02    0.874  0.145 
 9 sigma_08                 2.38           2.33        0.0542    0.357  NA     NA     
10 sigma_10                 2.06           2.05        0.0104    0.286  NA     NA     
11 sigma_12                 2.52           2.54       -0.0182    0.345  NA     NA     
12 sigma_14                 2.21           2.23       -0.0241    0.302  NA     NA     
13 correlation_08_10        0.445          0.571      -0.126     0.152  NA     NA     
14 correlation_08_12        0.541          0.661      -0.120     0.139  NA     NA     
15 correlation_10_12        0.420          0.563      -0.143     0.150  NA     NA     
16 correlation_08_14        0.377          0.522      -0.145     0.157  NA     NA     
17 correlation_10_14        0.620          0.726      -0.107     0.117  NA     NA     
18 correlation_12_14        0.620          0.728      -0.108     0.118  NA     NA

So yes, I agree the results look close overall. Also, relative to mmrm, brms.mmrm seems to systematically underestimate the correlations and overestimate the standard errors of fixed effects.

Any outstanding issues with this? Do you think it is a good time to start a vignette with these real data comparisons?

When role = "response" in brm_data(), I understood from the documentation that no level_baseline = ... is required, which makes sense to me. However, I received an error when trying level_baseline = NULL. (With effect_baseline = FALSE and effect_time = TRUE in brm_formula() the parameterization looked good.)

Sorry about the confusion. The error behavior is correct and the documentation is incorrect. When the outcome variable is the raw response, we need to identify a baseline level so we can produce marginal estimates of change from baseline. Should be fixed in 25bfd4e.

Thanks very much, @wlandau! The documentation is clear now. Still, I wonder what, from a user perspective, would trigger a model of the 'raw response', apart from the convenience if this is how the data are available (which isn't a good rationale). If I am interested in marginal estimates of change from baseline, wouldn't I then model 'change from baseline' directly? And if I decide to model 'raw response' longitudinally, am I then really interested in change from baseline? In my experience, in a clinical trial setting mostly 'change from baseIine' is modeled. I would expect someone who decides to model 'raw response' to have an interest in the treatment effect and the marginal means of the response. If it is just an option to also obtain marginal means of the change from baseline, this is fine of course. Apologies, in case I am missing something (technical) here.

Related: There is a general equivalance between both options w.r.t. treatment effect estimation, as mentioned e.g. in this EMA guideline [section 5.6.]. From the guideline:

Note that when the baseline is included as a covariate in a standard linear model, the estimated treatment effects are identical for both ‘change from baseline’ (on an additive scale) and the ‘raw outcome’ analysis. Consequently if the appropriate adjustment is done, then the choice of endpoint becomes solely an issue of interpretability.

Maybe it would be helpful and not too prescriptive if we add a note to the documentation that adjustment for baseline values is generally recommended, and perhaps that both options, 'raw response' and 'change from baseline', when adjusted for the baseline yield the same treatment effect estimates (the latter to be verified).

Any outstanding issues with this? Do you think it is a good time to start a vignette with these real data comparisons?

Sounds good to me. Let's perhaps align on Wednesday on datasets and (broadly) model specifications.

Still, I wonder what, from a user perspective, would trigger a model of the 'raw response', apart from the convenience if this is how the data are available (which isn't a good rationale).

When modeling the raw response, the formula can be specified so that placebo and treatment are constrained to have the same mean at baseline. This extra assumption is reasonable in many trials, and it may improve precision in some cases. My colleagues use this kind of model and refer to it as a "cLDA" (constrained longitudinal data analysis) model.

If I am interested in marginal estimates of change from baseline, wouldn't I then model 'change from baseline' directly?

Maybe, but treatment differences are contrasts as well, which still leaves us in a situation where we are modeling x and reporting contrasts of x.

And if I decide to model 'raw response' longitudinally, am I then really interested in change from baseline? In my experience, in a clinical trial setting mostly 'change from baseIine' is modeled.

Yes, and this is what the Novartis vignettes recommend as well.

I would expect someone who decides to model 'raw response' to have an interest in the treatment effect and the marginal means of the response. If it is just an option to also obtain marginal means of the change from baseline, this is fine of course.

Good news is that the package supports all of the above. If raw response is modeled, then the package gets you marginals of the response, change from baseline, and treatment differences of change from baseline. If change from baseline is modeled, then you get marginals of change from baseline and treatment differences of change from baseline.

Maybe it would be helpful and not too prescriptive if we add a note to the documentation that adjustment for baseline values is generally recommended, and perhaps that both options, 'raw response' and 'change from baseline', when adjusted for the baseline yield the same treatment effect estimates (the latter to be verified).

I agree it is worth discussing how opinionated we want brms.mmrm to be, and which opinions it should have, even if it tries to make all the major options possible and convenient. I agree that baseline should be a covariate when change from baseline is the outcome variable, and citing the findings and recommendations of regulatory guidances (that EMA one in particular) will be helpful.

@wlandau, great answers, thanks!!

This issue can be closed.

openpharma / brms.mmrm

Comparison of brms.mmrm and mmrm: orthodont dataset #66

> 1. First, the estimates are quite close, also the standard errors! The covariances look different, but maybe I missed something here.