Mixed effect Models with Repeated Measures (MMRM) are often used as the primary analysis of continuous longitudinal endpoints in clinical trials. In this setting, an MMRM is a specific linear mixed effects model that includes as fixed effects the variables: treatment arm, categorical visit, treatment by visit interaction, and other covariates for adjustment (e.g. age, gender). The covariance structure of the residuals can take on different forms. Often, an unstructured (i.e. saturated parametrization) covariance matrix is assumed which can be represented by random effects in the mixed model.

All of this has been implemented in proprietary software such as `SAS`

, whose `PROC MIXED`

routine is generally considered the gold standard for mixed models. However, this does not allow the use of interactive web applications to explore the clinical study data in a flexible way. Therefore, we wanted to implement MMRM in R in a way which could easily be embedded into a `shiny`

application. See the `teal.modules.clinical`

package for more details about using this code inside a `shiny`

application.

Descriptive statistics for a relevant analysis population (e.g. those patients with at least one post baseline visit) can be obtained by the functions described in Section @ref(baseline-tables).

Under the linear mixed model (`LMM`

) framework, an outcome vector \({\bf Y}\) is modeled as \[\bf Y = X \boldsymbol \beta + Z b + {\boldsymbol \varepsilon}\] where \(\bf X\) is the matrix for fixed effects, \(\bf Z\) is the matrix for random effects, \(\mathbf b \sim N(0, \mathbf D(\boldsymbol \theta))\), and \({\boldsymbol \varepsilon} \sim N(0, \mathbf R(\boldsymbol \theta))\). Letting \(\mathbf V = \mathbf Z \mathbf D \mathbf Z^T + \mathbf R\), the marginal and conditional models are then given by \(\mathbf Y \sim N(\bf X \boldsymbol \beta, V)\) and \(\mathbf Y | \mathbf b \sim N(\bf X \boldsymbol \beta + Z b, R)\), respectively.

An MMRM is a special case of a `LMM`

such that \(\bf Y\) is a collection of measurements made on a set of individuals over time, i.e. \(\mathbf Y = ({\bf Y}_1^T, {\bf Y}_2^T, ...)^T\) where \(\mathbf Y_i = \mathbf X_i \boldsymbol \beta + \mathbf Z_i \mathbf b_i + {\boldsymbol \varepsilon}_i\). In our context of clinical trials, individuals are patients identified by their unique subject `id`

, measurements are of treatment `response`

, and fixed effects include treatment `arm`

, categorical `visit`

, treatment by visit interaction, and potentially other `covariates`

(e.g. age, gender).

The parameters \(\boldsymbol \beta\), \(\bf b\), and \(\boldsymbol \theta\) can be estimated by maximizing the penalized and restricted maximum likelihoods. The average treatment response at each visit is often of particular interest. However, simply comparing predicted marginal responses \(E(Y_{ij})\) averaged across treatment x visit groups will not account for potential imbalance in covariates. (Imbalance may occur even in randomized trials due to patients dying or dropping out over time.) For a more fair comparison, **least-squares (LS) mean**:

- establishes a reference grid where each cell represents a unique combination of the factor (covariate) levels,
- calculates the predicted marginal response for each cell, and then
- takes a weighted average of the predicted marginal responses.

The following simple example illustrates the concept of LS means. Suppose we have a longitudinal clinical trial where three factors are considered: treatment (A and B), visit (1, 2, …), and gender (M and F). The marginal predicted response (sample size) for each reference cell is as follows.

A1 | B1 | … | |
---|---|---|---|

M | 100 (20) | 90 (35) | |

F | 50 (30) | 40 (15) |

The average predicted responses for arms A and B at Visit 1 are \((100 \times 20 + 50 \times 30) / 50 = 70\) and \((90 \times 35 + 40 \times 15) / 50 = 75\), respectively. The overall mean is higher in arm B than in arm A even though the mean is lower within each gender category. This seeming contradiction is caused by the imbalance in the data. LS mean calculates the weighted average across cells of the same treatment and visit. In this example, this is equivalent to taking the weighted average over each column. One may assign *equal* weights to each cell, i.e. \(0.5 \times 100 + 0.5 \times 50 = 75\) and \(0.5 \times 90 + 0.5 \times 40 = 65\). Alternatively, one may assign weights *proportional* to the observed frequencies of the factor combinations, i.e. \(0.55 \times 100 + 0.45 \times 50 = 77.5\) and \(0.55 \times 90 + 0.45 \times 40 = 67.5\). In both cases, the LS mean of response is lower in arm B than in arm A.

