Introduction

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

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

Methodology

Statistical model

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

Estimation

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:

  1. establishes a reference grid where each cell represents a unique combination of the factor (covariate) levels,
  2. calculates the predicted marginal response for each cell, and then
  3. 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.

Inference

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.

Covariance structure

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\]

Model fitting

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.

Setup

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