hlm.qmd

---
title: "Hierarchical Linear Modeling"
---

# Preamble

## Install Libraries

```{r}
#install.packages("remotes")
#remotes::install_github("DevPsyLab/petersenlab")
```

## Load Libraries

```{r}
library("petersenlab")
library("lme4")
library("nlme")
library("lmerTest")
library("MASS")
library("MCMCglmm")
library("mgcv")
library("parallelly")
library("performance")
library("ggplot2")
library("ggforce")
library("tidyverse")
```

## Import Data

```{r}
#| eval: false

mydata <- read.csv("https://osf.io/cqn3d/download")
mydata_long <- read.csv("./data/mydata_long.csv")
```

```{r}
#| include: false

mydata <- read.csv("./data/nlsy_math_long.csv") #https://osf.io/cqn3d/download
mydata_long <- read.csv("./data/mydata_long.csv")
```

## Simulate Data

```{r}
set.seed(52242)

mydata$outcome <- rpois(nrow(mydata), 4)
```

# Terms

These models go by a variety of different terms:

- hierarchical linear model (HLM)
- multilevel model (MLM)
- mixed effects model
- mixed model

# Overview

- <https://isaactpetersen.github.io/Principles-Psychological-Assessment/reliability.html#sec-mixedModels>
- <https://isaactpetersen.github.io/Fantasy-Football-Analytics-Textbook/mixed-models.html>

# 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.

```{r}
# Timepoint
mydata <- mydata |>
  arrange(id, age) |> 
  group_by(id) |>
  mutate(timepoint = 1:n())

# Age
mydata$ageYears <- mydata$age / 12
mydata$ageMonthsCentered <- mydata$age - min(mydata$age, na.rm = TRUE)

mydata$ageYearsCentered <- mydata$ageMonthsCentered / 12
mydata$ageYearsCenteredSquared <- mydata$ageYearsCentered ^ 2

# Age at T1
ageAtT1 <- mydata |> 
  ungroup() |> 
  filter(timepoint == 1) |> 
  mutate( # if we had date of birth, we could compute birth cohort (e.g., birth year), if we had time of measurement, we could compute period (e.g., year) at T1
    ageAtT1 = ageYears, # use age at baseline (T1) or birth year to reduce bias from attrition-related missing data
    ageAtT1Centered = ageAtT1 - min(ageAtT1, na.rm = TRUE),
    ageAtT1CenteredSquared = ageAtT1Centered ^ 2) |>
  select(id, ageAtT1, ageAtT1Centered, ageAtT1CenteredSquared)

mydata <- left_join(
  mydata,
  ageAtT1,
  by = c("id")
)

# Time in Study
mydata$timeInStudy <- mydata$ageYears - mydata$ageAtT1
mydata$timeInStudySquared <- mydata$timeInStudy ^ 2

# Sex
mydata$sex <- factor(
  mydata$female,
  levels = c(0, 1),
  labels = c("male", "female")
)
```

# Estimator: ML or REML

For small sample sizes, restricted maximum likelihood (REML) is preferred over maximum likelihood (ML).
ML preferred when there is a small number (< 4) of fixed effects; REML is preferred when there are more (> 4) fixed effects.
The greater the number of fixed effects, the greater the difference between REML and ML estimates.
Likelihood ratio (LR) tests for REML require exactly the same fixed effects specification in both models.
So, to compare models with different fixed effects with an LR test (to determine whether to include a particular fixed effect), ML must be used.
In contrast to the maximum likelihood estimation, REML can produce unbiased estimates of variance and covariance parameters, variance estimates are larger in REML than ML.
To compare whether an effect should be fixed or random, use REML.
To simultaneously compare fixed and random effects, use ML.

# Plot Observed Growth Curves

## Spaghetti Plot

```{r}
ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = math,
    group = id)) +
  geom_line() +
  labs(
    x = "Age (years)",
    y = "Math Score"
  ) +
  theme_classic()
```

## Faceted Plot

```{r}
facetedPlots <- mydata |>
  ggplot(
    aes(
      x = ageYears,
      y = math,
      group = id
    )
  ) +
  geom_point() +
  geom_smooth(
    method = "lm"
  ) +
  coord_cartesian(
    ylim = c(min(mydata$math, na.rm = TRUE), max(mydata$math, na.rm = TRUE))) +
  labs(
    x = "Age (years)",
    y = "Math Score"
  ) +
  theme_classic() +
  ggforce::facet_wrap_paginate(
    ~ id,
    ncol = 4,
    nrow = 3)

num_pages <- ggforce::n_pages(facetedPlots)
```

