Repeated measures within participant

Question

I'm trying to learn linear mixed effects models and how to estimate them using the R package lme4 and I am confused about some aspects.

I have a dataset where a number of participants completed 6 trials for each condition in a 2 × 2 factorial design (factor1 and factor2). I want to model intercept, slope, and interaction between factors as random effects. Following the software paper by Bates et al. (2015) for the lme4 package and using their example as a starting point, I came up with this model:

response ~ 1 + factor1 * factor2 + (1 + factor1 * factor2 | participant)

What I'm unclear about is how to account for the repeated-measures character of the responses; within participant, responses for different conditions may be correlated across trials. There is no trial variable, just 6 different observations for each condition and participant. E.g. for participant 1:

Question: Do I have to add something to my model to account for the correlation within participant between conditions across trials, and if yes, what?

Adding to my confusion, the conditional model in Eq. (2), i.e. within-participant in my application, $$ (\mathcal{Y} | \mathcal{B} = b) \mathrel{\sim} \mathcal{N} \left( X \beta + Z b, \sigma^2 W^{-1} \right) $$ includes a covariance matrix $\sigma^2 W^{-1}$, but the description says that $W$ is diagonal, i.e. it assumes observations within participant are uncorrelated.

In principle, the example analysis of sleepstudy in the paper has a similar problem: The reaction time measurements form a time series across days. However, in that case there is only one time series per subject, so there is no way to estimate these correlations (at least not without further assumptions like an AR process) – so I guess they don't need to be modeled?

Answer 1

Question: Do I have to add something to my model to account for the correlation within participant between conditions across trials, and if yes, what?

No, you do not need to. Random intercepts or coefficients by group are equivalent to compound symmetry of the residual variance-covariance matrix by group. Therefore, a correlation structure of the residuals is often unnecessary in presence of random effects. Otherwise, the algorithm will be confused about how to split the correlation pattern between random effects and correlated residuals, a dilemma analogous to perfect multicollinearity among predictors.

More important than the possible residual correlation is the limited responses. Your responses appear to be a 1-9 Likert scale. The residuals might be spuriously correlated due to this limited range if you apply linear models. You should consider the package ordinal for random effects in cumulative link models. See Liddell, T. M., & Kruschke, J. K. (2018). Analyzing ordinal data with metric models: What could possibly go wrong? Journal of Experimental Social Psychology, 79, 328–348. https://doi.org/10.1016/j.jesp.2018.08.009.

By (1 + factor1 * factor2 | participant), you will have 4 + 3 + 2 + 1 = 10 extra parameters to estimate for the error variance component (a variance-covariance matrix) of random effects. Because all random effects are grouped by the participant, the error term (which includes random intercepts, random coefficients, and residuals) among 24 measurements within each participants are already correlated, which results in an intra-class correlation (ICC) by $\sigma^2_\text{random} / (\sigma^2_\text{random} + \sigma^2_\text{residual})$. See https://en.wikipedia.org/wiki/Intraclass_correlation. Here, the denominator should be the sum of the 4 x 4 matrix (sum of 16 elements), whereas the numerator should be the sum of a 3 x 3 matrix less the row and column of the residuals (sum of 9 elements). We can consider that residuals are random effects grouped by each observation. Nevertheless, this unrestricted 10-unique-element matrix (a general positive-definite symmetric matrix) might be too complex. If you see any component to be close to zero or have a large standard error, you can restrict this variance component matrix to follow a simpler pattern (e.g., diagonal for only 4 parameters).

In your study, each participant experienced four conditions out of two two-level factors, levels A and B in factor1 and levels 1 and 2 in factor2. Under each condition, a participant was measured six times. Therefore, each participants have 24 measurements. Just like if someone measures my height six times, these measurements can be independent although they are all of the same person.

Theoretically, the trials will be independent if their error term (the difference between predicted and observed response) is not related to each other. There is no way to prove it because the error term cannot be directly observed, but this error independence needs to be justified logically. Does one overestimation in one trial leads on average an overestimation or underestimation in the next trial, after unobserved individual effects are controlled for by random effects?

