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) |

introduction.R

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)
#> Loading required package: tern
#> Loading required package: rtables
#> Loading required package: formatters
#> Loading required package: magrittr
#>
#> Attaching package: 'rtables'
#> The following object is masked from '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
```