```{r}
#| eval: false

for(i in 1:num_pages){
  print(facetedPlots + 
    ggforce::facet_wrap_paginate(
      ~ id,
      ncol = 4,
      nrow = 3,
      page = i))
}
```

```{r}
#| echo: false

for(i in 1:2){
  print(facetedPlots + 
    ggforce::facet_wrap_paginate(
      ~ id,
      ncol = 4,
      nrow = 3,
      page = i))
}
```

# Linear Mixed Models {#sec-linear}

The following models are models that are fit in a linear mixed modeling framework.

## Growth Curve Models {#sec-gcm}

### Linear Growth Curve Model {#sec-linearGCM}

#### `lme4`

```{r}
linearMixedModel <- lmer(
  math ~ sex + ageYearsCentered + sex:ageYearsCentered + (1 + ageYearsCentered | id), # random intercepts and slopes; sex as a fixed-effect predictor of the intercepts and slopes
  data = mydata,
  REML = FALSE, #for ML
  na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa"))

summary(linearMixedModel)

print(effectsize::standardize_parameters(
  linearMixedModel,
  method = "refit"),
  digits = 2)

performance::model_performance(linearMixedModel)
performance::r2(linearMixedModel)
```

##### Prototypical Growth Curve

```{r}
newData <- expand.grid(
  female = c(0, 1),
  ageYears = c(
    min(mydata$ageYears, na.rm = TRUE),
    max(mydata$ageYears, na.rm = TRUE))
)

newData$ageYearsCentered <- newData$ageYears - min(newData$ageYears)

newData$sex <- NA
newData$sex[which(newData$female == 0)] <- "male"
newData$sex[which(newData$female == 1)] <- "female"
newData$sex <- as.factor(newData$sex)

newData$predictedValue <- predict( # predict.merMod
  linearMixedModel,
  newdata = newData,
  re.form = NA
)

ggplot(
  data = newData,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    color = sex)) +
  geom_line() +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

##### Individuals' Growth Curves

```{r}
mydata$predictedValue <- predict(
  linearMixedModel,
  newdata = mydata,
  re.form = NULL
)

ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id,
    color = sex)) +
  geom_line(
    alpha = 0.4
  ) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

##### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id)) +
  geom_line( # individuals' trajectories
    color = "gray",
    alpha = 0.5
  ) +
  geom_line( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = ageYears,
      y = predictedValue,
      group = sex,
      color = sex),
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

##### Extract Random Effects {#sec-extractRandomEffects}

```{r}
#| output: false

ranef(linearMixedModel)
```

```{r}
#| echo: false

linearMixedModel_re <- ranef(linearMixedModel)

linearMixedModel_re$id <- head(linearMixedModel_re$id)
linearMixedModel_re
```

##### Extract Participant-Specific Coefficients {#sec-extractParticipantCoefficients}

```{r}
#| output: false

coef(linearMixedModel)
```

```{r}
#| echo: false

linearMixedModel_coef <- coef(linearMixedModel)

linearMixedModel_coef$id <- head(linearMixedModel_coef$id)
linearMixedModel_coef
```

#### `nlme`

```{r}
linearMixedModel_nlme <- lme(
  math ~ sex + ageYearsCentered + sex:ageYearsCentered, # sex as a fixed-effect predictor of the intercepts and slopes
  random = ~ 1 + ageYearsCentered|id, # random intercepts and slopes
  data = mydata,
  method = "ML",
  na.action = na.exclude)

summary(linearMixedModel_nlme)

print(effectsize::standardize_parameters(
  linearMixedModel_nlme,
  method = "refit"),
  digits = 2)

performance::model_performance(linearMixedModel_nlme)
performance::r2(linearMixedModel_nlme)
```

#### R-Squared {#sec-rSquared}

The marginal $R^2$ is the proportion of variance in the outcome variable that is accounted for by the fixed effects in the model, whereas the conditional $R^2$ is the proportion of variance in the outcome variable that is accounted for by the fixed and random effects in the model.

```{r}
performance::r2(linearMixedModel)
performance::r2(linearMixedModel_nlme)
```

#### Intraclass Correlation Coefficent {#sec-icc}

```{r}
icc(linearMixedModel)
icc(linearMixedModel_nlme)
```

### Growth Curve Model with Timepoint-Specific Errors {#sec-timepointSpecificErrorsGCM}

