In this chapter, we are going to combine some of the concepts that we covered in Section 2, on mixed-factorial designs, with Section 3, on times series data. In combining these sections, we can see how we would analyze data when we have two time series that come from same person being measured at different times.

In rehabilitation, for instance we might be interested in the rate at which a person with Parkinson’s disease acquires a skill when they are on or off dopamine replacement medication. Alternately, might be interested in hemibody differences between the left and right arms, or left and right legs, over time.

As before however, we will work with our hypothetical dataset in which younger and older adults walked in a virtual environment on four different trials under three different conditions (arbitrary conditions A, B, and C).

## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'tibble' was built under R version 4.3.3
## Warning: package 'purrr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.3

1. Modeling Changes in Speed across 4 Trials, in 2 Different Conditions

To account for the fact that participants were tested under each of the three different conditions, we need a random slope for time within condition (rather than only within subject), since each subject’s trajectory across trials can differ from one condition to the next.

In these models, we will mean-center time to create a new variable time.c (\(time_i - \overline{time}\)). Centering time around it’s mean changes the interpretation of the intercept from being the predicted value on Trial 0 (which does exist) to the predicted value on average over time. (Note that the “average trial” also doesn’t exist, but has a more useful interpretation.)

# Fixed slope random intercepts model ---- 
mod01<-lmer(speed~ 
                  # Fixed-effects 
                  1+time.c*condition+time.c.sq*condition+ 
                  # Random-effects 
                  (1|subID)+(1+time.c+time.c.sq|condition:subID),
                  data=DATA, REML=TRUE,
                control=lmerControl(optimizer="Nelder_Mead",
                                    optCtrl=list(maxfun=2e5)))
## boundary (singular) fit: see help('isSingular')
anova(mod01)
summary(mod01)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: speed ~ 1 + time.c * condition + time.c.sq * condition + (1 |  
##     subID) + (1 + time.c + time.c.sq | condition:subID)
##    Data: DATA
## Control: lmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: -423.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0476 -0.4643 -0.0816  0.3121  5.2157 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev. Corr       
##  condition:subID (Intercept) 0.001517 0.03894             
##                  time.c      0.001779 0.04217   0.16      
##                  time.c.sq   0.001532 0.03914  -0.15 -1.00
##  subID           (Intercept) 0.019438 0.13942             
##  Residual                    0.012659 0.11251             
## Number of obs: 480, groups:  condition:subID, 120; subID, 40
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)            0.797500   0.026956  66.081478  29.586  < 2e-16 ***
## time.c                -0.088800   0.010381 130.524209  -8.554 2.78e-14 ***
## conditionB             0.110578   0.021939 100.774437   5.040 2.05e-06 ***
## conditionC             0.049875   0.021939 100.774437   2.273  0.02513 *  
## time.c.sq              0.036000   0.010836 142.598721   3.322  0.00113 ** 
## time.c:conditionB     -0.085825   0.014681 130.524209  -5.846 3.81e-08 ***
## time.c:conditionC     -0.032750   0.014681 130.524209  -2.231  0.02740 *  
## conditionB:time.c.sq   0.006188   0.015324 142.598723   0.404  0.68699    
## conditionC:time.c.sq   0.000500   0.015324 142.598723   0.033  0.97402    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) time.c cndtnB cndtnC tm.c.s tm.c:B tm.c:C cnB:..
## time.c       0.023                                                 
## conditionB  -0.407 -0.028                                          
## conditionC  -0.407 -0.028  0.500                                   
## time.c.sq   -0.358 -0.367  0.440  0.440                            
## tm.c:cndtnB -0.016 -0.707  0.040  0.020  0.259                     
## tm.c:cndtnC -0.016 -0.707  0.020  0.040  0.259  0.500              
## cndtnB:tm..  0.253  0.259 -0.622 -0.311 -0.707 -0.367 -0.183       
## cndtnC:tm..  0.253  0.259 -0.311 -0.622 -0.707 -0.183 -0.367  0.500
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

A Note (2024 revision): This Model is Already a Singular Fit

If you run the model above, you will see a “boundary (singular) fit” warning. Looking at VarCorr(mod01) shows why: the estimated correlation between the random linear slope and the random quadratic slope of time (within condition:subID) is -1.00 – an impossible-to-improve-on boundary value, meaning the model is trying to estimate more structure (a full 3x3 variance-covariance matrix, including three separate correlation parameters) than these data can actually support.

This is exactly the kind of situation discussed in Section 3, Part 5.3 of this repository: the “maximal” random-effects structure is not automatically the right choice just because it is more complete on paper. One straightforward alternative is to remove the correlations among the random slopes/intercept using the || operator in lme4:

# Zero-correlation alternative ----
mod01b<-lmer(speed~ 
                  1+time.c*condition+time.c.sq*condition+ 
                  (1|subID)+(1+time.c+time.c.sq||condition:subID),
                  data=DATA, REML=TRUE,
                control=lmerControl(optimizer="Nelder_Mead",
                                    optCtrl=list(maxfun=2e5)))
isSingular(mod01b)
## [1] FALSE
anova(mod01b)

This resolves the convergence warning (confirmed: isSingular(mod01b) returns FALSE) while reaching the same qualitative conclusions about which effects are significant. Removing the correlation parameters assumes that individual differences in the intercept, linear slope, and quadratic slope are independent of one another, which is a real (and testable, via model comparison) assumption – but it is often a very reasonable one to make when, as here, the fully-correlated model cannot even be estimated without hitting a boundary.