Practically, modeling correlation across trials within the same participant requires prespecified patterns of the residuals. Are the trials done by different surveyors, scales, or equipment, so they can be grouped by Surveyor 1 ... Surveyor 6? Are the trials done sequentially in a specific order or at different time like Day 1, Day 2, ..., Day 6? Did the trials take place at specific locations that can be geographically related by their coordinates? If the study did not record any information regarding how different trials should be associated with each other, a researcher cannot retrieve the residual-correlation pattern by trial even if one suspects any.

Assuming the order of or association pattern among repeated measurements are available as a factor or continuous variable time, you will need to decide whether it should be coded as 1–6, 1–24, or some other values where the distance between consecutive trials differ. AR(1) might be a good choice, but other alternative correlation structures require considering. For the trial index time, we need its fixed effect as either a continuous or categorical variable.

The package {nlme} allows modeling random intercepts and coefficients along with AR process of the residual term, but it may not be necessary or significant for certain data. In the sleepstudy data for example, if the variation by participant is controlled for, the residuals that represent the difference between predicted and observed reaction time may not present enough temporal correlation to necessitate error-correlation consideration. You can also use the package {glmmTMB} for some special error-correlation patterns.

Answer 2

In addition to my general answer on the distinction between random effects and correlated residuals, a practical example follows using the classical reaction-time data and the R package {nlme}. For a study that addresses random coefficients and correlated residuals simultaneously, see a Master's thesis Qilong Yi (2000) Random effects and AR(1) models in longitudinal data analysis https://tspace.library.utoronto.ca/bitstream/1807/15042/1/MQ49731.pdf. For sample codes of fitting both random coefficients and AR(1) in {glmmTMB}: ar1 covariance with nested random effects https://github.com/glmmTMB/glmmTMB/issues/329 need glmmTMB((1|id)+ ar1(0 + factor(time)|id), disp = ~ 1).

As a primer, maximum likelihood (ML) and restricted maximum likelihood (REML) estimators make different estimates of error variances. As the following example shows in generalized least squares (GLS), ML underestimates the error variance, whereas REML estimate is identical to OLS, which produces unbiased estimates of error variances. Note that glm() uses OLS instead of ML. When comparing models, we must compare those estimated by ML against others estimated by ML and those estimated by REML against others estimated by REML with the same fixed-effect component (same formula other than the random effects, residual correlation, or variance function). This is because goodness-of-fit statistics including log likelihoods (thus AIC and BIC) are calculated in different ways under ML and REML. By default, gls() and lme() in {nlme} use REML, whereas {glmmTMB} uses ML.

Because of added complexity, fitting linear models with mixed effects can occasionally run into numerical problems. When using {nlme}, it is important to use its function intervals() to assess whether estimates and confidence intervals are at parameter boundaries. When using {glmmTMB}, use glmmTMB()$sdr$pdHess and diagnose(glmmTMB()) to assess fitting problems.

data("sleepstudy", package = "lme4")
library(nlme)
library(glmmTMB)

# OLS
summary(Model1 <- lm(
  Reaction ~ Days, data = sleepstudy))
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  251.405      6.610  38.033  < 2e-16 ***
Days          10.467      1.238   8.454 9.89e-15 ***
Residual standard error: 47.71 on 178 degrees of freedom
Multiple R-squared:  0.2865,    Adjusted R-squared:  0.2825 
F-statistic: 71.46 on 1 and 178 DF,  p-value: 9.894e-15"
logLik(Model1)
"'log Lik.' -950.1465 (df=3)"
AIC(Model1)
"1906.293"
sigma(Model1)
"47.71472"
coef(Model1)
"(Intercept)        Days 
  251.40510    10.46729"

# glm()
summary(Model2 <- glm(
  Reaction ~ Days, data = sleepstudy))
"            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  251.405      6.610  38.033  < 2e-16 ***
Days          10.467      1.238   8.454 9.89e-15 ***
(Dispersion parameter for gaussian family taken to be 2276.694)
    Null deviance: 567954  on 179  degrees of freedom
Residual deviance: 405252  on 178  degrees of freedom
AIC: 1906.3"
logLik(Model2) == logLik(Model1) # TRUE
sigma(Model2) == sigma(Model1) # TRUE
AIC(Model2) == AIC(Model1) # TRUE
coef(Model2) == coef(Model1) # TRUE        TRUE
"glm() gives OLS estimates of error variance when family = gaussian"