Adapted from Usami & Murayama (2018):

```{r}
timepointSpecificErrorsMixedModel <- lmer(
  math ~ sex + ageYearsCentered + sex:ageYearsCentered + (1 | id) + (1 | ageYearsCentered), # timepoint-specific errors: observations are cross-classified with person and timepoint; sex as a fixed-effect predictor of the intercepts and slopes
  data = mydata,
  REML = FALSE, #for ML
  na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa"))

summary(timepointSpecificErrorsMixedModel)

print(effectsize::standardize_parameters(
  timepointSpecificErrorsMixedModel,
  method = "refit"),
  digits = 2)

performance::model_performance(timepointSpecificErrorsMixedModel)
performance::r2(timepointSpecificErrorsMixedModel)
```

### Quadratic Growth Curve Model {#sec-quadraticGCM}

When using higher-order polynomials, we could specify contrast codes for time to reduce multicollinearity between the linear and quadratic growth factors: <https://tdjorgensen.github.io/SEM-in-Ed-compendium/ch27.html#saturated-growth-model> (archived at <https://perma.cc/88YL-USLC>)

```{r}
factorLoadings <- poly(
  x = c(0,1,2,3), # times (can allow unequal spacing)
  degree = 2)

factorLoadings

linearLoadings <- factorLoadings[,1]
quadraticLoadings <- factorLoadings[,2]

linearLoadings
quadraticLoadings
```

#### Fit Model

```{r}
quadraticGCM <- lmer(
  math ~ sex + ageYearsCentered + ageYearsCenteredSquared + sex:ageYearsCentered + sex:ageYearsCenteredSquared + (1 + ageYearsCentered | id), # random intercepts and linear slopes; fixed quadratic slopes; sex as a fixed-effect predictor of the intercepts and slopes
  data = mydata,
  REML = FALSE, #for ML
  na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa"))

summary(quadraticGCM)

print(effectsize::standardize_parameters(
  quadraticGCM,
  method = "refit"),
  digits = 2)

performance::model_performance(quadraticGCM)
performance::r2(quadraticGCM)
```

This is equivalent to:

```{r}
quadraticGCM <- lmer(
  math ~ sex + ageYearsCentered + I(ageYearsCentered^2) + sex:ageYearsCentered + sex:I(ageYearsCentered^2) + (1 + ageYearsCentered | id), # random intercepts and slopes; sex as a fixed-effect predictor of the intercepts and slopes
  data = mydata,
  REML = FALSE, #for ML
  na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa"))

summary(quadraticGCM)

print(effectsize::standardize_parameters(
  quadraticGCM,
  method = "refit"),
  digits = 2)

performance::model_performance(quadraticGCM)
performance::r2(quadraticGCM)
```

#### Prototypical Growth Curve

```{r}
newData <- expand.grid(
  female = c(0, 1),
  ageYears = seq(from = min(mydata$ageYears, na.rm = TRUE), to = max(mydata$ageYears, na.rm = TRUE), length.out = 10000))

newData$ageYearsCentered <- newData$ageYears - min(newData$ageYears)
newData$ageYearsCenteredSquared <- newData$ageYearsCentered ^ 2

newData$sex <- NA
newData$sex[which(newData$female == 0)] <- "male"
newData$sex[which(newData$female == 1)] <- "female"
newData$sex <- as.factor(newData$sex)

newData$predictedValue <- predict( # predict.merMod
  quadraticGCM,
  newdata = newData,
  re.form = NA
)

ggplot(
  data = newData,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    color = sex)) +
  geom_line() +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Individuals' Growth Curves

```{r}
mydata$predictedValue <- predict(
  quadraticGCM,
  newdata = mydata,
  re.form = NULL
)

ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id,
    color = sex)) +
  geom_line(
    stat = "smooth",
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 0.5,
    alpha = 0.4
  ) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id)) +
  geom_line( # individuals' trajectories
    stat = "smooth",
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 0.5,
    color = "gray",
    alpha = 0.4
  ) +
  geom_line( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = ageYears,
      y = predictedValue,
      group = sex,
      color = sex),
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Extract Random Effects

```{r}
#| output: false

ranef(quadraticGCM)
```

```{r}
#| echo: false

quadraticGCM_re <- ranef(quadraticGCM)

quadraticGCM_re$id <- head(quadraticGCM_re$id)
quadraticGCM_re
```

### Logarithmic Growth Curve Model {#sec-logGCM}

