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

## Import Data

```{r}
#| eval: false

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

```{r}
#| include: false

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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