Mixed-Effects Models

For an overview of mixed modeling, see my chapter on “Mixed Modeling” in the following textbook that is available for free online:

Petersen, I. T. (2025). Fantasy football analytics: Statistics, prediction, and empiricism using R. University of Iowa Libraries. https://isaactpetersen.github.io/Fantasy-Football-Analytics-Textbook

1 Preamble

1.1 Load Libraries

Code

1.2 Load Data

The data for this example were simulated:

Code
mydata_long <- read.csv("./data/mydata_long.csv")

2 Research Questions

  1. What is the form of adolescents’ depressive symptom trajectories across ages 14 to 17?
  2. Do the trajectories differ between boys and girls?
  3. Does how much sleep the adolescent receives predict their depressive symptoms over time?
  4. Does the association between sleep and depressive symptoms differ between boys and girls?
  5. What is the association between sleep and depressive symptoms at the within-individual level versus between-individual level?

3 Pre-Model Computation

It can be helpful to center the age/time variable so that the intercept in a growth curve model has meaning. For instance, we can subtract the youngest participant age to set the intercepts to be the earliest age in the sample.

Code
mydata_long$sexFactor <- factor(mydata_long$sex) # create a factor variable for sex

mydata_long$ageCentered <- mydata_long$age - min(mydata_long$age, na.rm = TRUE) # shift age so that intercepts represent the earliest age in the sample

mydata_long$ageCenteredSquared <- mydata_long$ageCentered ^ 2

4 Linear Growth Curve Model

4.1 Fit the Model

In estimating linear growth curves of depression symptoms, this model fits:

  1. random intercepts and slopes for subjects,
  2. sex as a time-invariant fixed-effect predictor of the intercepts and slopes, and
  3. sleep as a time-varying fixed-effect predictor of depression, with effects allowed to vary over age and by sex, including two‐ and three‐way interactions with age and sex.
Code
linearGCM <- lmerTest::lmer(
  depression ~ sexFactor + ageCentered + sexFactor:ageCentered + 
    sleep + sleep:ageCentered + sexFactor:sleep + sexFactor:sleep:ageCentered + 
    (1 + ageCentered | ID),
  data = mydata_long,
  na.action = na.exclude)

4.2 Model Summary

Code
summary(linearGCM)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: depression ~ sexFactor + ageCentered + sexFactor:ageCentered +  
    sleep + sleep:ageCentered + sexFactor:sleep + sexFactor:sleep:ageCentered +  
    (1 + ageCentered | ID)
   Data: mydata_long

REML criterion at convergence: 4763.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6405 -0.4712 -0.0088  0.4393  2.5703 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 ID       (Intercept) 14.144   3.761         
          ageCentered  2.342   1.530    0.60 
 Residual              3.029   1.740         
Number of obs: 931, groups:  ID, 250

Fixed effects:
                                 Estimate Std. Error        df t value
(Intercept)                      17.52934    1.36182 532.76805  12.872
sexFactormale                    -1.88370    2.02693 535.62207  -0.929
ageCentered                       1.79858    0.71361 476.29718   2.520
sleep                            -0.86981    0.16355 505.61680  -5.318
sexFactormale:ageCentered        -2.52264    1.08549 493.04238  -2.324
ageCentered:sleep                -0.10633    0.08654 466.35616  -1.229
sexFactormale:sleep              -0.04683    0.24585 508.01056  -0.191
sexFactormale:ageCentered:sleep   0.10563    0.13247 482.46640   0.797
                                            Pr(>|t|)    
(Intercept)                     < 0.0000000000000002 ***
sexFactormale                                 0.3531    
ageCentered                                   0.0120 *  
sleep                                    0.000000157 ***
sexFactormale:ageCentered                     0.0205 *  
ageCentered:sleep                             0.2198    
sexFactormale:sleep                           0.8490    
sexFactormale:ageCentered:sleep               0.4256    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) sxFctr agCntr sleep  sxFc:C agCnt: sxFct:
sexFactorml -0.672                                          
ageCentered -0.536  0.360                                   
sleep       -0.968  0.650  0.573                            
sxFctrml:gC  0.353 -0.529 -0.657 -0.377                     
agCntrd:slp  0.566 -0.380 -0.978 -0.584  0.643              
sxFctrml:sl  0.644 -0.967 -0.381 -0.665  0.566  0.389       
sxFctrml:C: -0.370  0.564  0.639  0.382 -0.979 -0.653 -0.583
Code
print(effectsize::standardize_parameters(
  linearGCM,
  method = "refit"),
  digits = 2)
# Standardization method: refit

Parameter                                | Std. Coef. |         95% CI
----------------------------------------------------------------------
(Intercept)                              |       0.31 | [ 0.17,  0.44]
sexFactor [male]                         |      -0.71 | [-0.91, -0.50]
ageCentered                              |       0.16 | [ 0.11,  0.21]
sleep                                    |      -0.18 | [-0.22, -0.13]
sexFactor [male] × ageCentered           |      -0.28 | [-0.36, -0.21]
ageCentered × sleep                      |      -0.02 | [-0.05,  0.01]
sexFactor [male] × sleep                 |       0.02 | [-0.05,  0.09]
(sexFactor [male] × ageCentered) × sleep |       0.02 | [-0.03,  0.07]
Code
performance::r2(linearGCM)
# R2 for Mixed Models

  Conditional R2: 0.928
     Marginal R2: 0.175

4.3 Plots

