How to handle the residual covariance structure in a mixed model with repeated measurements?

Question

In clinical trials and other areas of applied statistics, we often need to model longitudinal data.

In most introductions to modelling longitudinal data with mixed effects models that I have seen, there seems to be very little focus on modelling the residual covariance structure. In R the most popular package for fitting mixed models seems to be lme4 and yet, as far as I know, it uses a compound symmetric (CS) structure with no option to specify anyything else. I know there are other packages in R (glmmTMB, mmrm and nlme for example) and other software entirely that allows other structures.

There is an expectation that not only will within-person repeated measures be correlated, but the correlations will not be equal accross all time points, but rather dimishish, the further apart the measurements are. Such a structure can be modelled by an autoregressive AR(1) structure. Another popular choice is of an unstructured (US) covariance matrix, where all the individual correlations are estimated seperately. I can see that with an US covariance matrix the software will have to estimate many more parameters that with AR(1) or CS structure, so it is far more flexble at the cost of being less efficient.

In many discussions on here and elsewhere, I have seen the recommendation to use AR(1), such as this answer by Frank Harrell:

What standard error should I use with correlated clusters in maximum likelihood estimation of multinomial logit

It is important to explicitly model the autocorrelation. You might model cluster effects as random intercepts, but within-cluster effects need to be handled with something like a continuous-time AR(1) correlation structure.

But then I read today this answer by Björn, who recommends an unstructured covariance matrix, and argued against AR(1):

Data from 2 eyes and repeated measures across time

you flexibly estimate how correlated different timepoints (and the eyes in the same person) are from the data without imposing something like a AR(1) structure (that tends to be very wrong for any real data, but if you assumed a more structured covariance structure that is appropriate you might have a gain in efficiency)

So how should we choose between, say AR(1) and US (or any other the other structures that some software supports such as Toeplitz, Heterogeneous AR(1) and Heterogeneous Toeplitz) ? Is it simply a matter of fitting several structures and using a significance test to choose between them ? If so then I assume there should be an adjustment for multiple testing ?

Answer 1

It is difficult, if not impossible, to determine the correct covariance structure in advance of fitting the model. Personally I would try unstructured and AR(1) first. The issue with an unstructured covariance matrix is that the software is required to estimate many parameters: $\frac{t(t+1)}{2}$ where $t$ is the number of time points, so this can obviously become quite large where there are many time points and can lead to problems with model convergence. On the other hand, CS and AR(1) both estimate just 2 parameters. So, practically speaking US might not really be an option. If that is the case then I would use AR(1) since we know that CS is likely not appropriate.

As for comparing models, I would tend to use AIC or BIC depending on the purpose of the model. As per Chakrabarti & Ghosh (2011) :

The Bayesian Information Criterion (BIC) is more useful in selecting a correct model while the AIC is more appropriate in finding the best model for predicting future observations.

So let's take a look at an example. Here we use the results of an RCT that consists of a completely randomised experimental design with data collected at 7 equally spaced points in time, with two factors, Program and Time. The response variable is strength and the repeated measures are within the Subj factor. The data are available in the SASmixed R package (it's called "Weights" there.

We fit 3 models, one with compound symmetry, one with AR(1) and one with an unstructured covariance matrix:

m0_cs <- mmrm(strength ~ Program + Time + Program:Time + cs(Time | Subj), data = dt)

m0_ar1 <- mmrm(strength ~ Program + Time + Program:Time + ar1(Time | Subj),  data = dt)

m0_us <- mmrm(strength ~ Program + Time + Program:Time + us(Time | Subj),  data = dt)

These converged normally and from them we obtain:

#            AIC      BIC   logLik deviance 
#
# CS      1424.8   1428.9   -710.4   1420.8 
# AR(1)   1270.8   1274.9   -633.4   1266.8 
# US      1290.9   1348.1   -617.4   1234.9 

So, we see that the AR(1) model would be selected by both AIC and BIC, so that makes life easy for us.

Let's also take a look at the covariance matrices that are estimated. First for compound symmetry:

         1       3       5       7       9      11      13
1  10.8002  9.6033  9.6033  9.6033  9.6033  9.6033  9.6033
3   9.6033 10.8002  9.6033  9.6033  9.6033  9.6033  9.6033
5   9.6033  9.6033 10.8002  9.6033  9.6033  9.6033  9.6033
7   9.6033  9.6033  9.6033 10.8002  9.6033  9.6033  9.6033
9   9.6033  9.6033  9.6033  9.6033 10.8002  9.6033  9.6033
11  9.6033  9.6033  9.6033  9.6033  9.6033 10.8002  9.6033
13  9.6033  9.6033  9.6033  9.6033  9.6033  9.6033 10.8002

Now for unstructured:

        1      3       5       7       9      11      13
1  8.7910 8.7685  8.9783  8.2104  8.6912  8.2333  8.4303
3  8.7685 9.4848  9.4762  8.5810  9.2148  8.7440  8.7013
5  8.9783 9.4762 10.7224  9.9401 10.6809 10.0846 10.2287
7  8.2104 8.5810  9.9401 10.0900 10.6132  9.9120 10.0570
9  8.6912 9.2148 10.6809 10.6132 12.1097 11.3586 11.3785
11 8.2333 8.7440 10.0846  9.9120 11.3586 11.7696 11.6642
13 8.4303 8.7013 10.2287 10.0570 11.3785 11.6642 12.7243

and lastly for AR(1)

         1       3       5       7       9      11      13
1  10.7600 10.2411  9.7473  9.2772  8.8298  8.4040  7.9988
3  10.2411 10.7600 10.2411  9.7473  9.2772  8.8298  8.4040
5   9.7473 10.2411 10.7600 10.2411  9.7473  9.2772  8.8298
7   9.2772  9.7473 10.2411 10.7600 10.2411  9.7473  9.2772
9   8.8298  9.2772  9.7473 10.2411 10.7600 10.2411  9.7473
11  8.4040  8.8298  9.2772  9.7473 10.2411 10.7600 10.2411
13  7.9988  8.4040  8.8298  9.2772  9.7473 10.2411 10.7600

References:

Chakrabarti, A., & Ghosh, J. K. (2011). AIC, BIC and recent advances in model selection. Philosophy of statistics, 583-605.