# gls() REML
summary(Model3 <- gls(
  Reaction ~ Days, data = sleepstudy))
"Generalized least squares fit by REML
       AIC      BIC    logLik
  1899.664 1909.209 -946.8318
                Value Std.Error  t-value p-value
(Intercept) 251.40510  6.610154 38.03317       0
Days         10.46729  1.238195  8.45366       0
AIC differs because fitted by REML, default method"
intervals(Model3)
"                lower      est.     upper
(Intercept) 238.360753 251.40510 264.44946
Days          8.023855  10.46729  12.91072
 Residual standard error:
   lower     est.    upper 
43.23105 47.71472 53.24421"
logLik(Model3) == logLik(Model1) # FALSE
logLik(Model3) - logLik(Model1) # 3.314696 REML vs. OLS, not really comparable
sigma(Model3) == sigma(Model1) # FALSE
sigma(Model3) - sigma(Model1) # -7.105427e-15 within floating-point precision 
AIC(Model3) == AIC(Model1) # FALSE
AIC(Model3) - AIC(Model1) # -6.629393 REML vs. OLS
coef(Model3) == coef(Model1) # FALSE        TRUE
coef(Model3) - coef(Model1) # 2.842171e-14 0.000000e+00

# gls() ML
summary(Model4 <- gls(
  Reaction ~ Days, data = sleepstudy, method = "ML"))
"       AIC      BIC    logLik
  1906.293 1915.872 -950.1465
                Value Std.Error  t-value p-value
(Intercept) 251.40510  6.610154 38.03317       0
Days         10.46729  1.238195  8.45366       0"
intervals(Model4)
"                 lower      est.     upper
(Intercept) 238.360753 251.40510 264.44946
Days          8.023855  10.46729  12.91072
 Residual standard error:
   lower     est.    upper 
43.01267 47.44890 52.91347
REML/ML affects error variance: ML reports smaller sigma and 95% CI than OLS
Because OLS gives unbiased variance sum/(n - 1) whereas ML sum/n"
logLik(Model4) == logLik(Model1) # FALSE
logLik(Model4) - logLik(Model1) # -1.136868e-13 floating-point precision
sigma(Model4) == sigma(Model1) # FALSE
sigma(Model4) - sigma(Model1) # -0.2658222 due to ML vs. OLS error variance
AIC(Model4) == AIC(Model1) # FALSE
AIC(Model4) - AIC(Model1) # 2.273737e-13 floating-point precision
coef(Model4) == coef(Model1) # FALSE        TRUE
coef(Model4) - coef(Model1) # 2.842171e-14 0.000000e+00
"gls(method = 'ML') gives true ML estimate in contrast to glm()
gls(method = 'ML') gls(method = 'REML') give identical coef and SE
  but differ in sigma: sigma(ML) > sigma(REML) = sigma(glm) = sigm(lm)
  and likelihood value: logLik(REML) != logLik(ML) = logLik(glm) = logLik(lm)"

A Gaussian linear model Model5 with compound-symmetry correlation structure is equivalent to Model6 with random intercepts. The two models have identical restricted loglikelihood and coefficient estimates. Their standard errors and confidence intervals differ by only a tiny bit. In contrast, the residual standard error estimates are widely apart because mixed-effects models separate random effects from residuals. They both fit better than Model3 by lower AIC and BIC and higher logLik, showing that incorporating repeated measurements within subjects is beneficial.

# gls() corCompSymm
summary(Model5 <- gls(
  Reaction ~ Days, data = sleepstudy,
  correlation = corCompSymm(form = ~ 1 | Subject)))
"       AIC      BIC    logLik
  1794.465 1807.192 -893.2325
Correlation Structure: Compound symmetry
 Formula: ~1 | Subject 
 Parameter estimate(s):
      Rho 
0.5893103 
                Value Std.Error  t-value p-value
(Intercept) 251.40510  9.746736 25.79378       0
Days         10.46729  0.804221 13.01543       0
Residual standard error: 48.3595 