Code
newData <- expand.grid(
  sexFactor = c("female", "male"),
  age = c(
    min(mydata_long$age, na.rm = TRUE),
    max(mydata_long$age, na.rm = TRUE)),
  sleepFactor = c("low sleep", "average sleep", "high sleep")
)

newData$ageCentered <- newData$age - min(newData$age)

newData$sleep <- NA
newData$sleep[which(newData$sexFactor == "female" & newData$age == 14 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age == 14 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age == 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age == 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

newData$sleep[which(newData$sexFactor == "female" & newData$age == 14 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age == 14 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age == 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age == 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

newData$sleep[which(newData$sexFactor == "female" & newData$age == 14 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age == 14 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age == 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age == 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

newData$modelImpliedDepression <- predict( # predict.merMod
  linearGCM,
  newdata = newData,
  re.form = NA # include no random effects
)

4.3.1 Prototypical Growth Curve By Sex

Code
ggplot(
  data = newData |> filter(sleepFactor == "average sleep"),
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    color = sexFactor)) +
  geom_line() + # fit lines to the model-implied points
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex

4.3.2 Individuals’ Growth Curves

Code
mydata_long$modelImpliedDepression <- predict( # predict.merMod
  linearGCM,
  newdata = mydata_long,
  re.form = NULL # include all random effects
)

ggplot(
  data = mydata_long,
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    group = ID,
    color = sexFactor)) +
  geom_smooth( # fit smooth to each participant's points
    method = "lm",
    se = FALSE,
    linewidth = 0.5
  ) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Individuals’ Model-Implied Growth Curves of Depression Symptoms

Individuals’ Model-Implied Growth Curves of Depression Symptoms

4.3.3 Individuals’ Trajectories Overlaid with Prototypical Trajectory By Sex

Code
ggplot(
  data = mydata_long,
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    group = ID)) +
  geom_smooth( # fit smooth to each participant's points
    method = "lm",
    se = FALSE,
    linewidth = 0.5,
    color = "gray",
    alpha = 0.30
  ) +
  geom_line( # fit lines to the model-implied points
    data = newData |> filter(sleepFactor == "average sleep"),
    mapping = aes(
      x = age,
      y = modelImpliedDepression,
      group = sexFactor,
      color = sexFactor),
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Individuals’ Model-Implied Growth Curves of Depression Symptoms Overlaid With Prototypical Growth Curves By Sex

Individuals’ Model-Implied Growth Curves of Depression Symptoms Overlaid With Prototypical Growth Curves By Sex

4.3.4 Prototypical Trajectory By Sex and Sleep

Code
ggplot(
  data = newData |> filter(sleepFactor != "average sleep"),
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    color = sexFactor,
    linetype = sleepFactor)) +
  geom_line() + # fit lines to the model-implied points
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex",
    linetype = "Sleep"
  ) +
  theme_classic()

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex and Sleep

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex and Sleep

5 Polynomial (Quadratic) Growth Curve Model

5.1 Fit the Model

In estimating polynomial (in this case, quadratic) growth curves of depression symptoms, this model fits:

  1. random intercepts and random linear slopes and fixed quadratic slopes for subjects,
  2. sex as a time-invariant fixed-effect predictor of the intercepts and slopes, and
  3. sleep as a time-varying fixed-effect predictor of depression, with effects allowed to vary over age and by sex, including two‐ and three‐way interactions with age and sex.
Code
quadraticGCM <- lmerTest::lmer(
  depression ~ sexFactor + ageCentered + ageCenteredSquared + 
    sexFactor:ageCentered + sexFactor:ageCenteredSquared + sleep + 
    sleep:ageCentered + sleep:ageCenteredSquared + sexFactor:sleep + 
    sexFactor:sleep:ageCentered + sexFactor:sleep:ageCenteredSquared + 
    (1 + ageCentered | ID), # model failed to converge with random quadratic slopes (ageCenteredSquared)
  data = mydata_long,
  na.action = na.exclude)

5.2 Model Summary

Code
summary(quadraticGCM)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
depression ~ sexFactor + ageCentered + ageCenteredSquared + sexFactor:ageCentered +  
    sexFactor:ageCenteredSquared + sleep + sleep:ageCentered +  
    sleep:ageCenteredSquared + sexFactor:sleep + sexFactor:sleep:ageCentered +  
    sexFactor:sleep:ageCenteredSquared + (1 + ageCentered | ID)
   Data: mydata_long

REML criterion at convergence: 4644.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6202 -0.4282 -0.0301  0.4115  3.1856 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 ID       (Intercept) 14.808   3.848         
          ageCentered  2.522   1.588    0.53 
 Residual              2.237   1.496         
Number of obs: 931, groups:  ID, 250

Fixed effects:
                                        Estimate Std. Error        df t value
(Intercept)                             17.21489    1.39365 606.93322  12.352
sexFactormale                           -0.87170    2.06548 605.43014  -0.422
ageCentered                              1.96233    1.69142 578.56873   1.160
ageCenteredSquared                      -0.09582    0.52663 488.62087  -0.182
sleep                                   -0.76115    0.16713 574.77856  -4.554
sexFactormale:ageCentered               -7.72049    2.53109 579.67774  -3.050
sexFactormale:ageCenteredSquared         1.83789    0.80054 497.69539   2.296
ageCentered:sleep                       -0.35167    0.20756 571.64643  -1.694
ageCenteredSquared:sleep                 0.08849    0.06493 488.72438   1.363
sexFactormale:sleep                     -0.17418    0.25079 572.16459  -0.695
sexFactormale:ageCentered:sleep          0.75518    0.31258 576.74511   2.416
sexFactormale:ageCenteredSquared:sleep  -0.22759    0.09916 500.05445  -2.295
                                                   Pr(>|t|)    
(Intercept)                            < 0.0000000000000002 ***
sexFactormale                                       0.67315    
ageCentered                                         0.24646    
ageCenteredSquared                                  0.85569    
sleep                                            0.00000642 ***
sexFactormale:ageCentered                           0.00239 ** 
sexFactormale:ageCenteredSquared                    0.02210 *  
ageCentered:sleep                                   0.09075 .  
ageCenteredSquared:sleep                            0.17356    
sexFactormale:sleep                                 0.48762    
sexFactormale:ageCentered:sleep                     0.01600 *  
sexFactormale:ageCenteredSquared:sleep              0.02213 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) sxFctr agCntr agCntS sleep  sxFc:C sxF:CS agCnt: agCnS:
sexFactorml -0.675                                                        
ageCentered -0.557  0.376                                                 
agCntrdSqrd  0.372 -0.251 -0.915                                          
sleep       -0.968  0.653  0.577 -0.377                                   
sxFctrml:gC  0.372 -0.544 -0.668  0.612 -0.386                            
sxFctrml:CS -0.245  0.358  0.602 -0.658  0.248 -0.913                     
agCntrd:slp  0.563 -0.380 -0.989  0.912 -0.580  0.661 -0.600              
agCntrdSqr: -0.366  0.247  0.906 -0.992  0.377 -0.606  0.652 -0.918       
sxFctrml:sl  0.645 -0.967 -0.385  0.251 -0.666  0.564 -0.363  0.387 -0.251
sxFctrml:C: -0.374  0.557  0.657 -0.605  0.385 -0.989  0.908 -0.664  0.609
sxFctrm:CS:  0.240 -0.359 -0.594  0.649 -0.247  0.905 -0.992  0.601 -0.655
            sxFct: sxF:C:
sexFactorml              
ageCentered              
agCntrdSqrd              
sleep                    
sxFctrml:gC              
sxFctrml:CS              
agCntrd:slp              
agCntrdSqr:              
sxFctrml:sl              
sxFctrml:C: -0.575       
sxFctrm:CS:  0.370 -0.915
Code
print(effectsize::standardize_parameters(
  quadraticGCM,
  method = "refit"),
  digits = 2)
# Standardization method: refit

Parameter                                       | Std. Coef. |         95% CI
-----------------------------------------------------------------------------
(Intercept)                                     |       0.31 | [ 0.17,  0.44]
sexFactor [male]                                |      -0.70 | [-0.91, -0.49]
ageCentered                                     |      -0.15 | [-0.23, -0.06]
ageCenteredSquared                              |       0.32 | [ 0.25,  0.39]
sleep                                           |      -0.17 | [-0.21, -0.13]
sexFactor [male] × ageCentered                  |      -0.28 | [-0.40, -0.16]
sexFactor [male] × ageCenteredSquared           |   5.95e-03 | [-0.10,  0.11]
ageCentered × sleep                             |      -0.07 | [-0.15,  0.01]
ageCenteredSquared × sleep                      |       0.05 | [-0.02,  0.13]
sexFactor [male] × sleep                        |       0.03 | [-0.03,  0.09]
(sexFactor [male] × ageCentered) × sleep        |       0.14 | [ 0.03,  0.26]
(sexFactor [male] × ageCenteredSquared) × sleep |      -0.13 | [-0.25, -0.02]
Code
performance::r2(quadraticGCM)
# R2 for Mixed Models

  Conditional R2: 0.947
     Marginal R2: 0.180

5.3 Plots

Code
newData <- expand.grid(
  sexFactor = c("female", "male"),
  age = 14:17,
  sleepFactor = c("low sleep", "average sleep", "high sleep")
)

newData$ageCentered <- newData$age - min(newData$age)
newData$ageCenteredSquared <- newData$ageCentered ^ 2

newData$sleep <- NA
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

newData$sleep[which(newData$sexFactor == "female" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

newData$sleep[which(newData$sexFactor == "female" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

newData$modelImpliedDepression <- predict( # predict.merMod
  quadraticGCM,
  newdata = newData,
  re.form = NA # include no random effects
)

5.3.1 Prototypical Growth Curve By Sex

Code
ggplot(
  data = newData |> filter(sleepFactor == "average sleep"),
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    color = sexFactor)) +
  geom_smooth( # fit quadratic smooth to the model-implied points; use geom_line() to fit a non-smooth line
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex

5.3.2 Individuals’ Growth Curves

Code
mydata_long$modelImpliedDepression <- predict( # predict.merMod
  quadraticGCM,
  newdata = mydata_long,
  re.form = NULL # include all random effects
)

ggplot(
  data = mydata_long,
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    group = ID,
    color = sexFactor)) +
  geom_smooth( # fit quadratic smooth to each participant's points
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 0.5
  ) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Individuals’ Model-Implied Growth Curves of Depression Symptoms

Individuals’ Model-Implied Growth Curves of Depression Symptoms

5.3.3 Individuals’ Trajectories Overlaid with Prototypical Trajectory By Sex

Code
ggplot(
  data = mydata_long,
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    group = ID)) +
  geom_smooth( # fit quadratic smooth to each participant's points
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 0.5,
    color = "gray",
    alpha = 0.30
  ) +
  geom_smooth( # fit quadratic smooth to the model-implied points; use geom_line() to fit a non-smooth line
    data = newData |> filter(sleepFactor == "average sleep"),
    mapping = aes(
      x = age,
      y = modelImpliedDepression,
      group = sexFactor,
      color = sexFactor),
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Individuals’ Model-Implied Growth Curves of Depression Symptoms Overlaid With Prototypical Growth Curves By Sex

Individuals’ Model-Implied Growth Curves of Depression Symptoms Overlaid With Prototypical Growth Curves By Sex

5.3.4 Prototypical Trajectory By Sex and Sleep

Code
ggplot(
  data = newData |> filter(sleepFactor != "average sleep"),
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    color = sexFactor,
    linetype = sleepFactor)) +
  geom_smooth( # fit quadratic smooth to the model-implied points; use geom_line() to fit a non-smooth line
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex",
    linetype = "Sleep"
  ) +
  guides(
    linetype = guide_legend(
      override.aes = list(color = "black")
    )
  ) +
  theme_classic()

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex and Sleep

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex and Sleep

6 Generalized Additive Mixed Model

Generalized additive mixed models (GAMMs) can be estimated using the mgcv::gamm() and gamm4::gamm4() functions.

6.1 Fit the Model

Code
num_cores <- parallelly::availableCores() - 1
num_true_cores <- parallelly::availableCores(logical = FALSE) - 1

In estimating generalized additive growth models of depression symptoms, this model includes:

  1. random intercepts and random linear and quadratic slopes for subjects,
  2. sex‐specific smooth functions of age, sex and sleep as fixed effects with their interaction, and
  3. sex‐specific smooth interactions between age and sleep, along with subject‐specific smooth deviations over age.
Code
gamGCM <- mgcv::gamm(
  depression ~ sexFactor + sleep + sexFactor:sleep +
    
    # sex-specific trajectories
    s(ageCentered, by = sexFactor, k = 4) +
    
    # sex-specific effects of sleep on trajectories
    ti(ageCentered, sleep, by = sexFactor, k = 4),
    
    # random effects
    random = list(ID = ~ 1 + ageCentered),
  data = mydata_long,
  nthreads = num_cores)

gamGCM_standardized <- mgcv::gamm(
  scale(depression) ~ sexFactor + scale(sleep) + sexFactor:scale(sleep) +
    
    # sex-specific trajectories
    s(scale(ageCentered), by = sexFactor, k = 4) +
    
    # sex-specific effects of sleep on trajectories
    ti(scale(ageCentered), scale(sleep), by = sexFactor, k = 4),
    
    # random effects
    random = list(ID = ~ 1 + scale(ageCentered)),
  data = mydata_long,
  nthreads = num_cores)

6.2 Model Summary

Code
gamGCM

Call:
NULL

--- 
Linear mixed-effects model fit by maximum likelihood
  Data: strip.offset(mf) 
  Log-likelihood: -2318.183
  Fixed: y ~ X - 1 
                             X(Intercept) 
                               19.7858582 
                           XsexFactormale 
                               -5.8872940 
                                   Xsleep 
                               -0.9818995 
                     XsexFactormale:sleep 
                                0.1574553 
       Xs(ageCentered):sexFactorfemaleFx1 
                                1.0448345 
         Xs(ageCentered):sexFactormaleFx1 
                               -0.7819394 
Xti(ageCentered,sleep):sexFactorfemaleFx1 
                                0.7604740 
  Xti(ageCentered,sleep):sexFactormaleFx1 
                                0.1158859 

Random effects:
 Formula: ~Xr - 1 | g
 Structure: pdIdnot
             Xr1      Xr2
StdDev: 3.912446 3.912446

 Formula: ~Xr.0 - 1 | g.0 %in% g
 Structure: pdIdnot
           Xr.01    Xr.02
StdDev: 3.985657 3.985657

 Formula: ~Xr.1 - 1 | g.1 %in% g.0 %in% g
 Structure: pdTens
           Xr.11    Xr.12    Xr.13    Xr.14    Xr.15    Xr.16    Xr.17    Xr.18
StdDev: 44618.83 21980.19 16659.22 41523.98 41392.15 14714.95 14338.73 3306.106

 Formula: ~Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g
 Structure: pdTens
           Xr.21    Xr.22    Xr.23    Xr.24    Xr.25    Xr.26    Xr.27    Xr.28
StdDev: 55490.18 55490.18 55490.18 11012.31 4.804365 11012.31 1.664289 11012.31

 Formula: ~1 + ageCentered | ID %in% g.2 %in% g.1 %in% g.0 %in% g
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr  
(Intercept) 3.827861 (Intr)
ageCentered 1.577057 0.531 
Residual    1.493607       

Number of Observations: 931
Number of Groups: 
                                   g                           g.0 %in% g 
                                   1                                    1 
                 g.1 %in% g.0 %in% g         g.2 %in% g.1 %in% g.0 %in% g 
                                   1                                    1 
ID %in% g.2 %in% g.1 %in% g.0 %in% g 
                                 250 
--- 

Family: gaussian 
Link function: identity 

Formula:
depression ~ sexFactor + sleep + sexFactor:sleep + s(ageCentered, 
    by = sexFactor, k = 4) + ti(ageCentered, sleep, by = sexFactor, 
    k = 4)

Estimated degrees of freedom:
2.81 2.77 1.00 2.27  total = 12.85
Code
gamGCM_standardized

Call:
NULL

--- 
Linear mixed-effects model fit by maximum likelihood
  Data: strip.offset(mf) 
  Log-likelihood: -559.9014
  Fixed: y ~ X - 1 
                                           X(Intercept) 
                                             0.30652577 
                                         XsexFactormale 
                                            -0.69948925 
                                          Xscale(sleep) 
                                            -0.16808561 
                            XsexFactormale:scale(sleep) 
                                             0.02695381 
              Xs(scale(ageCentered)):sexFactorfemaleFx1 
                                             0.15806708 
                Xs(scale(ageCentered)):sexFactormaleFx1 
                                            -0.11829512 
Xti(scale(ageCentered),scale(sleep)):sexFactorfemaleFx1 
                                             0.11504791 
  Xti(scale(ageCentered),scale(sleep)):sexFactormaleFx1 
                                             0.01753180 

Random effects:
 Formula: ~Xr - 1 | g
 Structure: pdIdnot
              Xr1       Xr2
StdDev: 0.5918887 0.5918887

 Formula: ~Xr.0 - 1 | g.0 %in% g
 Structure: pdIdnot
            Xr.01     Xr.02
StdDev: 0.6029715 0.6029715

 Formula: ~Xr.1 - 1 | g.1 %in% g.0 %in% g
 Structure: pdTens
           Xr.11    Xr.12    Xr.13    Xr.14    Xr.15    Xr.16    Xr.17    Xr.18
StdDev: 10014.16 4167.731 2462.413 9718.988 9706.694 3397.839 3362.515 488.6785

 Formula: ~Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g
 Structure: pdTens
           Xr.21    Xr.22    Xr.23    Xr.24     Xr.25    Xr.26     Xr.27
StdDev: 15279.14 15279.14 15279.14 3032.222 0.7268107 3032.222 0.2517759
           Xr.28
StdDev: 3032.222

 Formula: ~1 + scale(ageCentered) | ID %in% g.2 %in% g.1 %in% g.0 %in% g
 Structure: General positive-definite, Log-Cholesky parametrization
                   StdDev    Corr  
(Intercept)        0.8144423 (Intr)
scale(ageCentered) 0.2663372 0.798 
Residual           0.2259592       

Number of Observations: 931
Number of Groups: 
                                   g                           g.0 %in% g 
                                   1                                    1 
                 g.1 %in% g.0 %in% g         g.2 %in% g.1 %in% g.0 %in% g 
                                   1                                    1 
ID %in% g.2 %in% g.1 %in% g.0 %in% g 
                                 250 
--- 

Family: gaussian 
Link function: identity 

Formula:
scale(depression) ~ sexFactor + scale(sleep) + sexFactor:scale(sleep) + 
    s(scale(ageCentered), by = sexFactor, k = 4) + ti(scale(ageCentered), 
    scale(sleep), by = sexFactor, k = 4)

Estimated degrees of freedom:
2.81 2.77 1.00 2.27  total = 12.85
Code
performance::r2(gamGCM)
  R2: 0.211

6.3 Plots

Code
newData <- expand.grid(
  sexFactor = c("female", "male"),
  age = seq(
    from = min(mydata_long$age, na.rm = TRUE),
    to = max(mydata_long$age, na.rm = TRUE),
    length.out = 10000
  ),
  sleepFactor = c("low sleep", "average sleep", "high sleep")
)

newData$ageCentered <- newData$age - min(newData$age)

newData$sleep <- NA
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 17 & newData$sleepFactor == "average sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

newData$sleep[which(newData$sexFactor == "female" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 17 & newData$sleepFactor == "low sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE) - sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

newData$sleep[which(newData$sexFactor == "female" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 14 & newData$age < 15 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 14)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 15 & newData$age < 16 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 15)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 16 & newData$age < 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 16)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "female" & newData$age >= 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "female" & mydata_long$age == 17)], na.rm = TRUE)
newData$sleep[which(newData$sexFactor == "male" & newData$age >= 17 & newData$sleepFactor == "high sleep")] <- mean(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE) + sd(mydata_long$sleep[which(mydata_long$sexFactor == "male" & mydata_long$age == 17)], na.rm = TRUE)

# add a dummy id
newData$ID <- mydata_long$ID[1]

newData$modelImpliedDepression <- predict( # predict.gam
  gamGCM,
  newdata = newData
)

6.3.1 Prototypical Growth Curve By Sex

Code
ggplot(
  data = newData |> filter(sleepFactor == "average sleep"),
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    color = sexFactor)) +
  geom_smooth( # fit GAM smooth to the model-implied points; use geom_line() to fit a non-smooth line
    method = "gam",
    se = FALSE,
    linewidth = 0.5
  ) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex

6.3.2 Individuals’ Growth Curves

Code
mydata_long$modelImpliedDepression <- predict( # predict.gam
  gamGCM,
  newdata = mydata_long
)

gamMethod <- function(formula, data, weights, ...) {
  if(nrow(data) >= 4){ # gam smooth if 4+ timepoints
    mgcv::gam(y ~ s(x, bs = "cs", k = 3), data = data, ...)
  } else if(nrow(data) == 3){  # quadratic smooth if 3 timepoints
    lm(y ~ x + I(x^2), data = data, ...)
  } else{ # linear if fewer than 3 timepoints
    lm(y ~ x, data = data, ...)
  }
}

ggplot(
  data = mydata_long,
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    group = ID,
    color = sexFactor)) +
  geom_smooth( # fit GAM smooth to each participant's points
    method = gamMethod,
    se = FALSE,
    linewidth = 0.5
  ) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Individuals’ Model-Implied Growth Curves of Depression Symptoms

Individuals’ Model-Implied Growth Curves of Depression Symptoms

6.3.3 Individuals’ Trajectories Overlaid with Prototypical Trajectory By Sex

Code
ggplot(
  data = mydata_long,
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    group = ID)) +
  geom_smooth( # fit GAM smooth to each participant's points
    method = gamMethod,
    se = FALSE,
    linewidth = 0.5,
    color = "gray",
    alpha = 0.30
  ) +
  geom_smooth( # fit GAM smooth to the model-implied points; use geom_line() to fit a non-smooth line
    data = newData |> filter(sleepFactor == "average sleep"),
    mapping = aes(
      x = age,
      y = modelImpliedDepression,
      group = sexFactor,
      color = sexFactor),
    method = "gam",
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex"
  ) +
  theme_classic()