LS means are calculated in `tern.mmrm`

via the R package `emmeans`

. Users have the option to weigh marginal predicted responses with either `equal`

weights or `proportional`

weights. Note that for proportional weights, the weights are calculated at each visit by taking into account the observed frequencies of factor combinations at that time. Therefore, even though covariate imbalance may vary over time, LS mean provides an adjusted analysis of treatment response at all visits.

Performing inference on estimated parameters (e.g. calculating p-values) is less straightforward for MMRM. This is because the exact null distributions for parameter effects are unknown. `SAS`

addresses this issue by utilizing Satterthwaite’s method to approximate the adjusted degrees of freedom for \(F\) and \(t\) tests. `lme4`

and `lmerTest`

have also implemented Satterthwaite’s method. Unfortunately we found that these are not robust in their convergence behavior. Compared to `lme4`

, the R package `nlme`

can consider more flexible covariance structures. However, we have chosen not to use this package because it does not provide exact Satterthwaite adjusted degrees of freedom and the available approximate degrees of freedom can differ substantially. Therefore we built the new package `mmrm`

. With `mmrm`

, `tern.mmrm`

is able to replicate outputs from `SAS`

.

Users of `tern.mmrm`

have currently the following options for the covariance structure \(\mathbf V_i\):

- Unstructured: \[V_{ij} = \theta_{ij}\]
- Homogeneous AR(1): \[(\mathbf V_i)_{jk} = \sigma^2 \rho^{|j-k|}\]
- Heterogeneous AR(1): \[(\mathbf V_i)_{jk} = \sigma_j \sigma_k \rho^{|j-k|}\]
- Heterogeneous Toeplitz: \[(\mathbf V_i)_{jk}=\sigma_j \sigma_k \theta_{|j-k|}\]
- Heterogeneous Ante-Dependence: \[(\mathbf V_i)_{jk} = \sigma_j \sigma_k \prod_{i=j}^{k}\rho_i\]

In this section, we illustrate how to fit a MMRM with `tern.mmrm`

and how to fit a MMRM manually in R. We then compare with `SAS`

to show that the results match.

Our example dataset consists of several variables: subject ID (`USUBJID`

), visit number (`AVISIT`

), treatment (`ARMCD`

= `TRT`

or `PBO`

), 3-category race, sex, `FEV1`

at baseline (%), and `FEV1`

at study visits (%). `FEV1`

(forced expired volume in one second) is a measure of how quickly the lungs can be emptied. Low levels of `FEV1`

may indicate chronic obstructive pulmonary disease (`COPD`

). The scientific question at hand is whether treatment leads to an increase in `FEV1`

over time after adjusting for baseline covariates.

```
library(tern.mmrm)
#> Lade nötiges Paket: tern
#> Lade nötiges Paket: rtables
#> Lade nötiges Paket: formatters
#> Lade nötiges Paket: magrittr
#>
#> Attache Paket: 'rtables'
#> Das folgende Objekt ist maskiert 'package:utils':
#>
#> str
#> Registered S3 method overwritten by 'tern':
#> method from
#> tidy.glm broom
data(mmrm_test_data)
head(mmrm_test_data)
#> # A tibble: 6 × 7
#> USUBJID AVISIT ARMCD RACE SEX FEV1_BL FEV1
#> <fct> <fct> <fct> <fct> <fct> <dbl> <dbl>
#> 1 PT1 VIS1 TRT Black or African American Female 25.3 NA
#> 2 PT1 VIS2 TRT Black or African American Female 25.3 40.0
#> 3 PT1 VIS3 TRT Black or African American Female 25.3 NA
#> 4 PT1 VIS4 TRT Black or African American Female 25.3 20.5
```