“an exponential pattern of change—in which change appears to ‘level off’ over time—can be approximated through linear (and potentially quadratic) slopes for a natural-log-transformed” version of time (Hoffman, 2025).

#### Fit Model

```{r}
logGCM <- lmer(
  math ~ sex + log(ageYearsCentered + 1) + sex:log(ageYearsCentered + 1) + (1 + log(ageYearsCentered + 1) | id), # random intercepts and logarithmic slopes; sex as a fixed-effect predictor of the intercepts and slopes
  data = mydata,
  REML = FALSE, #for ML
  na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa"))

summary(logGCM)

print(effectsize::standardize_parameters(
  logGCM,
  method = "refit"),
  digits = 2)

performance::model_performance(logGCM)
performance::r2(logGCM)
```

#### Prototypical Growth Curve

```{r}
newData <- expand.grid(
  female = c(0, 1),
  ageYears = seq(from = min(mydata$ageYears, na.rm = TRUE), to = max(mydata$ageYears, na.rm = TRUE), length.out = 10000))

newData$ageYearsCentered <- newData$ageYears - min(newData$ageYears)

newData$sex <- NA
newData$sex[which(newData$female == 0)] <- "male"
newData$sex[which(newData$female == 1)] <- "female"
newData$sex <- as.factor(newData$sex)

newData$predictedValue <- predict( # predict.merMod
  logGCM,
  newdata = newData,
  re.form = NA
)

ggplot(
  data = newData,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    color = sex)) +
  geom_line() +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Individuals' Growth Curves

```{r}
mydata$predictedValue <- predict(
  logGCM,
  newdata = mydata,
  re.form = NULL
)

ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id,
    color = sex)) +
  geom_line(
    stat = "smooth",
    method = "lm",
    formula = y ~ log(x + 1),
    se = FALSE,
    linewidth = 0.5,
    alpha = 0.4
  ) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id)) +
  geom_line( # individuals' trajectories
    stat = "smooth",
    method = "lm",
    formula = y ~ log(x + 1),
    se = FALSE,
    linewidth = 0.5,
    color = "gray",
    alpha = 0.4
  ) +
  geom_line( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = ageYears,
      y = predictedValue,
      group = sex,
      color = sex),
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Extract Random Effects

```{r}
#| output: false

ranef(logGCM)
```

```{r}
#| echo: false

logGCM_re <- ranef(logGCM)

logGCM_re$id <- head(logGCM_re$id)
logGCM_re
```

### Spline Growth Curve Model {#sec-splineGCM}

#### Create Knot

```{r}
mydata$knot <- NA
mydata$knot[which(mydata$ageYears <= 11)] <- 0
mydata$knot[which(mydata$ageYears > 11)] <- 1
```

#### Fit Model

```{r}
splineGCM <- lmer(
  math ~ sex + ageYearsCentered + sex:ageYearsCentered + knot + knot:ageYearsCentered + sex:knot + sex:knot:ageYearsCentered + (1 + ageYearsCentered | id), # random intercepts and linear slopes; fixed quadratic slopes; sex as a fixed-effect predictor of the intercepts and slopes
  data = mydata,
  REML = FALSE, #for ML
  na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa"))

summary(splineGCM)

print(effectsize::standardize_parameters(
  splineGCM,
  method = "refit"),
  digits = 2)

performance::model_performance(splineGCM)
performance::r2(splineGCM)
```

#### Prototypical Growth Curve

```{r}
newData <- expand.grid(
  female = c(0, 1),
  ageYears = seq(from = min(mydata$ageYears, na.rm = TRUE), to = max(mydata$ageYears, na.rm = TRUE), length.out = 10000))

newData$ageYearsCentered <- newData$ageYears - min(newData$ageYears)

newData$knot <- NA
newData$knot[which(newData$ageYears <= 11)] <- 0
newData$knot[which(newData$ageYears > 11)] <- 1

newData$sex <- NA
newData$sex[which(newData$female == 0)] <- "male"
newData$sex[which(newData$female == 1)] <- "female"
newData$sex <- as.factor(newData$sex)

newData$predictedValue <- predict( # predict.merMod
  splineGCM,
  newdata = newData,
  re.form = NA
)

ggplot(
  data = newData,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    color = sex)) +
  geom_line() +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Individuals' Growth Curves

```{r}
mydata$predictedValue <- predict(
  splineGCM,
  newdata = mydata,
  re.form = NULL
)

ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id,
    color = sex)) +
  geom_line(
    alpha = 0.4
  ) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id)) +
  geom_line( # individuals' trajectories
    color = "gray",
    alpha = 0.5
  ) +
  geom_line( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = ageYears,
      y = predictedValue,
      group = sex,
      color = sex),
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Extract Random Effects