Individuals’ Model-Implied Growth Curves of Depression Symptoms Overlaid With Prototypical Growth Curves By Sex

Individuals’ Model-Implied Growth Curves of Depression Symptoms Overlaid With Prototypical Growth Curves By Sex

6.3.4 Prototypical Trajectory By Sex and Sleep

Code
ggplot(
  data = newData |> filter(sleepFactor != "average sleep"),
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    color = sexFactor,
    linetype = sleepFactor)) +
  geom_smooth( # fit GAM smooth to the model-implied points; use geom_line() to fit a non-smooth line
    method = "gam",
    se = FALSE) +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex",
    linetype = "Sleep"
  ) +
  guides(
    linetype = guide_legend(
      override.aes = list(color = "black")
    )
  ) +
  theme_classic()

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex and Sleep

Prototypical Model-Implied Growth Curves of Depression Symptoms By Sex and Sleep

7 Compare Models

Code
anova(
  linearGCM,
  quadraticGCM
)
Code
AIC(
  linearGCM,
  quadraticGCM,
  gamGCM
)
Code
MuMIn::AICc(
  linearGCM,
  quadraticGCM,
  gamGCM
)

8 Disaggregation of Within-Person and Between-Person Effects

8.1 Pre-Model Computations

Code
mydata_long <- mydata_long |>
  group_by(ID) |>
  mutate(
    sleep_personMean = mean(sleep, na.rm = TRUE),
    sleep_personMeanCentered = sleep - sleep_personMean
  ) |>
  ungroup()

