I was recently made aware of the release of the mmrm package in R. It has been developed by a group of programmers and statisticians at a number of pharmaceutical companies, led by Daniel Sabanes Bove at Roche, as part of the ASA Biopharmaceutical Section Software Engineering Working Group. I’ve written previously about fitting mixed models for repeated measures (MMRM) using R, Stata and SAS. In R, this can be done using the gls function in the nlme package, but there are a number of limitations with this approach. For example, it is difficult (or impossible) to fit models where you allow the covariance parameters to be distinct between treatment groups. In this post, I’ll take a very quick look at the new mmrm package in R.

To demonstrate the typical use of the mmrm package, I reload the simulated dataset generated in this previous post. This is a dataset with a baseline measurement y0 and three follow-up measurements, as shown in the following summary of the wide version of the data:

summary(wideData) id trt y0 y1 y2 y3 Min. : 1.0 Min. :0.000 Min. :-3.00794 Min. :-3.2529 Min. :-3.0906 Min. :-2.8766 1st Qu.:125.8 1st Qu.:0.000 1st Qu.:-0.85620 1st Qu.:-0.6638 1st Qu.:-0.4164 1st Qu.:-0.2294 Median :250.5 Median :0.000 Median :-0.08929 Median : 0.2174 Median : 0.5147 Median : 0.7423 Mean :250.5 Mean :0.488 Mean :-0.03882 Mean : 0.1731 Mean : 0.4559 Mean : 0.7058 3rd Qu.:375.2 3rd Qu.:1.000 3rd Qu.: 0.86261 3rd Qu.: 1.0390 3rd Qu.: 1.4003 3rd Qu.: 1.6568 Max. :500.0 Max. :1.000 Max. : 4.40421 Max. : 3.6834 Max. : 4.1396 Max. : 4.6730 NA's :42 NA's :75 NA's :113

## MMRM using the nlme package

To fit this model using the older nlme package, in my previous code I used the following syntax:

library(nlme) mmrm <- nlme::gls(y~y0*factor(time)+trt*factor(time), na.action=na.omit, data=longData, correlation=nlme::corSymm(form=~time | id), weights=nlme::varIdent(form=~1|time))

The values passed to the correlation and weights arguments are needed to specify that we wanted an unstructured variance covariance matrix across the time points (visits).

## MMRM using the mmrm package

To fit the same model using the new mmrm package, I discovered I first needed to convert the id, time and trt variables to factors. I then fitted the MMRM model using:

library(mmrm) mmrm2 <- mmrm(y~y0*time+trt*time+us(time|id), data=longData)

Note the simpler syntax for specifying an unstructured covariance matrix. Comparing the estimates (not shown here) I obtained from this fit to those from gls, they were essentially identical, with the tiny differences presumably due to differences in the way the two functions/packages are maximising the likelihood function.

## Separate covariance matrices by treatment arm

One of the things that is not possible to do with the nlme/gls route is to allow the parameters of the residual covariance matrix to differ by treatment arms. This is something that is very easy to do in the mmrm package by modifying the us() argument as follows:

mmrm3 <- mmrm(y~y0*time+trt*time+us(time| trt/id), data=longData)

I have only just brushed the surface of what the package is capable I think, but it looks like a fantastic development for anyone who is interested in fitting mixed models to repeated measures data. For more on the package, see this blog post by the working group who developed the package and of course the package's page at CRAN, from which there are links to various vignettes and the package manual.

Plus it gives the exact Satterthwaite degrees of freedom rather than the approximate ones given by the emmeans (it was good, when there’s nothing else!). Soon also the Kenward-Roger ones will be implemented, what haven’t happened for the nlme::gls. The authors also announces plans on adding the robust estimators of the (co)variance, based on the clubSandwich package. Finally we got a tool that reproduces SAS properly!

Thanks Cairo! Yes, meanwhile Kenward-Roger adjustment is implemented (see version 0.2.2 on CRAN). Feel free to try it out and give us feedback!

Cheers, Daniel