```{r}
#| output: false

ranef(splineGCM)
```

```{r}
#| echo: false

splineGCM_re <- ranef(splineGCM)

splineGCM_re$id <- head(splineGCM_re$id)
splineGCM_re
```

### Accelerated Longitudinal Study: Accounting for Age at T1 and Time in Study {#sec-acceleratedLongitudinal}

Per Hoffman ([here](https://www.cilvr.umd.edu/Conference2010/CILVRposts2010/Hoffman_post.pdf), archived at <https://perma.cc/6X8R-EX44>; and [here](https://www.lesahoffman.com/PSQF7375_AdvLong/PSQF7375_AdvLong_Lecture3_AltTime.pdf), archived at <https://perma.cc/25JY-JE8P>), use age at baseline (i.e., T1; or birth year) instead of mean age to lessen bias from attrition‐related missing data.

#### Linear

```{r}
acceleratedLongitudinalMixedModel_linear <- lmer(
  math ~ timeInStudy*sex*ageAtT1Centered + (1 + timeInStudy | id),
  data = mydata,
  REML = FALSE, #for ML
  na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa"))

summary(acceleratedLongitudinalMixedModel_linear)

print(effectsize::standardize_parameters(
  acceleratedLongitudinalMixedModel_linear,
  method = "refit"),
  digits = 2)

performance::model_performance(acceleratedLongitudinalMixedModel_linear)
performance::r2(acceleratedLongitudinalMixedModel_linear)
```

##### Prototypical Growth Curve

```{r}
newData <- unique(mydata[,c("ageYears","ageAtT1","timeInStudy")])

newData$ageAtT1Centered <- newData$ageAtT1 - min(mydata$ageAtT1, na.rm = TRUE)

newData <- tidyr::expand_grid(
  newData,
  female = c(0, 1)
)

newData$sex <- NA
newData$sex[which(newData$female == 0)] <- "male"
newData$sex[which(newData$female == 1)] <- "female"
newData$sex <- as.factor(newData$sex)

newData$predictedValue <- predict( # predict.merMod
  acceleratedLongitudinalMixedModel_linear,
  newdata = newData,
  re.form = NA
)

ggplot(
  data = newData,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    color = sex)) +
  geom_smooth() + #for jagged line: geom_line()
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

##### Individuals' Growth Curves

```{r}
mydata$predictedValue <- predict(
  acceleratedLongitudinalMixedModel_linear,
  newdata = mydata,
  re.form = NULL
)

ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id,
    color = sex)) +
  geom_line(
    alpha = 0.4
  ) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

##### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id)) +
  geom_line( # individuals' trajectories
    color = "gray",
    alpha = 0.5
  ) +
  geom_smooth( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = ageYears,
      y = predictedValue,
      group = sex,
      color = sex),
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

#### Quadratic

```{r}
acceleratedLongitudinalMixedModel_quadratic <- lmer(
  math ~ timeInStudy*sex*ageAtT1Centered + timeInStudySquared*sex*ageAtT1Centered + ageAtT1CenteredSquared + (1 + timeInStudy | id),
  data = mydata,
  REML = FALSE, #for ML
  na.action = na.exclude,
  control = lmerControl(optimizer = "bobyqa"))

summary(acceleratedLongitudinalMixedModel_quadratic)

print(effectsize::standardize_parameters(
  acceleratedLongitudinalMixedModel_quadratic,
  method = "refit"),
  digits = 2)

performance::model_performance(acceleratedLongitudinalMixedModel_quadratic)
performance::r2(acceleratedLongitudinalMixedModel_quadratic)
```

##### Prototypical Growth Curve