8.2 Fit the Model

In disaggregating within- and between-person effects, this model fits:

  1. random intercepts and random linear slopes for subjects,
  2. sex and sleep_personMean as time-invariant fixed-effect predictors of the intercepts and slopes, and
  3. sleep_personMeanCentered as a time-varying fixed-effect predictor of depression.
Code
withinBetweenPersonSleepModel_basic <- lmerTest::lmer(
  depression ~ sleep_personMean + sexFactor +              # time-invariant predictors of the intercepts
    sleep_personMeanCentered + ageCentered +               # time-varying predictors
    sleep_personMean:ageCentered + sexFactor:ageCentered + # time-invariant predictors of the slopes
    (1 + ageCentered | ID),
  data = mydata_long,
  na.action = na.exclude)

However, based on our modeling in Section 5 and our model comparison in Section 7, we determined that a quadratic model fit significantly better than a linear model, so we can extend the model to a quadratic model.

In disaggregating within- and between-person effects, this model fits:

  1. random intercepts and random linear slopes and fixed quadratic slopes for subjects,
  2. sex and sleep_personMean (and their interaction) as time-invariant fixed-effect predictors of the intercepts and slopes, and
  3. sleep_personMeanCentered as a time-varying fixed-effect predictor of depression, with effects allowed to vary over age and by sex, including two‐ and three‐way interactions with age and sex.