Here restricted AIC < Model3 1899.664, residual SE larger"
intervals(Model5)
"                 lower      est.     upper
(Intercept) 232.171083 251.40510 270.63913
Days          8.880251  10.46729  12.05432
 Correlation structure:
        lower      est.     upper
Rho 0.3960776 0.5893103 0.7510617
 Residual standard error:
   lower     est.    upper 
38.96263 48.35950 60.02268"

# lme() random intercept
summary(Model6 <- lme(
  Reaction ~ Days, data = sleepstudy,
  random = ~ 1 | Subject))
"       AIC      BIC    logLik
  1794.465 1807.192 -893.2325
Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:    37.12383 30.99123
Fixed effects:  Reaction ~ Days 
                Value Std.Error  DF  t-value p-value
(Intercept) 251.40510  9.746716 161 25.79383       0
Days         10.46729  0.804221 161 13.01543       0"
intervals(Model6)
"                 lower      est.     upper
(Intercept) 232.157211 251.40510 270.65300
Days          8.879103  10.46729  12.05547
 Random Effects:
  Level: Subject 
                   lower     est.    upper
sd((Intercept)) 25.90982 37.12383 53.19137
 Within-group standard error:
   lower     est.    upper 
27.78454 30.99123 34.56802"
logLik(Model6) == logLik(Model5) # FALSE
logLik(Model6) - logLik(Model5) # 1.032276e-10 floating-point precision
sigma(Model6) == sigma(Model5) # FALSE
sigma(Model6) - sigma(Model5) # -17.36827 difference components of sigma
AIC(Model6) == AIC(Model5) # FALSE
AIC(Model6) - AIC(Model5) # -2.064553e-10 floating-point precision
summary(Model6)$tTable[ , 1] == coef(Model5)
summary(Model6)$tTable[ , 1] - coef(Model5) # 1.136868e-13 5.329071e-15

Adding random effects of both the intercept and the slope further improves model fit. But allowing these two random effects to correlate appears unnecessary, as the estimated correlation coefficient between these random terms is very small, with a confidence interval encompassing zero and worse information criteria.

# lme() uncorrelated random intercept + slope
summary(Model7 <- lme(
  Reaction ~ Days, data = sleepstudy |> transform(ID = Subject),
  random = list(~ 1 | Subject, ~ 0 + Days | ID)))
"       AIC      BIC    logLik
  1753.669 1769.578 -871.8346
Random effects:
 Formula: ~1 | Subject
        (Intercept)
StdDev:    25.05133
 Formula: ~0 + Days | ID %in% Subject
            Days Residual
StdDev: 5.988172 25.56529

Fixed effects:  Reaction ~ Days 
                Value Std.Error  DF  t-value p-value
(Intercept) 251.40510  6.885381 161 36.51288       0
Days         10.46729  1.559566 161  6.71167       0"
intervals(Model7)
"                 lower      est.     upper
(Intercept) 237.807798 251.40510 265.00241
Days          7.387443  10.46729  13.54713
  Level: Subject 
                   lower     est.    upper
sd((Intercept)) 16.08536 25.05133 39.01492
  Level: ID 
            lower     est.    upper
sd(Days) 4.025233 5.988172 8.908353
 Within-group standard error:
   lower     est.    upper 
22.79178 25.56529 28.67630

If Subject is not copied and renamed, the print method of intervals()
will mistake these two components and report
 Random Effects:
  Level: Subject 
            lower     est.    upper
sd(Days) 16.08536 25.05133 39.01492
  Level: Subject 
            lower     est.    upper
sd(Days) 16.08536 25.05133 39.01492"

# lme() correlated random intercept + slope
summary(Model8 <- lme(
  Reaction ~ Days, data = sleepstudy,
  random = ~ 1 + Days | Subject))
"A worse fit than random intercept + AR1
       AIC      BIC    logLik
  1755.628 1774.719 -871.8141
Random effects:
 Formula: ~1 + Days | Subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 24.740241 (Intr)
Days         5.922103 0.066 
Residual    25.591843       
Fixed effects:  Reaction ~ Days 
                Value Std.Error  DF  t-value p-value
(Intercept) 251.40510  6.824516 161 36.83853       0
Days         10.46729  1.545783 161  6.77151       0"
intervals(Model8)
"                 lower      est.     upper
(Intercept) 237.927995 251.40510 264.88221
Days          7.414662  10.46729  13.51991
                           lower        est.      upper