```{r}
newData <- unique(mydata[,c("ageYears","ageAtT1","timeInStudy")])

newData$ageAtT1Centered <- newData$ageAtT1 - min(mydata$ageAtT1, na.rm = TRUE)
newData$timeInStudySquared <- newData$timeInStudy ^ 2
newData$ageAtT1CenteredSquared <- newData$ageAtT1Centered ^2

newData <- tidyr::expand_grid(
  newData,
  female = c(0, 1)
)

newData$sex <- NA
newData$sex[which(newData$female == 0)] <- "male"
newData$sex[which(newData$female == 1)] <- "female"
newData$sex <- as.factor(newData$sex)

newData$predictedValue <- predict( # predict.merMod
  acceleratedLongitudinalMixedModel_quadratic,
  newdata = newData,
  re.form = NA
)

ggplot(
  data = newData,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    color = sex)) +
  geom_smooth() + #for jagged line: geom_line()
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

##### Individuals' Growth Curves

```{r}
mydata$predictedValue <- predict(
  acceleratedLongitudinalMixedModel_quadratic,
  newdata = mydata,
  re.form = NULL
)

ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id,
    color = sex)) +
  geom_smooth(
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 0.5,
    alpha = 0.4
  ) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

##### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id)) +
  geom_smooth( # individuals' trajectories
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 0.4,
    color = "gray",
    alpha = 0.5
  ) +
  geom_smooth( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = ageYears,
      y = predictedValue,
      group = sex,
      color = sex),
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

# Generalized Linear Mixed Models {#sec-generalized}

<https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html> (archived at <https://perma.cc/9RFS-BCE7>; source code: <https://github.com/bbolker/mixedmodels-misc/blob/master/glmmFAQ.rmd>)

## `lmer`

```{r}
generalizedLinearMixedModel <- glmer(
  outcome ~ female + ageYearsCentered + (ageYearsCentered | id),
  family = poisson(link = "log"),
  data = mydata,
  na.action = na.exclude)

summary(generalizedLinearMixedModel)

print(effectsize::standardize_parameters(
  generalizedLinearMixedModel,
  method = "refit"),
  digits = 2)

performance::model_performance(generalizedLinearMixedModel)
performance::r2(generalizedLinearMixedModel)
```

## `MASS`

```{r}
glmmPQLmodel <- glmmPQL(
  outcome ~ female + ageYearsCentered,
  random = ~ 1 + ageYearsCentered|id,
  family = poisson(link = "log"),
  data = mydata)

summary(glmmPQLmodel)

print(effectsize::standardize_parameters(
  glmmPQLmodel,
  method = "refit"),
  digits = 2)
```

## `MCMCglmm`

```{r}
MCMCglmmModel <- MCMCglmm(
  outcome ~ female + ageYearsCentered,
  random = ~ us(ageYearsCentered):id,
  family = "poisson",
  data = na.omit(mydata[,c("id","outcome","female","ageYearsCentered")]))

summary(MCMCglmmModel)
```

# Nonlinear Mixed Models {#sec-nonlinear}

## Generalized Additive Mixed Model {#sec-gamm}

### Fit Model

```{r}
#| label: parallel-cores-mixed-models
#| include: false

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

```{r}
#| eval: false

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

```{r}
gamModel <- mgcv::gam(
  math ~ s(ageYearsCentered, by = sex) + s(id, ageYearsCentered, bs = "re"),
  data = mydata,
  nthreads = num_cores
)

summary(gamModel)

print(effectsize::standardize_parameters(
  gamModel,
  method = "refit"),
  digits = 2)

performance::model_performance(gamModel)
performance::r2(gamModel)
```

### Prototypical Growth Curve

```{r}
newData <- expand.grid(
  female = c(0, 1),
  ageYears = seq(from = min(mydata$ageYears, na.rm = TRUE), to = max(mydata$ageYears, na.rm = TRUE), length.out = 10000))

newData$ageYearsCentered <- newData$ageYears - min(newData$ageYears)
newData$ageYearsCenteredSquared <- newData$ageYearsCentered ^ 2

newData$sex <- NA
newData$sex[which(newData$female == 0)] <- "male"
newData$sex[which(newData$female == 1)] <- "female"
newData$sex <- as.factor(newData$sex)

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

newData$predictedValue <- predict(
  gamModel,
  newdata = newData,
  exclude = "s(ageYearsCentered,id)"
)

ggplot(
  data = newData,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    color = sex)) +
  geom_line() +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

### Individuals' Growth Curves

```{r}
mydata$predictedValue <- predict(
  gamModel,
  newdata = mydata
)

gamMethod <- function(formula, data, weights, ...) { # from here: https://stackoverflow.com/a/76762792/2029527 (archived at https://perma.cc/85C6-FKJ4)
  if(nrow(data) < 4) return(lm(y ~ x, data = data, ...))
  else  mgcv::gam(y ~ s(x, bs = "cs", k = 3), data = data, ...)
}

ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id,
    color = sex)) +
  geom_line(
    stat = "smooth",
    method = gamMethod,
    se = FALSE,
    linewidth = 0.5,
    alpha = 0.4
  ) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot(
  data = mydata,
  mapping = aes(
    x = ageYears,
    y = predictedValue,
    group = id)) +
  geom_line( # individuals' trajectories
    stat = "smooth",
    method = gamMethod,
    se = FALSE,
    linewidth = 0.5,
    color = "gray",
    alpha = 0.4
  ) +
  geom_line( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = ageYears,
      y = predictedValue,
      group = sex,
      color = sex),
    linewidth = 2) +
  labs(
    x = "Age (years)",
    y = "Math Score",
    color = "Sex"
  ) +
  theme_classic()