The * symbol (unlike the : symbol) estimates all lower-order main effects and interactions, to simplify the syntax.

Code
withinBetweenPersonSleepModel <- lmerTest::lmer(
  depression ~ sexFactor*sleep_personMean*ageCentered + 
    sexFactor*sleep_personMean*ageCenteredSquared + 
    sexFactor*sleep_personMeanCentered*ageCentered + 
    sexFactor*sleep_personMeanCentered*ageCenteredSquared + 
    (1 + ageCentered | ID), # model fails to converge with a random effect for ageCenteredSquared
  data = mydata_long,
  na.action = na.exclude)

8.3 Model Summary

Code
summary(withinBetweenPersonSleepModel)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: depression ~ sexFactor * sleep_personMean * ageCentered + sexFactor *  
    sleep_personMean * ageCenteredSquared + sexFactor * sleep_personMeanCentered *  
    ageCentered + sexFactor * sleep_personMeanCentered * ageCenteredSquared +  
    (1 + ageCentered | ID)
   Data: mydata_long

REML criterion at convergence: 4626.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5691 -0.4181 -0.0233  0.3799  3.1684 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 ID       (Intercept) 14.331   3.786         
          ageCentered  2.404   1.551    0.53 
 Residual              2.237   1.496         
Number of obs: 931, groups:  ID, 250

Fixed effects:
                                                           Estimate Std. Error
(Intercept)                                                24.00570    2.85852
sexFactormale                                              -5.01340    4.32597
sleep_personMean                                           -1.60234    0.35245
ageCentered                                                 6.31714    2.04483
ageCenteredSquared                                         -0.70171    0.56569
sleep_personMeanCentered                                   -1.23490    0.31669
sexFactormale:sleep_personMean                              0.33208    0.53553
sexFactormale:ageCentered                                  -9.59764    3.11058
sleep_personMean:ageCentered                               -0.89401    0.25235
sexFactormale:ageCenteredSquared                            2.42488    0.86957
sleep_personMean:ageCenteredSquared                         0.16374    0.06980
sexFactormale:sleep_personMeanCentered                      0.03113    0.45428
ageCentered:sleep_personMeanCentered                        1.09341    0.55380
ageCenteredSquared:sleep_personMeanCentered                -0.39972    0.18201
sexFactormale:sleep_personMean:ageCentered                  0.99055    0.38523
sexFactormale:sleep_personMean:ageCenteredSquared          -0.30129    0.10778
sexFactormale:ageCentered:sleep_personMeanCentered         -0.47191    0.77720
sexFactormale:ageCenteredSquared:sleep_personMeanCentered   0.25192    0.25677
                                                                 df t value