sd((Intercept))       15.5175354 24.74024072 39.4443767
sd(Days)               3.9041696  5.92210282  8.9830375
cor((Intercept),Days) -0.5760831  0.06556383  0.6572156
 Within-group standard error:
   lower     est.    upper 
22.79706 25.59184 28.72925"

Using autoregressive residual correlation fits the model substantially better in both AIC and BIC than those of compound symmetry residuals or random effects. Combining autoregressive residuals and random intercepts improves AIC marginally but inflates BIC.

# gls() AR1
summary(Model9 <- gls(
  Reaction ~ Days, data = sleepstudy,
  correlation = corAR1(form = ~ Days | Subject)))
"       AIC      BIC   logLik
  1747.206 1759.933 -869.603
Correlation Structure: AR(1)
 Formula: ~Days | Subject      # here Days correct for missing days if any
 Parameter estimate(s):
     Phi 
0.799953
                Value Std.Error   t-value p-value
(Intercept) 253.73770 11.250640 22.553180       0
Days         10.46673  1.699413  6.159028       0
Residual standard error: 49.51497

Here restricted AIC < Model3 1794.465 < Model5 1899.664
residual/coef SE larger
corAR1 is better than corCompSymm for the data"
intervals(Model9)
"                 lower      est.     upper
(Intercept) 231.535905 253.73770 275.93950
Days          7.113145  10.46673  13.82032
 Correlation structure:
        lower     est.     upper
Phi 0.6995489 0.799953 0.8693829
 Residual standard error:
   lower     est.    upper 
40.59756 49.51497 60.39113"
varcomp_vcov(Model9)
"             cor_params     sigma_sq
cor_params  0.001717378     17.05477
sigma_sq   17.054771342 236905.06677"

# lme() AR1 + random intercept
summary(Model10 <- lme(
  Reaction ~ Days, data = sleepstudy,
  random = ~ 1 | Subject, 
  correlation = corAR1(form = ~ Days | Subject)))
"       AIC    BIC    logLik
  1746.191 1762.1 -868.0955
Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:    31.71358 37.63734
Correlation Structure: AR(1)
 Formula: ~Days | Subject    
 Parameter estimate(s):
      Phi 
0.6531333 
Fixed effects:  Reaction ~ Days 
                Value Std.Error  DF   t-value p-value
(Intercept) 252.84144  10.94093 161 23.109688       0
Days         10.46687   1.34408 161  7.787388       0

Here restricted AIC < Model7 1747.206
residual SE larger than lme() random intercepts only
re + corAR1 is better than corAR1 alone for the data"
intervals(Model10)
" Fixed effects:
                 lower      est.     upper
(Intercept) 231.235210 252.84144 274.44768
Days          7.812573  10.46687  13.12117
 Random Effects:
  Level: Subject 
                   lower     est.    upper
sd((Intercept)) 18.49463 31.71358 54.38072
 Correlation structure:
        lower      est.     upper
Phi 0.4341414 0.6531333 0.7992378
 Within-group standard error:
   lower     est.    upper 
29.19232 37.63734 48.52541 "
varcomp_vcov(Model10)
"                             Tau.Subject.var((Intercept))    cor_params
Tau.Subject.var((Intercept))                 311496.61500 -19.066132731
cor_params                                      -19.06613   0.008051599
sigma_sq                                     -77449.46722  29.062277325
                                 sigma_sq
Tau.Subject.var((Intercept)) -77449.46722
cor_params                       29.06228
sigma_sq                     129827.95559"

Adding both random intercepts and random slopes in presence of autoregressive errors further improves model fits when these two components are uncorrelated. Allowing the two random terms to correlate, however, confuses the nlminb algorithm. Switching the optimizer to optim converges but shows that correlation between the two random terms essentially cannot be identified.

# lme() AR1 + uncorrelated random intercept/slope
summary(Model11 <- lme(
  Reaction ~ Days, data = sleepstudy |> transform(ID = Subject),
  random = list(~ 1 | Subject, ~ 0 + Days | ID), 
  correlation = corAR1(form = ~ Days | Subject)))