```

## Asymptotic Growth {#sec-nonlinearAsymptotic}

```{r}
nonlinearModel <- nlme(
  height ~ SSasymp(age, Asym, R0, lrc),
  data = Loblolly,
  fixed = Asym + R0 + lrc ~ 1,
  random = Asym ~ 1)

summary(nonlinearModel)
```

# Disaggregation of Within-Person and Between-Person Effects {#sec-disaggregatingWithinBetween}

## Pre-Model Computations

```{r}
#| code-fold: true

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
```

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

## Fit the Model

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

1. random intercepts and random linear slopes for subjects,
1. `sex` and `sleep_personMean` as time-invariant fixed-effect predictors of the intercepts and slopes, and
1. `sleep_personMeanCentered` as a time-varying fixed-effect predictor of depression.

```{r}
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)
```

We can also extend the model to a quadratic model (or to any other functional form of change).

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

1. random intercepts and random linear slopes and fixed quadratic slopes for subjects,
1. `sex` and `sleep_personMean` (and their interaction) as time-invariant fixed-effect predictors of the intercepts and slopes, and
1. `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.

```{r}
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)
```

## Model Summary

```{r}
summary(withinBetweenPersonSleepModel)

print(effectsize::standardize_parameters(
  withinBetweenPersonSleepModel,
  method = "refit"),
  digits = 2)

performance::r2(withinBetweenPersonSleepModel)
```

## Plots

### Between-Person Effect of Sleep

```{r}
#| fig-cap: "Between-Person Association Between Sleep and Depression"

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()
```

### Within-Person Effect of Sleep

```{r}
#| fig-cap: "Within-Person Association Between Sleep and Depression"

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()
```

# Moderated Multilevel Models {#sec-moderation}