(Intercept)                                               261.31526   8.398
sexFactormale                                             261.53048  -1.159
sleep_personMean                                          261.37873  -4.546
ageCentered                                               667.32798   3.089
ageCenteredSquared                                        447.12110  -1.240
sleep_personMeanCentered                                  616.86197  -3.899
sexFactormale:sleep_personMean                            261.42255   0.620
sexFactormale:ageCentered                                 668.34418  -3.085
sleep_personMean:ageCentered                              667.55202  -3.543
sexFactormale:ageCenteredSquared                          453.87104   2.789
sleep_personMean:ageCenteredSquared                       446.91031   2.346
sexFactormale:sleep_personMeanCentered                    619.99241   0.069
ageCentered:sleep_personMeanCentered                      625.42997   1.974
ageCenteredSquared:sleep_personMeanCentered               644.97219  -2.196
sexFactormale:sleep_personMean:ageCentered                668.20437   2.571
sexFactormale:sleep_personMean:ageCenteredSquared         453.82812  -2.796
sexFactormale:ageCentered:sleep_personMeanCentered        641.12293  -0.607
sexFactormale:ageCenteredSquared:sleep_personMeanCentered 657.31531   0.981
                                                                     Pr(>|t|)
(Intercept)                                               0.00000000000000292
sexFactormale                                                        0.247552
sleep_personMean                                          0.00000835813054210
ageCentered                                                          0.002090
ageCenteredSquared                                                   0.215463
sleep_personMeanCentered                                             0.000107
sexFactormale:sleep_personMean                                       0.535731
sexFactormale:ageCentered                                            0.002116
sleep_personMean:ageCentered                                         0.000424
sexFactormale:ageCenteredSquared                                     0.005516
sleep_personMean:ageCenteredSquared                                  0.019425
sexFactormale:sleep_personMeanCentered                               0.945382
ageCentered:sleep_personMeanCentered                                 0.048780
ageCenteredSquared:sleep_personMeanCentered                          0.028436
sexFactormale:sleep_personMean:ageCentered                           0.010345
sexFactormale:sleep_personMean:ageCenteredSquared                    0.005402
sexFactormale:ageCentered:sleep_personMeanCentered                   0.543940
sexFactormale:ageCenteredSquared:sleep_personMeanCentered            0.326910
                                                             
(Intercept)                                               ***
sexFactormale                                                
sleep_personMean                                          ***
ageCentered                                               ** 
ageCenteredSquared                                           
sleep_personMeanCentered                                  ***
sexFactormale:sleep_personMean                               
sexFactormale:ageCentered                                 ** 
sleep_personMean:ageCentered                              ***
sexFactormale:ageCenteredSquared                          ** 
sleep_personMean:ageCenteredSquared                       *  
sexFactormale:sleep_personMeanCentered                       
ageCentered:sleep_personMeanCentered                      *  
ageCenteredSquared:sleep_personMeanCentered               *  
sexFactormale:sleep_personMean:ageCentered                *  
sexFactormale:sleep_personMean:ageCenteredSquared         ** 
sexFactormale:ageCentered:sleep_personMeanCentered           
sexFactormale:ageCenteredSquared:sleep_personMeanCentered    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
print(effectsize::standardize_parameters(
  withinBetweenPersonSleepModel,
  method = "refit"),
  digits = 2)
# Standardization method: refit

Parameter                                                          | Std. Coef. |         95% CI
------------------------------------------------------------------------------------------------
(Intercept)                                                        |       0.31 | [ 0.18,  0.45]
sexFactor [male]                                                   |      -0.71 | [-0.91, -0.51]
sleep personMean                                                   |      -0.34 | [-0.48, -0.21]
ageCentered                                                        |      -0.14 | [-0.23, -0.06]
ageCenteredSquared                                                 |       0.32 | [ 0.25,  0.39]
sleep personMeanCentered                                           |      -0.09 | [-0.11, -0.06]
sexFactor [male] × sleep personMean                                |       0.11 | [-0.09,  0.31]
sexFactor [male] × ageCentered                                     |      -0.28 | [-0.40, -0.16]
sleep personMean × ageCentered                                     |      -0.15 | [-0.23, -0.06]
sexFactor [male] × ageCenteredSquared                              |   3.62e-03 | [-0.10,  0.11]
sleep personMean × ageCenteredSquared                              |       0.08 | [ 0.01,  0.15]
sexFactor [male] × sleep personMeanCentered                        |       0.02 | [-0.02,  0.05]
ageCentered × sleep personMeanCentered                             |       0.11 | [ 0.00,  0.22]
ageCenteredSquared × sleep personMeanCentered                      |      -0.12 | [-0.23, -0.01]
(sexFactor [male] × sleep personMean) × ageCentered                |       0.16 | [ 0.04,  0.28]
(sexFactor [male] × sleep personMean) × ageCenteredSquared         |      -0.15 | [-0.26, -0.05]
(sexFactor [male] × ageCentered) × sleep personMeanCentered        |      -0.05 | [-0.20,  0.11]
(sexFactor [male] × ageCenteredSquared) × sleep personMeanCentered |       0.08 | [-0.08,  0.23]
Code
performance::r2(withinBetweenPersonSleepModel)
# R2 for Mixed Models

  Conditional R2: 0.950
     Marginal R2: 0.246

