# Can lme4 fit a repeated measures model with unstructured covariance? How does it differ from fitting using gls in nlme?

I recently read this paper by Donohue et al. in Pharmaceutical Statistics (Open access: https://doi.org/10.1002/pst.2285). Mixed model repeated measure (MMRM) analyses are discussed. In particular, one of the base models fit includes time points as categorical fixed effects (e.g., visit 1, 2, and 3 included as fixed effects), and their model adjusts for repeated measures by including no random intercept while allowing for different visit periods to be correlated for each subject, I think, using syntax such as (0 + Visit_Factor | Subject). This is not something I was used to seeing, and it indeed requires some over-ride argument for lmerControl, otherwise you get an error. I include a simplified and reproducible example that adapts their Supplemental Code to using the sleepstudy dataset included with lme4:

library(lme4)
data(sleepstudy)
sleep_data <- sleepstudy
sleep_data$Day_Factor <- as.factor(sleep_data$Days)
mod1 <- lmer(Reaction ~ Day_Factor + (0 + Day_Factor | Subject), data = sleep_data,
control = lmerControl(check.nobs.vs.nRE = "ignore"))
summary(mod1)
# Similar idea here, but Days no longer treated as categorical in the main fixed effect, but still has the same repeated measures adjustment
# mod1b <- lmer(Reaction ~ Days + (0 + Day_Factor | Subject), data = sleep_data,
#             control = lmerControl(check.nobs.vs.nRE = "ignore"))


The output is as follows:

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Day_Factor + (0 + Day_Factor | Subject)
Data: sleep_data
Control: lmerControl(check.nobs.vs.nRE = "ignore")

REML criterion at convergence: 1559.9

Scaled residuals:
Min       1Q   Median       3Q      Max
-1.09523 -0.21529  0.00571  0.18605  1.45012

Random effects:
Groups   Name        Variance Std.Dev. Corr
Subject  Day_Factor0 1009.68  31.776
Day_Factor1 1093.96  33.075   0.75
Day_Factor2  845.33  29.075   0.48 0.79
Day_Factor3 1480.94  38.483   0.47 0.75 0.89
Day_Factor4 1778.40  42.171   0.45 0.66 0.70 0.92
Day_Factor5 2648.64  51.465   0.37 0.53 0.50 0.73 0.86
Day_Factor6 3952.24  62.867   0.22 0.31 0.46 0.68 0.75 0.75
Day_Factor7 2480.48  49.804   0.50 0.48 0.60 0.60 0.70 0.69 0.71
Day_Factor8 3593.10  59.942   0.33 0.40 0.41 0.60 0.75 0.91 0.73 0.77
Day_Factor9 4458.02  66.768   0.52 0.55 0.43 0.57 0.72 0.84 0.46 0.66 0.89
Residual               20.32   4.508
Number of obs: 180, groups:  Subject, 18

Fixed effects:
Estimate Std. Error t value
(Intercept)  256.652      7.565  33.928
Day_Factor1    7.844      5.614   1.397
Day_Factor2    8.710      7.494   1.162
Day_Factor3   26.340      8.761   3.007
Day_Factor4   31.998      9.461   3.382
Day_Factor5   51.867     11.722   4.425
Day_Factor6   55.526     15.121   3.672
Day_Factor7   62.099     10.418   5.960
Day_Factor8   79.978     13.710   5.834
Day_Factor9   94.199     13.541   6.957

. . .

optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0073684 (tol = 0.002, component 1)


Questions: I think I understand that they are not simply using a random intercept in the model (i.e., (1 | Subject)) because that would imply equal correlation between each observation for a Subject (Day 1 and Day 9 for a given subject have the same correlation as Day 8 and Day 9), and an unstructured covariance assumes less inherently. Is the desired model sensible? Minor convergence issue aside, is this a valid approach to fitting the desired model? Without a random intercept and just adjusting for the repeated measures, it seems somewhat similar to GEE although restricted maximum likelihood is used here.

Using gls in the nlme package, I tried to fit what I think is an equivalent model.

library(nlme)
mod2 <- gls(Reaction ~ Day_Factor, data = sleep_data,
correlation = corSymm(form = ~ 1 | Subject))
summary(mod2)


With output:

Generalized least squares fit by REML
Model: Reaction ~ Day_Factor
Data: sleep_data
AIC      BIC    logLik
1698.281 1873.886 -793.1405

Correlation Structure: General
Formula: ~1 | Subject
Parameter estimate(s):
Correlation:
1      2      3      4      5      6      7      8      9
2   0.879
3   0.791  0.916
4   0.704  0.873  0.920
5   0.598  0.765  0.773  0.934
6   0.169  0.361  0.316  0.573  0.744
7   0.088  0.199  0.241  0.476  0.600  0.726
8   0.661  0.679  0.737  0.719  0.731  0.510  0.533
9  -0.172 -0.066 -0.081  0.157  0.350  0.775  0.683  0.342
10  0.300  0.376  0.290  0.459  0.607  0.784  0.465  0.535  0.776

Coefficients:
Value Std.Error   t-value p-value
(Intercept) 256.65181 12.341806 20.795320  0.0000
Day_Factor1   7.84395  6.065572  1.293192  0.1977
Day_Factor2   8.71009  7.988673  1.090306  0.2771
Day_Factor3  26.34021  9.503064  2.771759  0.0062
Day_Factor4  31.99762 11.069563  2.890594  0.0043
Day_Factor5  51.86665 15.906888  3.260641  0.0013
Day_Factor6  55.52645 16.666616  3.331597  0.0011
Day_Factor7  62.09878 10.155765  6.114633  0.0000
Day_Factor8  79.97770 18.896844  4.232331  0.0000
Day_Factor9  94.19942 14.604891  6.449854  0.0000

Correlation:
(Intr) Dy_Fc1 Dy_Fc2 Dy_Fc3 Dy_Fc4 Dy_Fc5 Dy_Fc6 Dy_Fc7 Dy_Fc8
Day_Factor1 -0.246
Day_Factor2 -0.324  0.773
Day_Factor3 -0.385  0.766  0.856
Day_Factor4 -0.448  0.653  0.662  0.915
Day_Factor5 -0.644  0.494  0.427  0.705  0.845
Day_Factor6 -0.675  0.349  0.414  0.658  0.755  0.844
Day_Factor7 -0.411  0.341  0.535  0.559  0.640  0.640  0.705
Day_Factor8 -0.766  0.302  0.303  0.530  0.673  0.901  0.855  0.676
Day_Factor9 -0.592  0.339  0.261  0.500  0.668  0.862  0.674  0.589  0.910

Standardized residuals:
Min           Q1          Med           Q3          Max
-2.169606808 -0.510143429  0.004249012  0.529962593  2.711585811

Residual standard error: 52.36185
Degrees of freedom: 180 total; 170 residual


Questions: Why are standard errors different between the gls() and lmer() fits? Is the model formulation different in some way that I am missing? Why is the estimated correlation matrix within subjects different between the two? Why might one prefer lmer over gls or even GEE when fitting this particular model?

Here are some other Stack Exchange posts that were a little relevant, but didn't exactly answer my question:

Which mixed model is closest to GLS with unstructured covariance? Random intercept + slopes or random slopes only?

Which mixed model in lme4 best mimics the unstructured covariance in GLS?

Linear mixed model with unstructured repeated measures

We have been working a lot in the past 2 years to fit MMRMs in R. Especially the unstructured covariance model is challenging. I would recommend to give the new R package mmrm a try. See https://openpharma.github.io/mmrm/main/ to get started. With that, in your example you can use:
library(lme4)
sleep_data$Day_Factor <- as.factor(sleep_data$Days)