[Simonsohn (2024)](https://doi.org/10.1177/25152459231207787)

```{r}
moderatedMultilevelModel <- mgcv::gam(
  math ~
    s(timeInStudy) +
    s(ageAtT1Centered) +
    ti(timeInStudy, ageAtT1Centered) +
    s(id, bs = "re") +
    s(id, by = timeInStudy, bs = "re"),
  data = mydata
)

summary(moderatedMultilevelModel)

print(effectsize::standardize_parameters(
  moderatedMultilevelModel,
  method = "refit"),
  digits = 2)

performance::model_performance(moderatedMultilevelModel)
performance::r2(moderatedMultilevelModel)
```

# Cross-Lagged Models

Multilevel cross-lagged models can be estimated using [multilevel structural equation modeling](mplus.qmd#sec-msem-lcm-sr).

# Robust Mixed Models

To evaluate the extent to which a finding could driven by outliers, this could be done in a number of different ways, such as:

- identifying and excluding influential observations based on DFBETAS or Cook’s distance (Nieuwenhuis, Grotenhuis, & Pelzer, 2012)
- fitting mixed models using rank-based estimation (Bilgic & Susmann, 2013; Finch, 2017) or robust estimating equations (Koller, 2016)
- estimating robust standard errors using a sandwich estimator (Wang & Merkle, 2018)

# Compare Models

For nested models:

```{r}
#| eval: false

anova(
  model1,
  model2,
  model3
)
```

For non-nested models:

```{r}
#| eval: false

performance::r2(model1)
performance::r2(model2)
performance::r2(model3)

performance::model_performance(model1)
performance::model_performance(model2)
performance::model_performance(model3)

AIC(
  model1,
  model2,
  model3
)

BIC(
  model1,
  model2,
  model3
)

bbmle::AICtab(
  model1,
  model2,
  model3
)

bbmle::AICctab(
  model1,
  model2,
  model3
)

bbmle::BICtab(
  model1,
  model2,
  model3
)

AICcmodavg::aictab(
  model1,
  model2,
  model3
)

AICcmodavg::bictab(
  model1,
  model2,
  model3
)

MuMIn::AICc(
  model1,
  model2,
  model3
)
```

# Assumptions

The within-group errors:

1. are independent
2. are identically normally distributed
3. have mean zero and variance sigma-squared
4. are independent of the random effects

The random effects:

5. are normally distributed
6. have mean zero and covariance matrix Psi (not depending on the group)
7. are independent for different groups

# Examining Model Assumptions

## Resources

Pinheiro and Bates (2000) book (p. 174, section 4.3.1)

<https://stats.stackexchange.com/questions/77891/checking-assumptions-lmer-lme-mixed-models-in-r> (archived at <https://perma.cc/J5GC-PCUT>)

## QQ Plots

Make QQ plots for each level of the random effects.
Vary the level from 0, 1, to 2 so that you can check the between- and within-subject residuals.

```{r}
qqnorm(
  linearMixedModel_nlme,
  ~ ranef(., level = 1))
```

## PP Plots

```{r}
ppPlot(linearMixedModel)
```

## QQ Plot of residuals

```{r}
qqnorm(resid(linearMixedModel))
qqline(resid(linearMixedModel))
```

## Plot residuals

```{r}
plot(linearMixedModel)
```

## Plot residuals by group (in the example below, level 2 represents the individual)

```{r}
plot(
  linearMixedModel,
  as.factor(id) ~ resid(.),
  abline = 0,
  xlab = "Residuals")
```

## Plot residuals by levels of a predictor

```{r}
plot(
  linearMixedModel_nlme,
  resid(., type = "p") ~ fitted(.) | female) #type = "p" specifies standardized residuals
```

## Can model heteroscedasticity of the within-group error with the weights argument

```{r}
linearMixedModel_nlmeVarStructure <- lme(
  math ~ female + ageYearsCentered,
  random = ~ 1 + ageYearsCentered|id,
  weights = varIdent(form = ~ 1 | female),
  method = "ML",
  data = mydata,
  na.action = na.exclude)

summary(linearMixedModel_nlmeVarStructure)
```

## Plot observed and fitted values

```{r}
plot(
  linearMixedModel,
  math ~ fitted(.))
```

## Plot QQ plot of residuals by levels of a predictor

```{r}
qqnorm(linearMixedModel_nlme, ~ resid(.) | female)
qqnorm(linearMixedModel_nlme, ~ resid(.))
```

## QQ plot of random effects

Make QQ plots for each level of the random effects.
Vary the level from 0, 1, to 2 so that you can check the between- and within-subject residuals.


```{r}
qqnorm(
  linearMixedModel_nlme,
  ~ ranef(., level = 0))

qqnorm(
  linearMixedModel_nlme,
  ~ ranef(., level = 1))

qqnorm(
  linearMixedModel_nlme,
  ~ ranef(., level = 2))
```

## QQ plot of random effects by levels of a predictor

```{r}
qqnorm(
  linearMixedModel_nlme,
  ~ ranef(., level = 1) | female)
```

## Pairs plot

```{r}
pairs(linearMixedModel_nlme)

pairs(
  linearMixedModel_nlme,
  ~ ranef(., level = 1) | female)
```

## Variance functions for modeling heteroscedasticity

- `varFixed`: fixed variance
- `varIdent`: different variances per stratum
- `varPower`: power of covariate
- `varExp`: exponential of covariate
- `varConstPower`: constant plus power of covariate
- `varComb`: combination of variance functions

## Correlation structures for modeling dependence

- `corCompSymm`: compound symmetry
- `corSymm`: general
- `corAR1`: autoregressive of order 1
- `corCAR1`: continuous-time AR(1)
- `corARMA`: autoregressive-moving average
- `corExp`: exponential
- `corGaus`: Gaussian
- `corLin`: linear
- `corRatio`: rational quadratic
- `corSpher`: spherical

# Power Analysis {#sec-powerAnalysis}

- `mlmpower`: <https://github.com/bkeller2/mlmpower> ([Enders et al., 2025](https://doi.org/10.1037/met0000614))
- <https://aguinis.shinyapps.io/ml_power/>
- <https://www.causalevaluation.org/power-analysis.html>
  - <https://powerupr.shinyapps.io/index/>
- <https://koumurayama.shinyapps.io/tmethod_mlm/>
- <https://webpower.psychstat.org/wiki/models/index>

# Session Info

```{r}
#| code-fold: true

sessionInfo()
```