"       AIC      BIC    logLik
  1738.196 1757.287 -863.0981
Random effects:
 Formula: ~1 | Subject
        (Intercept)
StdDev:    20.15631
 Formula: ~0 + Days | Subject %in% Subject
            Days Residual
StdDev: 5.666592   29.557
Correlation Structure: AR(1)
 Formula: ~Days | Subject/Subject 
 Parameter estimate(s):
      Phi 
0.4544552 
Fixed effects:  Reaction ~ Days 
                Value Std.Error  DF  t-value p-value
(Intercept) 252.15508  7.336304 161 34.37086       0
Days         10.46704  1.667154 161  6.27839       0
Warning message:
In lme.formula(Reaction ~ Days, data = sleepstudy, random = list(~1 |  :
  cannot use smaller level of grouping for 'correlation' than for 'random'. 
  Replacing the former with the latter."
intervals(Model11) # looks okay
"                 lower      est.     upper
(Intercept) 237.667293 252.15508 266.64288
Days          7.174727  10.46704  13.75934
  Level: Subject 
                   lower     est.    upper
sd((Intercept)) 9.386327 20.15631 43.28389
  Level: ID 
            lower     est.   upper
sd(Days) 3.516849 5.666592 9.13041
 Correlation structure:
        lower      est.     upper
Phi 0.2105079 0.4544552 0.6451227
 Within-group standard error:
   lower     est.    upper 
24.18822 29.55700 36.11744"

# lme() AR1 + correlated random intercept/slope
summary(Model12 <- lme(
  Reaction ~ Days, data = sleepstudy,
  random = ~ 1 + Days | Subject, 
  correlation = corAR1(form = ~ Days | Subject)))
"Error in lme.formula(Reaction ~ Days, data = sleepstudy, random = ~1 +  : 
  nlminb problem, convergence error code = 1
  message = iteration limit reached without convergence (10)"

summary(Model12 <- lme(
  Reaction ~ Days, data = sleepstudy,
  random = ~ 1 + Days | Subject, 
  correlation = corAR1(form = ~ Days | Subject),
  control = lmeControl(maxIter = 50000, msMaxIter = 1000, msMaxEval = 2000)))
"Error in lme.formula(Reaction ~ Days, data = sleepstudy, random = ~1 +  : 
  nlminb problem, convergence error code = 1
  message = singular convergence (7)"

summary(Model12 <- lme(
  Reaction ~ Days, data = sleepstudy,
  random = ~ 1 + Days | Subject, 
  correlation = corAR1(form = ~ Days | Subject), 
  control = lmeControl(opt = "optim")))
"       AIC     BIC    logLik
  1738.187 1760.46 -862.0935
Random effects:
 Formula: ~1 + Days | Subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 14.879196 (Intr)
Days         4.759861 0.897 
Residual    30.500694       
Correlation Structure: AR(1)
 Formula: ~Days | Subject 
 Parameter estimate(s):
      Phi 
0.4870368 
Fixed effects:  Reaction ~ Days 
                Value Std.Error  DF  t-value p-value
(Intercept) 252.24315  6.851137 161 36.81770       0
Days         10.46701  1.532133 161  6.83166       0"
intervals(Model12)
"                 lower      est.     upper
(Intercept) 238.713471 252.24315 265.77283
Days          7.441343  10.46701  13.49268
 Random Effects: Level: Subject 
                          lower       est.      upper
sd((Intercept))        5.494663 14.8791959 40.2919116
sd(Days)               2.644655  4.7598609  8.5668177
cor((Intercept),Days) -0.997741  0.8967291  0.9999933
 Correlation structure:
        lower      est.     upper
Phi 0.2791926 0.4870368 0.6512869
 Within-group standard error:
   lower     est.    upper 
25.45483 30.50069 36.54679 

cor((Intercept),Days) essentially [-1, 1], should be removed"

In summary, having both correlated random effects and correlated residuals may be too complex for the model fitting algorithms to converge. I have also run into occasions where widely apart ranges of predictors also result in noncoverages, and simple scaling of some predictors will resolve the optimization issue. We need to use converged models with estimates and confidence intervals away from theoretical boundaries and select parsimonious models based on AIC and BIC.