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 in lme4 best mimics the unstructured covariance in GLS?