8.4 Plots

8.4.1 Between-Person Effect of Sleep

Code
newData <- expand.grid(
  sexFactor = c("female", "male"),
  age = seq(
    from = min(mydata_long$age, na.rm = TRUE),
    to = max(mydata_long$age, na.rm = TRUE),
    length.out = 10000
  ),
  sleepFactor = c("low sleep", "high sleep")
)

newData$ageCentered <- newData$age - min(newData$age)
newData$ageCenteredSquared <- newData$ageCentered ^ 2

newData$sleep_personMean <- NA
newData$sleep_personMean[which(newData$sleepFactor == "low sleep" & newData$sexFactor == "female")] <- mean(mydata_long$sleep_personMean[which(mydata_long$sexFactor == "female")], na.rm = TRUE) - sd(mydata_long$sleep_personMean[which(mydata_long$sexFactor == "female")], na.rm = TRUE)

newData$sleep_personMean[which(newData$sleepFactor == "low sleep" & newData$sexFactor == "male")] <- mean(mydata_long$sleep_personMean[which(mydata_long$sexFactor == "male")], na.rm = TRUE) - sd(mydata_long$sleep_personMean[which(mydata_long$sexFactor == "male")], na.rm = TRUE)

newData$sleep_personMean[which(newData$sleepFactor == "high sleep" & newData$sexFactor == "female")] <- mean(mydata_long$sleep_personMean[which(mydata_long$sexFactor == "female")], na.rm = TRUE) + sd(mydata_long$sleep_personMean[which(mydata_long$sexFactor == "female")], na.rm = TRUE)

newData$sleep_personMean[which(newData$sleepFactor == "high sleep" & newData$sexFactor == "male")] <- mean(mydata_long$sleep_personMean[which(mydata_long$sexFactor == "male")], na.rm = TRUE) + sd(mydata_long$sleep_personMean[which(mydata_long$sexFactor == "male")], na.rm = TRUE)

newData$sleep_personMeanCentered <- 0

newData$modelImpliedDepression <- predict( # predict.merMod
  withinBetweenPersonSleepModel,
  newdata = newData,
  re.form = NA # include no random effects
)

ggplot(
  data = newData,
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    color = sexFactor,
    linetype = sleepFactor)) +
  geom_line() +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex",
    linetype = "Sleep"
  ) +
  theme_classic()

Between-Person Association Between Sleep and Depression

Between-Person Association Between Sleep and Depression

8.4.2 Within-Person Effect of Sleep

Code
newData <- expand.grid(
  sexFactor = c("female", "male"),
  age = seq(
    from = min(mydata_long$age, na.rm = TRUE),
    to = max(mydata_long$age, na.rm = TRUE),
    length.out = 10000
  ),
  sleepFactor = c("low sleep", "high sleep")
)

newData$ageCentered <- newData$age - min(newData$age)
newData$ageCenteredSquared <- newData$ageCentered ^ 2

newData$sleep_personMean <- mean(mydata_long$sleep_personMean, na.rm = TRUE)

newData$sleep_personMeanCentered <- NA
newData$sleep_personMeanCentered[which(newData$sleepFactor == "low sleep" & newData$sexFactor == "female")] <- mean(mydata_long$sleep_personMeanCentered[which(mydata_long$sexFactor == "female")], na.rm = TRUE) - sd(mydata_long$sleep_personMeanCentered[which(mydata_long$sexFactor == "female")], na.rm = TRUE)

newData$sleep_personMeanCentered[which(newData$sleepFactor == "low sleep" & newData$sexFactor == "male")] <- mean(mydata_long$sleep_personMeanCentered[which(mydata_long$sexFactor == "male")], na.rm = TRUE) - sd(mydata_long$sleep_personMeanCentered[which(mydata_long$sexFactor == "male")], na.rm = TRUE)

newData$sleep_personMeanCentered[which(newData$sleepFactor == "high sleep" & newData$sexFactor == "female")] <- mean(mydata_long$sleep_personMeanCentered[which(mydata_long$sexFactor == "female")], na.rm = TRUE) + sd(mydata_long$sleep_personMeanCentered[which(mydata_long$sexFactor == "female")], na.rm = TRUE)

newData$sleep_personMeanCentered[which(newData$sleepFactor == "high sleep" & newData$sexFactor == "male")] <- mean(mydata_long$sleep_personMeanCentered[which(mydata_long$sexFactor == "male")], na.rm = TRUE) + sd(mydata_long$sleep_personMeanCentered[which(mydata_long$sexFactor == "male")], na.rm = TRUE)

newData$modelImpliedDepression <- predict( # predict.merMod
  withinBetweenPersonSleepModel,
  newdata = newData,
  re.form = NA # include no random effects
)

ggplot(
  data = newData,
  mapping = aes(
    x = age,
    y = modelImpliedDepression,
    color = sexFactor,
    linetype = sleepFactor)) +
  geom_line() +
  labs(
    x = "Age (years)",
    y = "Depression Symptoms",
    color = "Sex",
    linetype = "Sleep"
  ) +
  theme_classic()

Within-Person Association Between Sleep and Depression

Within-Person Association Between Sleep and Depression



Developmental Psychopathology Lab