sem.qmd

---
title: "Structural Equation Modeling"
---

# Preamble

For a resource on using `lavaan` in `R` for structural equation modeling, see the following e-book: <https://tdjorgensen.github.io/SEM-in-Ed-compendium>

## Install Libraries

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

## Load Libraries

```{r}
library("lavaan")
library("semTools")
library("semPlot")
library("lavaanPlot")
library("lavaangui")
library("lcsm")
library("MBESS")
library("tidyverse")
library("ggforce")
```

## Options

```{r}
options(scipen = 999)
```

# Simulate Data

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

sampleSize <- 100

X <- rnorm(sampleSize)
M <- 0.5*X + rnorm(sampleSize)
Y <- 0.7*M + rnorm(sampleSize)

mydata <- data.frame(
  X = X,
  Y = Y,
  M = M)
```

# Import data

```{r}
longitudinalMI <- read.csv("./data/Bliese-Ployhart-2002-indicators-1.csv")
```

# Overview

<https://isaactpetersen.github.io/Principles-Psychological-Assessment/structural-equation-modeling.html>

# Analysis examples

<https://isaactpetersen.github.io/Principles-Psychological-Assessment/structural-equation-modeling.html#sec-semModelExample-sem>

# Plot Observed Growth Curves

## Spaghetti Plot

Transform data from wide to long format:

```{r}
Demo.growth$id <- 1:nrow(Demo.growth)

Demo.growth_long <- Demo.growth %>% 
  pivot_longer(
    cols = c(t1,t2,t3,t4),
    names_to = "variable",
    values_to = "value",
    names_pattern = "t(.)") %>% 
  rename(
    timepoint = variable,
    score = value
  )

Demo.growth_long$timepoint <- as.numeric(Demo.growth_long$timepoint)
```

Plot the observed trajectory for each participant:

```{r}
ggplot(
  data = Demo.growth_long,
  mapping = aes(
    x = timepoint,
    y = score,
    group = id)) +
  geom_line() +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
  theme_classic()
```

## Faceted Plot

```{r}
facetedPlots <- Demo.growth_long |>
  ggplot(
    aes(
      x = timepoint,
      y = score,
      group = id
    )
  ) +
  geom_point() +
  geom_smooth(
    method = "lm"
  ) +
  coord_cartesian(
    ylim = c(min(Demo.growth_long$score, na.rm = TRUE), max(Demo.growth_long$score, na.rm = TRUE))) +
  labs(
    x = "Timepoint",
    y = "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))
}
```

# Latent Growth Curve Model {#sec-lgcm}

See the chapter on latent growth curve modeling in `lavaan`: <https://tdjorgensen.github.io/SEM-in-Ed-compendium/ch27.html> (archived at <https://perma.cc/C79K-BCLD>).

For extensions, see the following resources:

- growth curves using latent indicators (rather than observed indicators): <https://tdjorgensen.github.io/SEM-in-Ed-compendium/ch27.html#latent-indicators-of-growth-factors> (archived at <https://perma.cc/7HMC-AGNP>)
- growth curves using discrete indicators (rather than continuous indicators): <https://tdjorgensen.github.io/SEM-in-Ed-compendium/ch27.html#discrete-indicators> (archived at <https://perma.cc/LT84-82EY>)
- parallel growth curves: <https://tdjorgensen.github.io/SEM-in-Ed-compendium/ch27.html#parallel-growth-curves> (archived at <https://perma.cc/A6WA-627J>)

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

### Model Syntax

#### Abbreviated

```{r}
lgcm1_syntax <- '
  # Intercept and slope
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  slope =~ 0*t1 + 1*t2 + 2*t3 + 3*t4

  # Regression paths
  intercept ~ x1 + x2
  slope ~ x1 + x2
  
  # Time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
'
```

#### Full

```{r}
lgcm2_syntax <- '
  # Intercept and slope
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  slope =~ 0*t1 + 1*t2 + 2*t3 + 3*t4

  # Regression paths
  intercept ~ x1 + x2
  slope ~ x1 + x2
  
  # Time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
  
  # Constrain observed intercepts to zero
  t1 ~ 0*1
  t2 ~ 0*1
  t3 ~ 0*1
  t4 ~ 0*1
  
  # Estimate mean of intercept and slope
  intercept ~ 1
  slope ~ 1
'
```

### Fit the Model

#### Abbreviated

```{r}
lgcm1_fit <- growth(
  lgcm1_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  int.ov.free = FALSE,
  int.lv.free = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

#### Full

```{r}
lgcm2_fit <- sem(
  lgcm2_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

### Summary Output

#### Abbreviated

```{r}
summary(
  lgcm1_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

#### Full

```{r}
summary(
  lgcm2_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Estimates of Model Fit

```{r}
fitMeasures(
  lgcm1_fit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "chisq.scaling.factor",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr",
    "rmsea.robust", "cfi.robust", "tli.robust"))
```

### Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(
  lgcm1_fit,
  type = "cor")
```

### Modification Indices

```{r}
modificationindices(
  lgcm1_fit,
  sort. = TRUE)
```

### Internal Consistency Reliability

```{r}
compRelSEM(lgcm1_fit)
```

### Path Diagram

```{r}
semPlot::semPaths(
  lgcm1_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

lavaanPlot::lavaanPlot(
  lgcm1_fit,
  coefs = TRUE,
  #covs = TRUE,
  stand = TRUE)

lavaanPlot::lavaanPlot2(
  lgcm1_fit,
  #stand = TRUE, # currently throws error; uncomment out when fixed: https://github.com/alishinski/lavaanPlot/issues/52
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(lgcm1_fit)
```

### Plot Trajectories

#### Prototypical Growth Curve

```{r}
lgcm1_intercept <- coef(lgcm1_fit)["intercept~1"]
lgcm1_slope <- coef(lgcm1_fit)["slope~1"]

timepoints <- 4

newData <- expand.grid(
  time = c(1, 4)
)

newData$predictedValue <- NA
newData$predictedValue[which(newData$time == 1)] <- lgcm1_intercept
newData$predictedValue[which(newData$time == 4)] <- lgcm1_intercept + (timepoints - 1)*lgcm1_slope

ggplot(
  data = newData,
  mapping = aes(
    x = time,
    y = predictedValue)) +
  geom_line(
    linewidth = 1
  ) +
  coord_cartesian(
    ylim = c(0, 5)) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
  theme_classic()
```

#### Individuals' Growth Curves

```{r}
newData1 <- as.data.frame(predict(lgcm1_fit))
newData1$id <- row.names(newData1)

newData1$t1 <- newData1$intercept
newData1$t4 <- newData1$intercept + (timepoints - 1)*newData1$slope

newData2 <- pivot_longer(
  newData1,
  cols = c(t1, t4)) %>% 
  select(-intercept, -slope)

newData2$time <- NA
newData2$time[which(newData2$name == "t1")] <- 1
newData2$time[which(newData2$name == "t4")] <- 4

ggplot(
  data = newData2,
  mapping = aes(
    x = time,
    y = value,
    group = factor(id))) +
  geom_line() +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
   theme_classic()
```

#### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot() +
  geom_line( # individuals' trajectories
    data = newData2,
    mapping = aes(
      x = time,
      y = value,
    group = factor(id)),
    color = "gray",
    alpha = 0.5
  ) +
  geom_line( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = time,
      y = predictedValue
    ),
    linewidth = 2
  ) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
   theme_classic()
```

## Latent Basis Growth Curve Model {#sec-latentBasisLGCM}

### Model Syntax

#### Abbreviated

```{r}
lbgcm1_syntax <- '
  # Intercept and slope
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  slope =~ 0*t1 + a*t2 + b*t3 + 1*t4 # freely estimate the loadings for t2 and t3

  # Regression paths
  intercept ~ x1 + x2
  slope ~ x1 + x2
  
  # Time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
'
```

#### Full

```{r}
lbgcm2_syntax <- '
  # Intercept and slope
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  slope =~ 0*t1 + a*t2 + b*t3 + 1*t4 # freely estimate the loadings for t2 and t3

  # Regression paths
  intercept ~ x1 + x2
  slope ~ x1 + x2
  
  # Time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
  
  # Constrain observed intercepts to zero
  t1 ~ 0*1
  t2 ~ 0*1
  t3 ~ 0*1
  t4 ~ 0*1
  
  # Estimate mean of intercept and slope
  intercept ~ 1
  slope ~ 1
'
```

### Fit the Model

#### Abbreviated

```{r}
lbgcm1_fit <- growth(
  lbgcm1_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  int.ov.free = FALSE,
  int.lv.free = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

#### Full

```{r}
lbgcm2_fit <- sem(
  lbgcm2_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

### Summary Output

#### Abbreviated

```{r}
summary(
  lbgcm1_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

#### Full

```{r}
summary(
  lbgcm2_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Estimates of Model Fit

```{r}
fitMeasures(
  lbgcm1_fit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "chisq.scaling.factor",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr",
    "rmsea.robust", "cfi.robust", "tli.robust"))
```

### Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(
  lbgcm1_fit,
  type = "cor")
```

### Modification Indices

```{r}
modificationindices(
  lbgcm1_fit,
  sort. = TRUE)
```

### Internal Consistency Reliability

```{r}
compRelSEM(lbgcm1_fit)
```

### Path Diagram

```{r}
semPaths(
  lbgcm1_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

lavaanPlot::lavaanPlot(
  lbgcm1_fit,
  coefs = TRUE,
  #covs = TRUE,
  stand = TRUE)

lavaanPlot::lavaanPlot2(
  lbgcm1_fit,
  stand = TRUE,
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(lbgcm1_fit)
```

### Plot Trajectories

#### Prototypical Growth Curve

```{r}
lbgcm1_intercept <- coef(lbgcm1_fit)["intercept~1"]
lbgcm1_slope <- coef(lbgcm1_fit)["slope~1"]
lbgcm1_slopeloadingt2 <- coef(lbgcm1_fit)["a"]
lbgcm1_slopeloadingt3 <- coef(lbgcm1_fit)["b"]

timepoints <- 4

newData <- data.frame(
  time = 1:4,
  slopeloading = c(0, lbgcm1_slopeloadingt2, lbgcm1_slopeloadingt3, 1)
)

newData$predictedValue <- NA
newData$predictedValue <- lbgcm1_intercept + lbgcm1_slope * newData$slopeloading

ggplot(
  data = newData,
  mapping = aes(
    x = time,
    y = predictedValue)) +
  geom_line(
    linewidth = 1
  ) +
  coord_cartesian(
    ylim = c(0, 5)) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
  theme_classic()
```

#### Individuals' Growth Curves

```{r}
person_factors <- as.data.frame(predict(lbgcm1_fit))
person_factors$id <- rownames(person_factors)

slope_loadings <- c(0, lbgcm1_slopeloadingt2, lbgcm1_slopeloadingt3, 1)

# Compute model-implied values for each person at each time point
individual_trajectories <- person_factors %>%
  rowwise() %>%
  mutate(
    t1 = intercept + slope * slope_loadings[1],
    t2 = intercept + slope * slope_loadings[2],
    t3 = intercept + slope * slope_loadings[3],
    t4 = intercept + slope * slope_loadings[4]
  ) %>%
  ungroup() %>%
  select(id, t1, t2, t3, t4) %>%
  pivot_longer(
    cols = t1:t4,
    names_to = "timepoint",
    values_to = "value") %>%
  mutate(
    time = as.integer(substr(timepoint, 2, 2)) # extract number from "t1", "t2", etc.
  )

ggplot(
  data = individual_trajectories,
  mapping = aes(
    x = time,
    y = value,
    group = factor(id))) +
  geom_line() +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
   theme_classic()
```

#### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot() +
  geom_line( # individuals' model-implied trajectories
    data = individual_trajectories,
    aes(
      x = time,
      y = value,
      group = id),
    color = "gray",
    alpha = 0.5
  ) +
  geom_line( # prototypical trajectory
    data = newData,
    aes(
      x = time,
      y = predictedValue),
    color = "blue",
    linewidth = 2
  ) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
   theme_classic()
```

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

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/2CJX-WLLZ>)

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

### Model Syntax

#### Abbreviated

```{r}
quadraticGCM1_syntax <- '
  # Intercept and slope
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  linear =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
  quadratic =~ 0*t1 + 1*t2 + 4*t3 + 9*t4

  # Regression paths
  intercept ~ x1 + x2
  linear ~ x1 + x2
  quadratic ~ x1 + x2
  
  # Time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
'
```

#### Full

```{r}
quadraticGCM2_syntax <- '
  # Intercept and slope
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  linear =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
  quadratic =~ 0*t1 + 1*t2 + 4*t3 + 9*t4

  # Regression paths
  intercept ~ x1 + x2
  linear ~ x1 + x2
  quadratic ~ x1 + x2
  
  # Time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
  
  # Constrain observed intercepts to zero
  t1 ~ 0*1
  t2 ~ 0*1
  t3 ~ 0*1
  t4 ~ 0*1
  
  # Estimate mean of intercept and slope
  intercept ~ 1
  linear ~ 1
  quadratic ~ 1
'
```

### Fit the Model

#### Abbreviated

```{r}
quadraticGCM1_fit <- growth(
  quadraticGCM1_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  int.ov.free = FALSE,
  int.lv.free = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

#### Full

```{r}
quadraticGCM2_fit <- sem(
  quadraticGCM2_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

### Summary Output

#### Abbreviated

```{r}
summary(
  quadraticGCM1_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

#### Full

```{r}
summary(
  quadraticGCM2_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Estimates of Model Fit

```{r}
fitMeasures(
  quadraticGCM1_fit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "chisq.scaling.factor",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr",
    "rmsea.robust", "cfi.robust", "tli.robust"))
```

### Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(
  quadraticGCM1_fit,
  type = "cor")
```

### Modification Indices

```{r}
modificationindices(
  quadraticGCM1_fit,
  sort. = TRUE)
```

### Internal Consistency Reliability

```{r}
compRelSEM(quadraticGCM1_fit)
```

### Path Diagram

```{r}
semPlot::semPaths(
  quadraticGCM1_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

lavaanPlot::lavaanPlot(
  quadraticGCM1_fit,
  coefs = TRUE,
  #covs = TRUE,
  stand = TRUE)

lavaanPlot::lavaanPlot2(
  quadraticGCM1_fit,
  #stand = TRUE, # currently throws error; uncomment out when fixed: https://github.com/alishinski/lavaanPlot/issues/52
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(quadraticGCM1_fit)
```

### Plot Trajectories

#### Prototypical Growth Curve

Calculated from intercept and slope parameters:

```{r}
quadraticGCM1_intercept <- coef(quadraticGCM1_fit)["intercept~1"]
quadraticGCM1_linear <- coef(quadraticGCM1_fit)["linear~1"]
quadraticGCM1_quadratic <- coef(quadraticGCM1_fit)["quadratic~1"]

timepoints <- 4

newData <- data.frame(
  time = 1:4,
  linearLoading = c(0, 1, 2, 3),
  quadraticLoading = c(0, 1, 4, 9)
)

newData$predictedValue <- NA
newData$predictedValue <- quadraticGCM1_intercept + (quadraticGCM1_linear * newData$linearLoading) + (quadraticGCM1_quadratic * newData$quadraticLoading)

ggplot(
  data = newData,
  mapping = aes(
    x = time,
    y = predictedValue)) +
  geom_smooth(
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 1
  ) +
  coord_cartesian(
    ylim = c(0, 5)) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
  theme_classic()
```

#### Individuals' Growth Curves

```{r}
person_factors <- as.data.frame(predict(quadraticGCM1_fit))
person_factors$id <- rownames(person_factors)

linear_loadings <- c(0, 1, 2, 3)
quadratic_loadings <- c(0, 1, 4, 9)

# Compute model-implied values for each person at each time point
individual_trajectories <- person_factors %>%
  rowwise() %>%
  mutate(
    t1 = intercept + (linear * linear_loadings[1]) + (quadratic * quadratic_loadings[1]),
    t2 = intercept + (linear * linear_loadings[2]) + (quadratic * quadratic_loadings[2]),
    t3 = intercept + (linear * linear_loadings[3]) + (quadratic * quadratic_loadings[3]),
    t4 = intercept + (linear * linear_loadings[4]) + (quadratic * quadratic_loadings[4])
  ) %>%
  ungroup() %>%
  select(id, t1, t2, t3, t4) %>%
  pivot_longer(
    cols = t1:t4,
    names_to = "timepoint",
    values_to = "value") %>%
  mutate(
    time = as.integer(substr(timepoint, 2, 2)) # extract number from "t1", "t2", etc.
  )

ggplot(
  data = individual_trajectories,
  mapping = aes(
    x = time,
    y = value,
    group = factor(id))) +
  geom_smooth(
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 0.5
  ) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
   theme_classic()
```

#### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot() +
  geom_smooth( # individuals' model-implied trajectories
    data = individual_trajectories,
    mapping = aes(
      x = time,
      y = value,
      group = id),
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 0.5,
    color = "gray",
    alpha = 0.30
  ) +
  geom_smooth( # prototypical trajectory
    data = newData,
    mapping = aes(
      x = time,
      y = predictedValue),
    inherit.aes = FALSE,
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    linewidth = 2
  ) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
   theme_classic()
```

## Spline (Piecewise) Growth Curve Model {#sec-splineLGCM}

### Model Syntax

#### Abbreviated

```{r}
splineGCM1_syntax <- '
  # Intercept and slope
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  slope =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
  knot =~ 0*t1 + 0*t2 + 1*t3 + 1*t4

  # Regression paths
  intercept ~ x1 + x2
  slope ~ x1 + x2
  knot ~ x1 + x2
  
  # Spline has no variance
  knot ~~ 0*knot

  # Spline does not covary with intercept and slope
  knot ~~ 0*intercept
  knot ~~ 0*slope
  
  # Time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
'
```

#### Full

```{r}
splineGCM2_syntax <- '
  # Intercept and slope
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  slope =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
  knot =~ 0*t1 + 0*t2 + 1*t3 + 1*t4

  # Regression paths
  intercept ~ x1 + x2
  slope ~ x1 + x2
  knot ~ x1 + x2
  
  # Spline has no variance
  knot ~~ 0*knot
  
  # Spline does not covary with intercept and slope
  knot ~~ 0*intercept
  knot ~~ 0*slope
  
  # Time-varying covariates
  t1 ~ c1
  t2 ~ c2
  t3 ~ c3
  t4 ~ c4
  
  # Constrain observed intercepts to zero
  t1 ~ 0*1
  t2 ~ 0*1
  t3 ~ 0*1
  t4 ~ 0*1
  
  # Estimate mean of intercept and slope
  intercept ~ 1
  slope ~ 1
  knot ~ 1
'
```

### Fit the Model

#### Abbreviated

```{r}
splineGCM1_fit <- growth(
  splineGCM1_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  int.ov.free = FALSE,
  int.lv.free = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

#### Full

```{r}
splineGCM2_fit <- sem(
  splineGCM2_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

### Summary Output

#### Abbreviated

```{r}
summary(
  splineGCM1_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

#### Full

```{r}
summary(
  splineGCM2_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Estimates of Model Fit

```{r}
fitMeasures(
  splineGCM1_fit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "chisq.scaling.factor",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr",
    "rmsea.robust", "cfi.robust", "tli.robust"))
```

### Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(
  splineGCM1_fit,
  type = "cor")
```

### Modification Indices

```{r}
modificationindices(
  splineGCM1_fit,
  sort. = TRUE)
```

### Internal Consistency Reliability

```{r}
compRelSEM(splineGCM1_fit)
```

### Path Diagram

```{r}
semPlot::semPaths(
  splineGCM1_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

lavaanPlot::lavaanPlot(
  splineGCM1_fit,
  coefs = TRUE,
  #covs = TRUE,
  stand = TRUE)

lavaanPlot::lavaanPlot2(
  splineGCM1_fit,
  #stand = TRUE, # currently throws error; uncomment out when fixed: https://github.com/alishinski/lavaanPlot/issues/52
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(splineGCM1_fit)
```

### Plot Trajectories

#### Prototypical Growth Curve

Calculated from intercept and slope parameters:

```{r}
splineGCM1_intercept <- coef(splineGCM1_fit)["intercept~1"]
splineGCM1_slope <- coef(splineGCM1_fit)["slope~1"]
splineGCM1_knot <- coef(splineGCM1_fit)["knot~1"]

timepoints <- 4

newData <- data.frame(
  time = 1:4,
  linearLoading = c(0, 1, 2, 3),
  knotLoading = c(0, 0, 1, 1)
)

newData$predictedValue <- NA
newData$predictedValue <- splineGCM1_intercept + (splineGCM1_slope * newData$linearLoading) + (splineGCM1_knot * newData$knotLoading)

ggplot(
  data = newData,
  mapping = aes(
    x = time,
    y = predictedValue)) +
  geom_line(
    linewidth = 1
  ) +
  coord_cartesian(
    ylim = c(0, 5)) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
  theme_classic()
```

#### Individuals' Growth Curves

```{r}
person_factors <- as.data.frame(predict(splineGCM1_fit))
person_factors$id <- rownames(person_factors)

slope_loadings <- c(0, 1, 2, 3)
knot_loadings <- c(0, 0, 1, 1)

# Compute model-implied values for each person at each time point
individual_trajectories <- person_factors %>%
  rowwise() %>%
  mutate(
    t1 = intercept + (slope * slope_loadings[1]) + (knot * knot_loadings[1]),
    t2 = intercept + (slope * slope_loadings[2]) + (knot * knot_loadings[2]),
    t3 = intercept + (slope * slope_loadings[3]) + (knot * knot_loadings[3]),
    t4 = intercept + (slope * slope_loadings[4]) + (knot * knot_loadings[4])
  ) %>%
  ungroup() %>%
  select(id, t1, t2, t3, t4) %>%
  pivot_longer(
    cols = t1:t4,
    names_to = "timepoint",
    values_to = "value") %>%
  mutate(
    time = as.integer(substr(timepoint, 2, 2)) # extract number from "t1", "t2", etc.
  )

ggplot(
  data = individual_trajectories,
  mapping = aes(
    x = time,
    y = value,
    group = factor(id))) +
  geom_line() +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
   theme_classic()
```

#### Individuals' Trajectories Overlaid with Prototypical Trajectory

```{r}
ggplot() +
  geom_line( # individuals' model-implied trajectories
    data = individual_trajectories,
    aes(
      x = time,
      y = value,
      group = id),
    linewidth = 0.5,
    color = "gray",
    alpha = 0.30
  ) +
  geom_line( # prototypical trajectory
    data = newData,
    aes(
      x = time,
      y = predictedValue),
    color = "blue",
    linewidth = 2
  ) +
  labs(
    x = "Timepoint",
    y = "Score"
  ) +
   theme_classic()
```

## Saturated Growth Curve Model {#sec-saturatedGCM}

<https://tdjorgensen.github.io/SEM-in-Ed-compendium/ch27.html#saturated-growth-model> (archived at <https://perma.cc/2CJX-WLLZ>)

# Latent Change Score Model {#sec-lcsm}

To generate the syntax for a latent change score model, we use the [lcsm](https://doi.org/10.32614/CRAN.package.lcsm) package.

## Model Syntax

```{r}
bivariateLCSM_syntax <- specify_bi_lcsm(
  timepoints = 3,
  var_x = "x",
  model_x = list(
    alpha_constant = TRUE, # alpha = intercept (constant change factor)
    beta = TRUE, # beta = proportional change factor (latent true score predicting its change score)
    phi = TRUE), # phi = autoregression of change scores
  var_y = "y",
  model_y = list(
    alpha_constant = TRUE, # alpha = intercept (constant change factor)
    beta = TRUE, # beta = proportional change factor (latent true score predicting its change score)
    phi = TRUE), # phi = autoregression of change scores
  coupling = list(
    delta_lag_xy = TRUE,
    delta_lag_yx = TRUE),
  change_letter_x = "g",
  change_letter_y = "j")

cat(bivariateLCSM_syntax)
```

## Fit the Model

```{r}
bivariateLCSM_fit <- fit_bi_lcsm(
  data = data_bi_lcsm,
  var_x = names(data_bi_lcsm)[2:4],
  var_y = names(data_bi_lcsm)[12:14],
  model_x = list(
    alpha_constant = TRUE, # alpha = intercept (constant change factor)
    beta = TRUE, # beta = proportional change factor (latent true score predicting its change score)
    phi = TRUE), # phi = autoregression of change scores
  model_y = list(
    alpha_constant = TRUE, # alpha = intercept (constant change factor)
    beta = TRUE, # beta = proportional change factor (latent true score predicting its change score)
    phi = TRUE), # phi = autoregression of change scores
  coupling = list(
    delta_lag_xy = TRUE,
    xi_lag_yx = TRUE),
  fixed.x = FALSE
  )
```

## Summary Output

```{r}
summary(
  bivariateLCSM_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

## Estimates of Model Fit

```{r}
fitMeasures(
  bivariateLCSM_fit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "chisq.scaling.factor",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr",
    "rmsea.robust", "cfi.robust", "tli.robust"))
```

## Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(
  bivariateLCSM_fit,
  type = "cor")
```

## Modification Indices

```{r}
modificationindices(
  bivariateLCSM_fit,
  sort. = TRUE)
```

## Path Diagram

```{r}
semPaths(
  bivariateLCSM_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

plot_lcsm(
  lavaan_object = bivariateLCSM_fit,
  lcsm = "bivariate",
  lavaan_syntax = bivariateLCSM_syntax,
  edge.label.cex = .9,
  lcsm_colours = TRUE)

#lavaanPlot::lavaanPlot( # throws error
#  bivariateLCSM_fit,
#  coefs = TRUE,
#  #covs = TRUE,
#  stand = TRUE)

lavaanPlot::lavaanPlot2(
  bivariateLCSM_fit,
  stand = TRUE,
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(bivariateLCSM_fit)
```

## Plot Trajectories

```{r}
plot_trajectories(
  data_bi_lcsm,
  id_var = "id",
  var_list = c(
    "x1", "x2", "x3", "x4", "x5", 
    "x6", "x7", "x8", "x9", "x10"),
  xlab = "Time",
  ylab = "X Score",
  connect_missing = FALSE)

  plot_trajectories(
  data_bi_lcsm,
  id_var = "id",
  var_list = c(
    "y1", "y2", "y3", "y4", "y5",
    "y6", "y7", "y8", "y9", "y10"),
  xlab = "Time",
  ylab = "Y Score",
  connect_missing = FALSE)
```

# Cross-Lagged Panel Model {#sec-clpm}

## Model Syntax

```{r}
clpm_syntax <- '
  # Autoregressive Paths
  t4 ~ t3
  t3 ~ t2
  t2 ~ t1
  
  c4 ~ c3
  c3 ~ c2
  c2 ~ c1
  
  # Concurrent Covariances
  t1 ~~ c1
  t2 ~~ c2
  t3 ~~ c3
  t4 ~~ c4
  
  # Cross-Lagged Paths
  t4 ~ c3
  t3 ~ c2
  t2 ~ c1
  
  c4 ~ t3
  c3 ~ t2
  c2 ~ t1
'
```

## Fit the Model

```{r}
clpm_fit <- sem(
  clpm_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  std.lv = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

## Summary Output

```{r}
summary(
  clpm_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

## Estimates of Model Fit

```{r}
fitMeasures(
  clpm_fit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "chisq.scaling.factor",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr",
    "rmsea.robust", "cfi.robust", "tli.robust"))
```

## Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(
  clpm_fit,
  type = "cor")
```

## Modification Indices

```{r}
modificationindices(
  clpm_fit,
  sort. = TRUE)
```

## Path Diagram

```{r}
semPaths(
  clpm_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

lavaanPlot::lavaanPlot(
  clpm_fit,
  coefs = TRUE,
  #covs = TRUE,
  stand = TRUE)

lavaanPlot::lavaanPlot2(
  clpm_fit,
  #stand = TRUE, # currently throws error; uncomment out when fixed: https://github.com/alishinski/lavaanPlot/issues/52
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(clpm_fit)
```

# Random Intercept Cross-Lagged Panel Model {#sec-riclpm}

## Model Syntax

### Abbreviated

Adapted from Mulder & Hamaker (2021): <https://doi.org/10.1080/10705511.2020.1784738>

<https://jeroendmulder.github.io/RI-CLPM/lavaan.html> (archived at <https://perma.cc/2K6A-WUJQ>)

```{r}
riclpm1_syntax <- '
  # Random Intercepts
  t =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  c =~ 1*c1 + 1*c2 + 1*c3 + 1*c4
  
  # Create Within-Person Centered Variables
  wt1 =~ 1*t1
  wt2 =~ 1*t2
  wt3 =~ 1*t3
  wt4 =~ 1*t4
  
  wc1 =~ 1*c1
  wc2 =~ 1*c2
  wc3 =~ 1*c3
  wc4 =~ 1*c4
  
  # Autoregressive Paths
  wt4 ~ wt3
  wt3 ~ wt2
  wt2 ~ wt1
  
  wc4 ~ wc3
  wc3 ~ wc2
  wc2 ~ wc1
  
  # Concurrent Covariances
  wt1 ~~ wc1
  wt2 ~~ wc2
  wt3 ~~ wc3
  wt4 ~~ wc4
  
  # Cross-Lagged Paths
  wt4 ~ wc3
  wt3 ~ wc2
  wt2 ~ wc1
  
  wc4 ~ wt3
  wc3 ~ wt2
  wc2 ~ wt1
  
  # Variance and Covariance of Random Intercepts
  t ~~ t
  c ~~ c
  t ~~ c
  
  # Variances of Within-Person Centered Variables
  wt1 ~~ wt1
  wt2 ~~ wt2
  wt3 ~~ wt3
  wt4 ~~ wt4
  
  wc1 ~~ wc1
  wc2 ~~ wc2
  wc3 ~~ wc3
  wc4 ~~ wc4
'
```

### Full

Adapted from Mund & Nestler (2017): <https://osf.io/a4dhk>

```{r}
riclpm2_syntax <- '
  # Random Intercepts
  t =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  c =~ 1*c1 + 1*c2 + 1*c3 + 1*c4
  
  # Create Within-Person Centered Variables
  wt1 =~ 1*t1
  wt2 =~ 1*t2
  wt3 =~ 1*t3
  wt4 =~ 1*t4
  
  wc1 =~ 1*c1
  wc2 =~ 1*c2
  wc3 =~ 1*c3
  wc4 =~ 1*c4
  
  # Autoregressive Paths
  wt4 ~ wt3
  wt3 ~ wt2
  wt2 ~ wt1
  
  wc4 ~ wc3
  wc3 ~ wc2
  wc2 ~ wc1
  
  # Concurrent Covariances
  wt1 ~~ wc1
  wt2 ~~ wc2
  wt3 ~~ wc3
  wt4 ~~ wc4
  
  # Cross-Lagged Paths
  wt4 ~ wc3
  wt3 ~ wc2
  wt2 ~ wc1
  
  wc4 ~ wt3
  wc3 ~ wt2
  wc2 ~ wt1
  
  # Variance and Covariance of Random Intercepts
  t ~~ t
  c ~~ c
  t ~~ c
  
  # Variances of Within-Person Centered Variables
  wt1 ~~ wt1
  wt2 ~~ wt2
  wt3 ~~ wt3
  wt4 ~~ wt4
  
  wc1 ~~ wc1
  wc2 ~~ wc2
  wc3 ~~ wc3
  wc4 ~~ wc4
  
  # Fix Error Variances of Observed Variables to Zero
  t1 ~~ 0*t1
  t2 ~~ 0*t2
  t3 ~~ 0*t3
  t4 ~~ 0*t4
  
  c1 ~~ 0*c1
  c2 ~~ 0*c2
  c3 ~~ 0*c3
  c4 ~~ 0*c4
  
  # Fix the Covariances Between the Random Intercepts and the Latents at T1 to Zero
  wt1 ~~ 0*t
  wt1 ~~ 0*c
  
  wc1 ~~ 0*t
  wc1 ~~ 0*c
  
  # Estimate Observed Intercepts
  t1 ~ 1
  t2 ~ 1
  t3 ~ 1
  t4 ~ 1
  
  c1 ~ 1
  c2 ~ 1
  c3 ~ 1
  c4 ~ 1
  
  # Fix the Means of the Latents to Zero
  wt1 ~ 0*1
  wt2 ~ 0*1
  wt3 ~ 0*1
  wt4 ~ 0*1
  
  wc1 ~ 0*1
  wc2 ~ 0*1
  wc3 ~ 0*1
  wc4 ~ 0*1
  
  t ~ 0*1
  c ~ 0*1
'
```

## Fit the Model

### Abbreviated

```{r}
riclpm1_fit <- lavaan(
  riclpm1_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  int.ov.free = TRUE,
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

### Full

```{r}
riclpm2_fit <- sem(
  riclpm2_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

## Summary Output

### Abbreviated

```{r}
summary(
  riclpm1_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Full

```{r}
summary(
  riclpm2_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

## Estimates of Model Fit

```{r}
fitMeasures(
  riclpm1_fit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "chisq.scaling.factor",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr",
    "rmsea.robust", "cfi.robust", "tli.robust"))
```

## Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(
  riclpm1_fit,
  type = "cor")
```

## Modification Indices

```{r}
modificationindices(
  riclpm1_fit,
  sort. = TRUE)
```

## Internal Consistency Reliability

```{r}
compRelSEM(riclpm1_fit)
```

## Path Diagram

```{r}
semPaths(
  riclpm1_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

lavaanPlot::lavaanPlot(
  riclpm1_fit,
  coefs = TRUE,
  #covs = TRUE,
  stand = TRUE)

lavaanPlot::lavaanPlot2(
  riclpm1_fit,
  #stand = TRUE, # currently throws error; uncomment out when fixed: https://github.com/alishinski/lavaanPlot/issues/52
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(riclpm1_fit)
```

# Latent Curve Model with Structured Residuals {#sec-lcm-sr}

A latent curve model with structured residuals (LCM-SR) is also called an autoregressive latent trajectory model with structured residuals (ALT-SR).
A LSM-SR model can also be estimated using [multilevel SEM](mplus.qmd#sec-msem-lcm-sr).

## Model Syntax

Adapted from Mund & Nestler (2017): <https://osf.io/a4dhk>

```{r}
lcmsr_syntax <- '
  # Define intercept and growth factors
  intercept.t =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
  slope.t =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
  
  intercept.c =~ 1*c1 + 1*c2 + 1*c3 + 1*c4
  slope.c =~ 0*c1 + 1*c2 + 2*c3 + 3*c4
  
  # Define phantom latent variables
  e.t1 =~ 1*t1
  e.t2 =~ 1*t2
  e.t3 =~ 1*t3
  e.t4 =~ 1*t4
  
  e.c1 =~ 1*c1
  e.c2 =~ 1*c2
  e.c3 =~ 1*c3
  e.c4 =~ 1*c4
  
  # Autoregressive paths
  e.t2 ~ a1*e.t1
  e.t3 ~ a1*e.t2
  e.t4 ~ a1*e.t3
  
  e.c2 ~ a2*e.c1
  e.c3 ~ a2*e.c2
  e.c4 ~ a2*e.c3
  
  # Cross-lagged paths
  e.c2 ~ c1*e.t1
  e.c3 ~ c1*e.t2
  e.c4 ~ c1*e.t3
  
  e.t2 ~ c2*e.c1
  e.t3 ~ c2*e.c2
  e.t4 ~ c2*e.c3
  
  # Some further constraints on the variance structure
  # 1. Set error variances of the observed variables to zero
  t1 ~~ 0*t1
  t2 ~~ 0*t2
  t3 ~~ 0*t3
  t4 ~~ 0*t4
  
  c1 ~~ 0*c1
  c2 ~~ 0*c2
  c3 ~~ 0*c3
  c4 ~~ 0*c4
  
  # 2. Let lavaan estimate the variance of the latent variables (residuals)
  e.t1 ~~ vart1*e.t1
  e.t2 ~~ vart2*e.t2
  e.t3 ~~ vart3*e.t3
  e.t4 ~~ vart4*e.t4
  
  e.c1 ~~ varc1*e.c1
  e.c2 ~~ varc2*e.c2
  e.c3 ~~ varc3*e.c3
  e.c4 ~~ varc4*e.c4
  
  # 3. We also want estimates of the intercept factor variances, the slope
  #    variances, and the covariances
  intercept.t ~~ varintercept.t*intercept.t
  intercept.c ~~ varintercept.c*intercept.c
  slope.t ~~ varslope.t*slope.t
  slope.c ~~ varslope.c*slope.c
  
  intercept.t ~~ covintercept*intercept.c
  slope.t ~~ covslope*slope.c
  
  intercept.t ~~ covintercept.tslope.t*slope.t
  intercept.t ~~ covintercept.tslope.c*slope.c
  intercept.c ~~ covintercept.cslope.t*slope.t
  intercept.c ~~ covintercept.cslope.c*slope.c
  
  # 4. We have to define that the covariance between the intercepts and
  #    the slopes and the latents of the first time point are zero
  e.t1 ~~ 0*intercept.t
  e.c1 ~~ 0*intercept.t
  e.t1 ~~ 0*slope.t
  e.c1 ~~ 0*slope.t
  e.t1 ~~ 0*intercept.c
  e.c1 ~~ 0*intercept.c
  e.t1 ~~ 0*slope.c
  e.c1 ~~ 0*slope.c
  
  # 5. Finally, we estimate the covariance between the latents of x and y
  #    of the first time point, the second time-point and so on. Note that
  #    for the second to fourth time point the correlation is constrained to
  #    the same value
  e.t1 ~~ cov1*e.c1
  e.t2 ~~ e1*e.c2
  e.t3 ~~ e1*e.c3
  e.t4 ~~ e1*e.c4
  
  # The model also contains a mean structure and we have to define some
  # constraints for this part of the model. The assumption is that we
  # only want estimates of the mean of the intercept factors. All other means
  # are defined to be zero:
  t1 ~ 0*1
  t2 ~ 0*1
  t3 ~ 0*1
  t4 ~ 0*1
  
  c1 ~ 0*1
  c2 ~ 0*1
  c3 ~ 0*1
  c4 ~ 0*1
  
  e.t1 ~ 0*1
  e.t2 ~ 0*1
  e.t3 ~ 0*1
  e.t4 ~ 0*1
  
  e.c1 ~ 0*1
  e.c2 ~ 0*1
  e.c3 ~ 0*1
  e.c4 ~ 0*1
  
  intercept.t ~ 1
  intercept.c ~ 1
  slope.t ~ 1
  slope.c ~ 1
'
```

## Fit the Model

```{r}
lcmsr_fit <- sem(
  lcmsr_syntax,
  data = Demo.growth,
  missing = "ML",
  estimator = "MLR",
  fixed.x = FALSE,
  em.h1.iter.max = 100000)
```

## Summary Output

```{r}
summary(
  lcmsr_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

## Estimates of Model Fit

```{r}
fitMeasures(
  lcmsr_fit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "chisq.scaling.factor",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr",
    "rmsea.robust", "cfi.robust", "tli.robust"))
```

## Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(
  lcmsr_fit,
  type = "cor")
```

## Modification Indices

```{r}
modificationindices(
  lcmsr_fit,
  sort. = TRUE)
```

## Internal Consistency Reliability

```{r}
compRelSEM(lcmsr_fit)
```

## Path Diagram

```{r}
semPaths(
  lcmsr_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

#lavaanPlot::lavaanPlot( # throws error
#  lcmsr_fit,
#  coefs = TRUE,
#  #covs = TRUE,
#  stand = TRUE)

lavaanPlot::lavaanPlot2(
  lcmsr_fit,
  #stand = TRUE, # currently throws error; uncomment out when fixed: https://github.com/alishinski/lavaanPlot/issues/52
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(lcmsr_fit)
```

# Mediation {#sec-mediation}

## Overview

It is important not to just consider mediation effects where there is a bivariate association between the predictor and the outcome.
Sometimes ***inconsistent mediation*** occurs where the indirect effect has an opposite sign (i.e., positive or negative) from the direct or total effect.
The idea is there there may be multiple mediating mechanisms, and that the predictor/cause may have both beneficial and harmful effects on the outcome, and that the mediating mechanisms can "cancel each other out" in terms of the total effect.
For instance, the predictor may help via mechanism A, but may hurt via mechanism B.
In this case, the total effect may be weak or small, even though there may be strong mediating effects.

## Model Syntax

```{r}
mediationModel <- '
  # direct effect (cPrime)
  Y ~ direct*X
  
  # mediator
  M ~ a*X
  Y ~ b*M
  
  # indirect effect = a*b
  indirect := a*b
  
  # total effect (c)
  total := direct + indirect
  totalAbs := abs(direct) + abs(indirect)
  
  # proportion mediated
  Pm := abs(indirect) / totalAbs
'
```

## Fit the Model

To get a robust estimate of the indirect effect, we obtain bootstrapped estimates from 1,000 bootstrap draws.
Typically, we would obtain bootstrapped estimates from 10,000 bootstrap draws, but this example uses only 1,000 bootstrap draws for a shorter runtime.

```{r}
mediationFit <- sem(
  mediationModel,
  data = mydata,
  se = "bootstrap",
  bootstrap = 1000, # generally use 10,000 bootstrap draws; this example uses 1,000 for speed
  parallel = "multicore", # parallelization for speed: use "multicore" for Mac/Linux; "snow" for PC
  iseed = 52242, # for reproducibility
  missing = "ML",
  estimator = "ML",
  # std.lv = TRUE, # for models with latent variables
  fixed.x = FALSE)
```

## Summary Output

```{r}
summary(
  mediationFit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

## Parameter Estimates

### Bias-Corrected Bootstrap

Adjusted bootstrap percentile (BCa) method, but with no correction for acceleration (only for bias):

```{r}
mediationFit_estimates_bca <- parameterEstimates(
  mediationFit,
  boot.ci.type = "bca.simple",
  standardized = TRUE)

mediationFit_estimates <- mediationFit_estimates_bca

mediationFit_estimates_bca
```

### Percentile Bootstrap

```{r}
mediationFit_estimates_perc <- parameterEstimates(
  mediationFit,
  boot.ci.type = "perc",
  standardized = TRUE)

mediationFit_estimates_perc
```

## Indirect Effect

### Parameter Estimate

Bias-Corrected Bootstrap:

```{r}
mediationFit_estimates_bca %>% 
  filter(label == "indirect")
```

Percentile Bootstrap:

```{r}
mediationFit_estimates_perc %>% 
  filter(label == "indirect")
```

### Effect Size {#sec-effectSizeMediation}

#### Standardized Estimate ($\beta$)

$$
\beta(ab) = ab \cdot \frac{SD_\text{X}}{SD_\text{Y}}
$$

```{r}
mediationFit_indirect <- mediationFit_estimates %>% 
  filter(label == "indirect") %>% 
  select(std.all) %>% 
  as.numeric

mediationFit_indirect
```

#### Proportion Mediated (*P*<sub>*M*</sub>) {#sec-proportionMediated}

$$
P_M = \frac{|ab|}{|c|} = \frac{|ab|}{|c'| + |ab|}
$$

Effect size: Proportion mediated (*P*<sub>*M*</sub>); i.e., the proportion of the total effect that is mediated; calculated by the magnitude of the indirect effect divided by the magnitude of the total effect:

```{r}
mediationFit_total <- mediationFit_estimates %>% 
  filter(label == "totalAbs") %>% 
  select(std.all) %>% 
  as.numeric

mediationFit_pm <- abs(mediationFit_indirect) / abs(mediationFit_total)
mediationFit_pm
```

In this case, the direct effect and indirect effect have opposite signs (negative and positive, respectively).
This is called *inconsistent mediation*, and can render the estimate of proportion mediated not a meaningful estimate of effect size (in this case, the estimate would exceed 1.0; Fairchild & McDaniel, 2017).
To address this, we can use the absolute value of the direct and indirect effects (in computation of the total effect), as recommended by MacKinnon et al. (2007).
When using the absolute value of the direct and indirect effects, the proportion mediated (*P*<sub>*M*</sub>) is `r petersenlab::apa(mediationFit_pm, decimals = 2)`, indicating that the mediator explained `r petersenlab::apa(mediationFit_pm * 100, decimals = 0)`% of the total effect.

#### Proportion of Variance in Y That is Explained by the Indirect Effect (*R*<sup>2</sup><sub>mediated</sub>) {#sec-rSquaredMediated}

Formulas from Lachowicz et al. (2018):

$$
\begin{aligned}
  R^2_\text{mediated} &= r^2_{\text{MY}} - (R^2_{\text{Y} \cdot \text{MX}} - r^2_{\text{XY}}) \\
  &= (\beta^2_{\text{YM} \cdot \text{X}} + \beta_{\text{YX} \cdot \text{M}} \cdot \beta_{\text{MX}}) ^2 - [\beta^2_{\text{YX}} + \beta^2_{\text{YM} \cdot \text{X}}(1 - \beta^2_{\text{MX}}) - \beta^2_{\text{YX}}]
\end{aligned}
$$

```{r}
rXY <- as.numeric(cor.test(
  ~ X + Y,
  data = mydata
)$estimate)

rMY <- as.numeric(cor.test(
  ~ M + Y,
  data = mydata
)$estimate)

RsquaredYmx <- summary(lm(
  Y ~ M + X,
  data = mydata))$r.squared

RsquaredMed1 <- (rMY^2) - (RsquaredYmx - (rXY^2))
RsquaredMed1

betaYMgivenX <- mediationFit_estimates %>% 
  filter(label == "b") %>% 
  select(std.all) %>% 
  as.numeric

betaYXgivenM <- mediationFit_estimates %>% 
  filter(label == "direct") %>% 
  select(std.all) %>% 
  as.numeric

betaMX <- mediationFit_estimates %>% 
  filter(label == "a") %>% 
  select(std.all) %>% 
  as.numeric

betaYX <- as.numeric(cor.test(
  ~ X + Y,
  data = mydata
)$estimate)

RsquaredMed2 <- ((betaYMgivenX + (betaYXgivenM * betaMX))^2) - ((betaYX^2) + (betaYMgivenX^2)*(1 - (betaMX^2)) - (betaYX^2))
RsquaredMed2
```

#### The Proportion of Variance in Y That is Accounted for Jointly by M and X (upsilon; $v$) {#sec-upsilon}

Formulas from Lachowicz et al. (2018):

$$
\begin{aligned}
  v &= (r_{\text{YM}} - \beta_{\text{MX}} \cdot \beta^2_{\text{YX} \cdot \text{M}}) ^ 2 - (R^2_{\text{Y} \cdot \text{MX}} - r^2_{\text{YX}})\\
  &= \beta^2_a \cdot \beta^2_b
\end{aligned}
$$

where $a$ is the $a$ path ($\beta^2_{\text{MX}}$), and $b$ is the $b$ path ($\beta^2_{\text{YM} \cdot \text{X}}$).

The estimate corrects for spurious correlation induced by the ordering of variables.

```{r}
upsilon1 <- ((rMY - (betaMX * (betaYXgivenM^2)))^2) - (RsquaredYmx - (rXY^2))
upsilon1

upsilon2 <- (betaYMgivenX^2) - (RsquaredYmx - (rXY^2))
upsilon2

upsilon3 <- mediationFit_indirect ^ 2
upsilon3

upsilon(
  x = mydata$X,
  mediator = mydata$M,
  dv = mydata$Y,
  bootstrap = FALSE
)
```

#### Ratio of the Indirect Effect Relative to Its Maximum Possible Value in the Data ($\kappa^2$) {#sec-kappaSquared}

$$
\kappa^2 = \frac{ab}{\text{MAX}(ab)}
$$

Kappa-squared ($\kappa^2$) is the ratio of the indirect effect relative to its maximum possible value in the data given the observed variability of X, Y, and M and their intercorrelations in the data.
This estimate is no longer recommended (Wen & Fan, 2015).

#### Other Effect Sizes {#sec-effectSizeMediationOther}

```{r}
mediation(
  x = mydata$X,
  mediator = mydata$M,
  dv = mydata$Y,
  bootstrap = FALSE
)
```

## Estimates of Model Fit

The model is saturated because it has as many estimated parameters as there are data points (i.e., in terms of means, variances, and covariances), so it has zero degrees of freedom.
Because the model is saturated, it has "perfect" fit.

```{r}
fitMeasures(
  mediationFit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr"))
```

## Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(mediationFit, type = "cor")
```

## Modification Indices

```{r}
modificationindices(mediationFit, sort. = TRUE)
```

## Internal Consistency Reliability

```{r}
compRelSEM(mediationFit)
```

## Path Diagram

```{r}
semPaths(
  mediationFit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

lavaanPlot::lavaanPlot(
  mediationFit,
  coefs = TRUE,
  #covs = TRUE,
  stand = TRUE)

lavaanPlot::lavaanPlot2(
  mediationFit,
  #stand = TRUE, # currently throws error; uncomment out when fixed: https://github.com/alishinski/lavaanPlot/issues/52
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(mediationFit)
```

# Moderation {#sec-moderation}

```{r}
states <- as.data.frame(state.x77)
names(states)[which(names(states) == "HS Grad")] <- "HS.Grad"
states$Income_rescaled <- states$Income/100
```

## Preparing the Predictors

Make sure to mean-center or orthogonalize predictors before computing the interaction term.

### Mean Center Predictors

```{r}
states$Illiteracy_centered <- scale(states$Illiteracy, scale = FALSE)
states$Murder_centered <- scale(states$Murder, scale = FALSE)
```

### Orthogonalized Predictors

Orthogonalizing is residual centering.

```{r}
states$interaction_notCentered <- states$Illiteracy * states$Murder

states$Illiteracy_orthogonalized <- resid(lm(
  data = states,
  interaction_notCentered ~ Illiteracy
))

states$Murder_orthogonalized <- resid(lm(
  data = states,
  interaction_notCentered ~ Murder
))
```

## Compute Interaction Term

```{r}
states$interaction <- states$Illiteracy_centered * states$Murder_centered # or: states$Illiteracy_orthogonalized * states$Murder_orthogonalized
```

## Model Syntax

```{r}
moderationModel <- '
  Income_rescaled ~ Illiteracy_centered + Murder_centered + interaction + HS.Grad
'
```

## Fit the Model

```{r}
moderationFit <- sem(
  moderationModel,
  data = states,
  missing = "ML",
  estimator = "MLR",
  std.lv = TRUE,
  fixed.x = FALSE)
```

## Summary Output

```{r}
summary(
  moderationFit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

## Estimates of Model Fit

The model is saturated because it has as many estimated parameters as there are data points (i.e., in terms of means, variances, and covariances), so it has zero degrees of freedom.
Because the model is saturated, it has "perfect" fit.

```{r}
fitMeasures(
  moderationFit,
  fit.measures = c(
    "chisq", "df", "pvalue",
    "baseline.chisq","baseline.df","baseline.pvalue",
    "rmsea", "cfi", "tli", "srmr"))
```

## Residuals of Observed vs. Model-Implied Correlation Matrix

```{r}
residuals(moderationFit, type = "cor")
```

## Modification Indices

```{r}
modificationindices(moderationFit, sort. = TRUE)
```

## Path Diagram

```{r}
semPaths(
  moderationFit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)

lavaanPlot::lavaanPlot(
  moderationFit,
  coefs = TRUE,
  #covs = TRUE,
  stand = TRUE)

lavaanPlot::lavaanPlot2(
  moderationFit,
  #stand = TRUE, # currently throws error; uncomment out when fixed: https://github.com/alishinski/lavaanPlot/issues/52
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(moderationFit)
```

## Interaction Plot {#sec-moderationInteractionPlot}

```{r}
# Created Model-Implied Predicted Data Object
modelImpliedPredictedData <- expand.grid(
  Illiteracy_factor = c("Low","Middle","High"),
  Murder_factor = c("Low","Middle","High"))

Illiteracy_mean <- mean(states$Illiteracy, na.rm = TRUE)
Illiteracy_sd <- sd(states$Illiteracy, na.rm = TRUE)

Murder_mean <- mean(states$Murder, na.rm = TRUE)
Murder_sd <- sd(states$Murder, na.rm = TRUE)

Illiteracy_centered_mean <- mean(states$Illiteracy_centered, na.rm = TRUE)
Illiteracy_centered_sd <- sd(states$Illiteracy_centered, na.rm = TRUE)

Murder_centered_mean <- mean(states$Murder_centered, na.rm = TRUE)
Murder_centered_sd <- sd(states$Murder_centered, na.rm = TRUE)

modelImpliedPredictedData <- modelImpliedPredictedData %>%
  mutate(
    Illiteracy = case_when(
      Illiteracy_factor == "Low" ~ Illiteracy_mean - Illiteracy_sd,
      Illiteracy_factor == "Middle" ~ Illiteracy_mean,
      Illiteracy_factor == "High" ~ Illiteracy_mean + Illiteracy_sd
    ),
    Illiteracy_centered = case_when(
      Illiteracy_factor == "Low" ~ Illiteracy_centered_mean - Illiteracy_centered_sd,
      Illiteracy_factor == "Middle" ~ Illiteracy_centered_mean,
      Illiteracy_factor == "High" ~ Illiteracy_centered_mean + Illiteracy_centered_sd
    ),
    Murder = case_when(
      Murder_factor == "Low" ~ Murder_mean - Murder_sd,
      Murder_factor == "Middle" ~ Murder_mean,
      Murder_factor == "High" ~ Murder_mean + Murder_sd
    ),
    Murder_centered = case_when(
      Murder_factor == "Low" ~ Murder_centered_mean - Murder_centered_sd,
      Murder_factor == "Middle" ~ Murder_centered_mean,
      Murder_factor == "High" ~ Murder_centered_mean + Murder_centered_sd
    ),
    interaction = Illiteracy_centered * Murder_centered,
    HS.Grad = mean(states$HS.Grad, na.rm = TRUE), # mean for covariates
    Income_rescaled = NA
  )

Murder_labels <- factor(
  modelImpliedPredictedData$Murder_factor,
  levels = c("High", "Middle", "Low"),
  labels = c("High (+1 SD)", "Middle (mean)", "Low (−1 SD)"))

modelImpliedPredictedData$Income_rescaled <- lavPredictY(
  moderationFit,
  newdata = modelImpliedPredictedData,
  ynames = "Income_rescaled"
) %>% 
  as.vector()

# Verify Computation Manually
moderationFit_parameters <- parameterEstimates(moderationFit)

moderationFit_parameters

intercept <- moderationFit_parameters[which(moderationFit_parameters$lhs == "Income_rescaled" & moderationFit_parameters$op == "~1"), "est"]
b_Illiteracy_centered <- moderationFit_parameters[which(moderationFit_parameters$lhs == "Income_rescaled" & moderationFit_parameters$rhs == "Illiteracy_centered"), "est"]
b_Murder_centered <- moderationFit_parameters[which(moderationFit_parameters$lhs == "Income_rescaled" & moderationFit_parameters$rhs == "Murder_centered"), "est"]
b_interaction <- moderationFit_parameters[which(moderationFit_parameters$lhs == "Income_rescaled" & moderationFit_parameters$rhs == "interaction"), "est"]
b_HS.Grad <- moderationFit_parameters[which(moderationFit_parameters$lhs == "Income_rescaled" & moderationFit_parameters$rhs == "HS.Grad"), "est"]

modelImpliedPredictedData <- modelImpliedPredictedData %>%
  mutate(
    Income_rescaled_calculatedManually = intercept + (b_Illiteracy_centered * Illiteracy_centered) + (b_Murder_centered * Murder_centered) + (b_interaction * interaction) + (b_HS.Grad * HS.Grad))

# Model-Implied Predicted Data
modelImpliedPredictedData

# Plot
ggplot(
  data = modelImpliedPredictedData,
  mapping = aes(
    x = Illiteracy,
    y = Income_rescaled,
    color = Murder_labels
  )
) +
  geom_line() +
  labs(color = "Murder")
```

## Simple Slopes and Regions of Significance {#sec-moderationRegionsOfSignificance}

<https://gabriellajg.github.io/EPSY-579-R-Cookbook-for-SEM/week6_1-lavaan-lab-4-mediated-moderation-moderated-mediation.html#step-5-johnson-neyman-interval> (archived at <https://perma.cc/6XR6-ZPSL>)

```{r}
# Find the min and max values of the moderator
Murder_centered_min <- min(modelImpliedPredictedData$Murder_centered, na.rm = TRUE)
Murder_centered_max <- max(modelImpliedPredictedData$Murder_centered, na.rm = TRUE)

Murder_centered_cutoff1 <- -1.5 # pick and titrate cutoff to help find the lower bound of the region of significance
Murder_centered_cutoff2 <- -1 # pick and titrate cutoff to help find the upper bound of the region of significance

Murder_centered_sd <- sd(modelImpliedPredictedData$Murder_centered, na.rm = TRUE)

Murder_centered_low <- mean(modelImpliedPredictedData$Murder_centered, na.rm = TRUE) - sd(modelImpliedPredictedData$Murder_centered, na.rm = TRUE)
Murder_centered_mean <- mean(modelImpliedPredictedData$Murder_centered, na.rm = TRUE)
Murder_centered_high <- mean(modelImpliedPredictedData$Murder_centered, na.rm = TRUE) + sd(modelImpliedPredictedData$Murder_centered, na.rm = TRUE)

# Extend the moderation model to compute the simple slopes and conditional effects at specific values of the moderator
moderationModelSimpleSlopes <- paste0('
  # Regression
  Income_rescaled ~ b1*Illiteracy_centered + b2*Murder_centered + b3*interaction + b4*HS.Grad
  
  # Simple Slopes
  SS_min := b1 + b3 * ', Murder_centered_min, '
  SS_cutoff1 := b1 + b3 * ', Murder_centered_cutoff1, '
  SS_cutoff2 := b1 + b3 * ', Murder_centered_cutoff2, '
  SS_low := b1 + b3 * ', Murder_centered_low, '
  SS_mean := b1 + b3 * ', Murder_centered_mean, '
  SS_high := b1 + b3 * ', Murder_centered_high, '
  SS_max := b1 + b3 * ', Murder_centered_max, '
')

# Fit the Model
set.seed(52242) # for reproducibility

moderationModelSimpleSlopes_fit <- sem(
  model = moderationModelSimpleSlopes, 
  data = states,
  missing = "ML",
  estimator = "ML",
  se = "bootstrap",
  bootstrap = 1000,
  std.lv = TRUE,
  fixed.x = FALSE)

summary(
  moderationModelSimpleSlopes_fit,
  #fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)

moderationModelSimpleSlopesFit_parameters <- parameterEstimates(
  moderationModelSimpleSlopes_fit,
  level = 0.95,
  boot.ci.type = "bca.simple")

moderationModelSimpleSlopesFit_parameters
```

A simple slope of the predictor on the outcome is considered significant at a given level of the moderator if the 95% confidence interval from the bootstrapped estimates of the simple slopes at that level of the moderator (i.e., [`ci.lower`,`ci.upper`]) does not include zero.
In this particular model, the predictor (`Illiteracy`) is not significant at any of the levels of the moderator (`Murder`), because the 95% confidence intervals of all simple slopes include zero, in this case, likely due to a small sample size ($N = 50$) and the resulting low power.

## Johnson-Neyman Plot {#sec-johnsonNeymanPlot}

As I noted above, the predictor is not significant at any levels of the moderator.
Nevertheless, I created a made up Johnson-Neyman plot by specifying the (fictitious) range of significance, for purposes of demonstration.
The band around the line indicates the 95% confidence interval of the simple slope of the predictor on the outcome as a function of different levels of the moderator.
In reality (unlike in this fictitious example), the regions of significance would only be regions where the 95% confidence interval of the simple slope does not include zero.

The standard error of the slope is the square root of the variance of the slope.
The forumula for computing the standard error of the slope is based on the formula for computing the variance of a weighted sum.

The slope of the predictor on the outcome at different levels of the moderator is calculated as (Jaccard & Turisi, 2003):

$$
\text{slope}_\text{predictor} = b_1 + b_3 \cdot Z
$$

The standard error of the slope of the predictor on the outcome at different levels of the moderator is calculated as (<https://stats.stackexchange.com/a/55973/20338>; archived at <https://perma.cc/V255-853Z>; Jaccard & Turisi, 2003):

$$
\begin{aligned}
SE(\text{slope}_\text{predictor}) &= \sqrt{Var(b_1) + Var(b_3) \cdot Z^2 + 2 \cdot Z \cdot Cov(b_1, b_3)} \\
SE(b_1 + b_3 \cdot Z) &=
\end{aligned}
$$

where:

- $b_1$ is the slope of the predictor on the outcome
- $b_3$ is the slope of the interaction term on the outcome
- $Z$ is the moderator

The variance of a weighted sum is:

$$
\begin{aligned}
Var(\text{slope}_\text{predictor}) &= Var(b_1) + Var(b_3) \cdot Z^2 + 2 \cdot Z \cdot Cov(b_1, b_3) \\
Var(b_1 + b_3 \cdot Z) &=
\end{aligned}
$$

The standard error is the square root of the variance.
The 95% confidence interval of the slope is $\pm$ `r qnorm(.975)` (i.e., `qnorm(.975)`) standard errors of the slope estimate.

```{r}
# Create a data frame for plotting
Murder_min <- min(states$Murder, na.rm = TRUE)
Murder_max <- max(states$Murder, na.rm = TRUE)

plot_data <- data.frame(
  Murder = seq(Murder_min, Murder_max, length.out = 10000)
)

plot_data$Murder_centered <- scale(plot_data$Murder, scale = FALSE)

# Calculate predicted slopes and confidence intervals
b1 <- moderationModelSimpleSlopesFit_parameters[which(moderationModelSimpleSlopesFit_parameters$label == "b1"), "est"]
b3 <- moderationModelSimpleSlopesFit_parameters[which(moderationModelSimpleSlopesFit_parameters$label == "b3"), "est"]

b1_se <- moderationModelSimpleSlopesFit_parameters[which(moderationModelSimpleSlopesFit_parameters$label == "b1"), "se"]
b3_se <- moderationModelSimpleSlopesFit_parameters[which(moderationModelSimpleSlopesFit_parameters$label == "b3"), "se"]

varianceCovarianceMatrix <- vcov(moderationFit)

b1_var <- varianceCovarianceMatrix["Income_rescaled~Illiteracy_centered","Income_rescaled~Illiteracy_centered"]
b3_var <- varianceCovarianceMatrix["interaction~~interaction","interaction~~interaction"]
cov_b1b3 <- varianceCovarianceMatrix["Income_rescaled~Illiteracy_centered","interaction~~interaction"]

#sqrt((b1_se^2) + ((b3_se^2) * plot_data$Murder_centered^2) + (2 * plot_data$Murder_centered * cov_b1b3))
#sqrt((b1_var) + ((b3_var) * plot_data$Murder_centered^2) + (2 * plot_data$Murder_centered * cov_b1b3))

plot_data$predicted_slopes <- b1 + b3 * plot_data$Murder_centered
plot_data$slope_se <- sqrt((b1_var) + ((b3_var) * plot_data$Murder_centered^2) + (2 * plot_data$Murder_centered * cov_b1b3))

# Calculated the 95% confidence interval around the simple slope
plot_data$lower_ci <- plot_data$predicted_slopes - qnorm(.975) * plot_data$slope_se
plot_data$upper_ci <- plot_data$predicted_slopes + qnorm(.975) * plot_data$slope_se

# Specify the significant range (based on the regions identified in the simple slopes analysis, see "Simple Slopes and Regions of Significance" section above)
plot_data$significant_slope <- FALSE
plot_data$significant_slope[which(plot_data$Murder_centered < -4.2 | plot_data$Murder_centered > 3.75)] <-TRUE # specify significant range

# Specify the significant region number (there are either 0, 1, or 2 significant regions; in such cases, there would be 1, 0 or 1 or 2, or 1 nonsignificant regions, respectively)--for instance, sig from 0-4, ns from 4-12, and sig from 12-16 would be 2 significant regions and 1 nonsignificant region
plot_data$significantRegionNumber <- NA
plot_data$significantRegionNumber[which(plot_data$Murder_centered < -4.2)] <- 1 # specify significant range 1
plot_data$significantRegionNumber[which(plot_data$Murder_centered > 3.75)] <- 2 # specify significant range 2

min(plot_data$Murder[which(plot_data$significant_slope == FALSE)])
max(plot_data$Murder[which(plot_data$significant_slope == FALSE)])

ggplot(plot_data, aes(x = Murder, y = predicted_slopes)) +
  geom_ribbon(
    data = plot_data %>% filter(significant_slope == FALSE),
    aes(ymin = lower_ci, ymax = upper_ci),
    fill = "#F8766D",
    alpha = 0.2) + 
  geom_ribbon(
    data = plot_data %>% filter(significantRegionNumber == 1),
    aes(ymin = lower_ci, ymax = upper_ci),
    fill = "#00BFC4",
    alpha = 0.2) + 
  geom_ribbon(
    data = plot_data %>% filter(significantRegionNumber == 2),
    aes(ymin = lower_ci, ymax = upper_ci),
    fill = "#00BFC4",
    alpha = 0.2) +
  geom_line(
    data = plot_data %>% filter(significant_slope == FALSE),
    aes(x = Murder, y = predicted_slopes),
    color = "#F8766D",
    linewidth = 2) +
  geom_line(
    data = plot_data %>% filter(significantRegionNumber == 1),
    aes(x = Murder, y = predicted_slopes),
    color = "#00BFC4",
    linewidth = 2) +
  geom_line(
    data = plot_data %>% filter(significantRegionNumber == 2),
    aes(x = Murder, y = predicted_slopes),
    color = "#00BFC4",
    linewidth = 2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = c(4.051215, 11.99938), linetype = 2, color = "#00BFC4") + # update based on numbers above
  labs(
    title = "Johnson-Neyman Plot",
    subtitle = "(blue = significant slope; pink = nonsignificant slope)",
    x = "Moderator (Murder)",
    y = "Simple Slope of Predictor (Illiteracy)") +
  theme_classic()
```

# Longitudinal Measurement Invariance {#sec-longitudinalMI}

The code examples below follow the approach suggested by Widaman et al. (2010).

Widaman, K. F., Ferrer, E., & Conger, R. D. (2010). Factorial invariance within longitudinal structural equation models: Measuring the same construct across time. *Child Development Perspectives*, *4*(1), 10–18. <https://doi.org/10.1111/j.1750-8606.2009.00110.x>

For the longitudinal measurement invariance models, we use simulated data from <https://mycourses.aalto.fi/mod/assign/view.php?id=1203047> (archived at <https://perma.cc/QTL2-ZHX2>).
As noted [here](https://mycourses.aalto.fi/mod/assign/view.php?id=1203047) (archived at <https://perma.cc/QTL2-ZHX2>), "The data are in wide format where the variables are coded as [VARIABLE NAME][time index][indicator index]."

## Configural Invariance

Evaluates whether there are the same number of latent factors across time and whether the indicators load onto the same latent factor(s) across time.

1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. For each latent construct, constrain the first indicator's factor loading to be the same across time
1. For each latent construct, constrain the first indicator's intercept to be the same across time

### Model Syntax

```{r}
configuralInvariance_syntax <- '
  # Factor Loadings
  jobsat_1 =~ NA*loadj1*JOBSAT11 + JOBSAT12 + JOBSAT13
  jobsat_2 =~ NA*loadj1*JOBSAT21 + JOBSAT22 + JOBSAT23
  jobsat_3 =~ NA*loadj1*JOBSAT31 + JOBSAT32 + JOBSAT33
  
  ready_1 =~ NA*loadr1*READY11 + READY12 + READY13
  ready_2 =~ NA*loadr1*READY21 + READY22 + READY23
  ready_3 =~ NA*loadr1*READY31 + READY32 + READY33
  
  commit_1 =~ NA*loadc1*COMMIT11 + COMMIT12 + COMMIT13
  commit_2 =~ NA*loadc1*COMMIT21 + COMMIT22 + COMMIT23
  commit_3 =~ NA*loadc1*COMMIT31 + COMMIT32 + COMMIT33
  
  # Factor Identification: Standardize Factors at T1
  
  ## Fix Factor Means at T1 to Zero
  jobsat_1 ~ 0*1
  ready_1 ~ 0*1
  commit_1 ~ 0*1
  
  ## Fix Factor Variances at T1 to One
  jobsat_1 ~~ 1*jobsat_1
  ready_1 ~~ 1*ready_1
  commit_1 ~~ 1*commit_1
  
  # Freely Estimate Factor Means at T2 and T3 (relative to T1)
  jobsat_2 ~ 1
  jobsat_3 ~ 1
  
  ready_2 ~ 1
  ready_3 ~ 1
  
  commit_2 ~ 1
  commit_3 ~ 1
  
  # Freely Estimate Factor Variances at T2 and T3 (relative to T1)
  jobsat_2 ~~ jobsat_2
  jobsat_3 ~~ jobsat_3
  
  ready_2 ~~ ready_2
  ready_3 ~~ ready_3
  
  commit_2 ~~ commit_2
  commit_3 ~~ commit_3
  
  # Fix Intercepts of Indicator 1 Across Time
  JOBSAT11 ~ intj1*1
  JOBSAT21 ~ intj1*1
  JOBSAT31 ~ intj1*1
  
  READY11 ~ intr1*1
  READY21 ~ intr1*1
  READY31 ~ intr1*1
  
  COMMIT11 ~ intc1*1
  COMMIT21 ~ intc1*1
  COMMIT31 ~ intc1*1
  
  # Free Intercepts of Remaining Manifest Variables
  JOBSAT12 ~ 1
  JOBSAT13 ~ 1
  JOBSAT22 ~ 1
  JOBSAT23 ~ 1
  JOBSAT32 ~ 1
  JOBSAT33 ~ 1
  
  READY12 ~ 1
  READY13 ~ 1
  READY22 ~ 1
  READY23 ~ 1
  READY32 ~ 1
  READY33 ~ 1
  
  COMMIT12 ~ 1
  COMMIT13 ~ 1
  COMMIT22 ~ 1
  COMMIT23 ~ 1
  COMMIT32 ~ 1
  COMMIT33 ~ 1
  
  # Estimate Residual Variances of Manifest Variables
  JOBSAT11 ~~ JOBSAT11
  JOBSAT12 ~~ JOBSAT12
  JOBSAT13 ~~ JOBSAT13
  JOBSAT21 ~~ JOBSAT21
  JOBSAT22 ~~ JOBSAT22
  JOBSAT23 ~~ JOBSAT23
  JOBSAT31 ~~ JOBSAT31
  JOBSAT32 ~~ JOBSAT32
  JOBSAT33 ~~ JOBSAT33
  
  READY11 ~~ READY11
  READY12 ~~ READY12
  READY13 ~~ READY13
  READY21 ~~ READY21
  READY22 ~~ READY22
  READY23 ~~ READY23
  READY31 ~~ READY31
  READY32 ~~ READY32
  READY33 ~~ READY33
  
  COMMIT11 ~~ COMMIT11
  COMMIT12 ~~ COMMIT12
  COMMIT13 ~~ COMMIT13
  COMMIT21 ~~ COMMIT21
  COMMIT22 ~~ COMMIT22
  COMMIT23 ~~ COMMIT23
  COMMIT31 ~~ COMMIT31
  COMMIT32 ~~ COMMIT32
  COMMIT33 ~~ COMMIT33
'
```

### Fit Model

```{r}
configuralInvariance_fit <- cfa(
  configuralInvariance_syntax,
  data = longitudinalMI,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  #std.lv = TRUE,
  fixed.x = FALSE)
```

### Model Summary

```{r}
summary(
  configuralInvariance_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Modification Indices

```{r}
modificationindices(
  configuralInvariance_fit,
  sort. = TRUE)
```

### Path Diagram

```{r}
#| fig-cap: "Path Diagram"

lavaanPlot::lavaanPlot2(
  configuralInvariance_fit,
  stand = TRUE,
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(configuralInvariance_fit)
```

## Configural Invariance With Correlated Residuals Within-Indicator Across Time

1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. For each latent construct, constrain the first indicator's factor loading to be the same across time
1. For each latent construct, constrain the first indicator's intercept to be the same across time
1. **Allow within-indicator residuals to covary across time**

### Model Syntax

```{r}
configuralInvarianceCorrelatedResiduals_syntax <- '
  # Factor Loadings
  jobsat_1 =~ NA*loadj1*JOBSAT11 + JOBSAT12 + JOBSAT13
  jobsat_2 =~ NA*loadj1*JOBSAT21 + JOBSAT22 + JOBSAT23
  jobsat_3 =~ NA*loadj1*JOBSAT31 + JOBSAT32 + JOBSAT33
  
  ready_1 =~ NA*loadr1*READY11 + READY12 + READY13
  ready_2 =~ NA*loadr1*READY21 + READY22 + READY23
  ready_3 =~ NA*loadr1*READY31 + READY32 + READY33
  
  commit_1 =~ NA*loadc1*COMMIT11 + COMMIT12 + COMMIT13
  commit_2 =~ NA*loadc1*COMMIT21 + COMMIT22 + COMMIT23
  commit_3 =~ NA*loadc1*COMMIT31 + COMMIT32 + COMMIT33
  
  # Factor Identification: Standardize Factors at T1
  
  ## Fix Factor Means at T1 to Zero
  jobsat_1 ~ 0*1
  ready_1 ~ 0*1
  commit_1 ~ 0*1
  
  ## Fix Factor Variances at T1 to One
  jobsat_1 ~~ 1*jobsat_1
  ready_1 ~~ 1*ready_1
  commit_1 ~~ 1*commit_1
  
  # Freely Estimate Factor Means at T2 and T3 (relative to T1)
  jobsat_2 ~ 1
  jobsat_3 ~ 1
  
  ready_2 ~ 1
  ready_3 ~ 1
  
  commit_2 ~ 1
  commit_3 ~ 1
  
  # Freely Estimate Factor Variances at T2 and T3 (relative to T1)
  jobsat_2 ~~ jobsat_2
  jobsat_3 ~~ jobsat_3
  
  ready_2 ~~ ready_2
  ready_3 ~~ ready_3
  
  commit_2 ~~ commit_2
  commit_3 ~~ commit_3
  
  # Fix Intercepts of Indicator 1 Across Time
  JOBSAT11 ~ intj1*1
  JOBSAT21 ~ intj1*1
  JOBSAT31 ~ intj1*1
  
  READY11 ~ intr1*1
  READY21 ~ intr1*1
  READY31 ~ intr1*1
  
  COMMIT11 ~ intc1*1
  COMMIT21 ~ intc1*1
  COMMIT31 ~ intc1*1
  
  # Free Intercepts of Remaining Manifest Variables
  JOBSAT12 ~ 1
  JOBSAT13 ~ 1
  JOBSAT22 ~ 1
  JOBSAT23 ~ 1
  JOBSAT32 ~ 1
  JOBSAT33 ~ 1
  
  READY12 ~ 1
  READY13 ~ 1
  READY22 ~ 1
  READY23 ~ 1
  READY32 ~ 1
  READY33 ~ 1
  
  COMMIT12 ~ 1
  COMMIT13 ~ 1
  COMMIT22 ~ 1
  COMMIT23 ~ 1
  COMMIT32 ~ 1
  COMMIT33 ~ 1
  
  # Estimate Residual Variances of Manifest Variables
  JOBSAT11 ~~ JOBSAT11
  JOBSAT12 ~~ JOBSAT12
  JOBSAT13 ~~ JOBSAT13
  JOBSAT21 ~~ JOBSAT21
  JOBSAT22 ~~ JOBSAT22
  JOBSAT23 ~~ JOBSAT23
  JOBSAT31 ~~ JOBSAT31
  JOBSAT32 ~~ JOBSAT32
  JOBSAT33 ~~ JOBSAT33
  
  READY11 ~~ READY11
  READY12 ~~ READY12
  READY13 ~~ READY13
  READY21 ~~ READY21
  READY22 ~~ READY22
  READY23 ~~ READY23
  READY31 ~~ READY31
  READY32 ~~ READY32
  READY33 ~~ READY33
  
  COMMIT11 ~~ COMMIT11
  COMMIT12 ~~ COMMIT12
  COMMIT13 ~~ COMMIT13
  COMMIT21 ~~ COMMIT21
  COMMIT22 ~~ COMMIT22
  COMMIT23 ~~ COMMIT23
  COMMIT31 ~~ COMMIT31
  COMMIT32 ~~ COMMIT32
  COMMIT33 ~~ COMMIT33
  
  # Residual Covariances Within Indicator Across Time
  JOBSAT11 ~~ JOBSAT21
  JOBSAT21 ~~ JOBSAT31
  JOBSAT11 ~~ JOBSAT31
  
  JOBSAT12 ~~ JOBSAT22
  JOBSAT22 ~~ JOBSAT32
  JOBSAT12 ~~ JOBSAT32
  
  JOBSAT13 ~~ JOBSAT23
  JOBSAT23 ~~ JOBSAT33
  JOBSAT13 ~~ JOBSAT33
  
  READY11 ~~ READY21
  READY21 ~~ READY31
  READY11 ~~ READY31
  
  READY12 ~~ READY22
  READY22 ~~ READY32
  READY12 ~~ READY32
  
  READY13 ~~ READY23
  READY23 ~~ READY33
  READY13 ~~ READY33
  
  COMMIT11 ~~ COMMIT21
  COMMIT21 ~~ COMMIT31
  COMMIT11 ~~ COMMIT31
  
  COMMIT12 ~~ COMMIT22
  COMMIT22 ~~ COMMIT32
  COMMIT12 ~~ COMMIT32
  
  COMMIT13 ~~ COMMIT23
  COMMIT23 ~~ COMMIT33
  COMMIT13 ~~ COMMIT33
'
```

### Fit Model

```{r}
configuralInvarianceCorrelatedResiduals_fit <- cfa(
  configuralInvarianceCorrelatedResiduals_syntax,
  data = longitudinalMI,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  #std.lv = TRUE,
  fixed.x = FALSE)
```

### Model Summary

```{r}
summary(
  configuralInvarianceCorrelatedResiduals_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Modification Indices

```{r}
modificationindices(
  configuralInvarianceCorrelatedResiduals_fit,
  sort. = TRUE)
```

### Path Diagram

```{r}
#| fig-cap: "Path Diagram"

lavaanPlot::lavaanPlot2(
  configuralInvarianceCorrelatedResiduals_fit,
  stand = TRUE,
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(configuralInvarianceCorrelatedResiduals_fit)
```

### Compare Models

```{r}
anova(
  configuralInvariance_fit,
  configuralInvarianceCorrelatedResiduals_fit
)
```

## Metric ("Weak Factorial") Invariance

Evaluates whether the items' factor loadings are the same across time.

1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. For each latent construct, constrain the first indicator's intercept to be the same across time
1. Allow within-indicator residuals to covary across time
1. **For each indicator, constrain its factor loading to be the same across time**

### Model Syntax

```{r}
metricInvariance_syntax <- '
  # Factor Loadings
  jobsat_1 =~ NA*loadj1*JOBSAT11 + loadj2*JOBSAT12 + loadj3*JOBSAT13
  jobsat_2 =~ NA*loadj1*JOBSAT21 + loadj2*JOBSAT22 + loadj3*JOBSAT23
  jobsat_3 =~ NA*loadj1*JOBSAT31 + loadj2*JOBSAT32 + loadj3*JOBSAT33
  
  ready_1 =~ NA*loadr1*READY11 + loadr2*READY12 + loadr3*READY13
  ready_2 =~ NA*loadr1*READY21 + loadr2*READY22 + loadr3*READY23
  ready_3 =~ NA*loadr1*READY31 + loadr2*READY32 + loadr3*READY33
  
  commit_1 =~ NA*loadc1*COMMIT11 + loadc2*COMMIT12 + loadc3*COMMIT13
  commit_2 =~ NA*loadc1*COMMIT21 + loadc2*COMMIT22 + loadc3*COMMIT23
  commit_3 =~ NA*loadc1*COMMIT31 + loadc2*COMMIT32 + loadc3*COMMIT33
  
  # Factor Identification: Standardize Factors at T1
  
  ## Fix Factor Means at T1 to Zero
  jobsat_1 ~ 0*1
  ready_1 ~ 0*1
  commit_1 ~ 0*1
  
  ## Fix Factor Variances at T1 to One
  jobsat_1 ~~ 1*jobsat_1
  ready_1 ~~ 1*ready_1
  commit_1 ~~ 1*commit_1
  
  # Freely Estimate Factor Means at T2 and T3 (relative to T1)
  jobsat_2 ~ 1
  jobsat_3 ~ 1
  
  ready_2 ~ 1
  ready_3 ~ 1
  
  commit_2 ~ 1
  commit_3 ~ 1
  
  # Freely Estimate Factor Variances at T2 and T3 (relative to T1)
  jobsat_2 ~~ jobsat_2
  jobsat_3 ~~ jobsat_3
  
  ready_2 ~~ ready_2
  ready_3 ~~ ready_3
  
  commit_2 ~~ commit_2
  commit_3 ~~ commit_3
  
  # Fix Intercepts of Indicator 1 Across Time
  JOBSAT11 ~ intj1*1
  JOBSAT21 ~ intj1*1
  JOBSAT31 ~ intj1*1
  
  READY11 ~ intr1*1
  READY21 ~ intr1*1
  READY31 ~ intr1*1
  
  COMMIT11 ~ intc1*1
  COMMIT21 ~ intc1*1
  COMMIT31 ~ intc1*1
  
  # Free Intercepts of Remaining Manifest Variables
  JOBSAT12 ~ 1
  JOBSAT13 ~ 1
  JOBSAT22 ~ 1
  JOBSAT23 ~ 1
  JOBSAT32 ~ 1
  JOBSAT33 ~ 1
  
  READY12 ~ 1
  READY13 ~ 1
  READY22 ~ 1
  READY23 ~ 1
  READY32 ~ 1
  READY33 ~ 1
  
  COMMIT12 ~ 1
  COMMIT13 ~ 1
  COMMIT22 ~ 1
  COMMIT23 ~ 1
  COMMIT32 ~ 1
  COMMIT33 ~ 1
  
  # Estimate Residual Variances of Manifest Variables
  JOBSAT11 ~~ JOBSAT11
  JOBSAT12 ~~ JOBSAT12
  JOBSAT13 ~~ JOBSAT13
  JOBSAT21 ~~ JOBSAT21
  JOBSAT22 ~~ JOBSAT22
  JOBSAT23 ~~ JOBSAT23
  JOBSAT31 ~~ JOBSAT31
  JOBSAT32 ~~ JOBSAT32
  JOBSAT33 ~~ JOBSAT33
  
  READY11 ~~ READY11
  READY12 ~~ READY12
  READY13 ~~ READY13
  READY21 ~~ READY21
  READY22 ~~ READY22
  READY23 ~~ READY23
  READY31 ~~ READY31
  READY32 ~~ READY32
  READY33 ~~ READY33
  
  COMMIT11 ~~ COMMIT11
  COMMIT12 ~~ COMMIT12
  COMMIT13 ~~ COMMIT13
  COMMIT21 ~~ COMMIT21
  COMMIT22 ~~ COMMIT22
  COMMIT23 ~~ COMMIT23
  COMMIT31 ~~ COMMIT31
  COMMIT32 ~~ COMMIT32
  COMMIT33 ~~ COMMIT33
  
  # Residual Covariances Within Indicator Across Time
  JOBSAT11 ~~ JOBSAT21
  JOBSAT21 ~~ JOBSAT31
  JOBSAT11 ~~ JOBSAT31
  
  JOBSAT12 ~~ JOBSAT22
  JOBSAT22 ~~ JOBSAT32
  JOBSAT12 ~~ JOBSAT32
  
  JOBSAT13 ~~ JOBSAT23
  JOBSAT23 ~~ JOBSAT33
  JOBSAT13 ~~ JOBSAT33
  
  READY11 ~~ READY21
  READY21 ~~ READY31
  READY11 ~~ READY31
  
  READY12 ~~ READY22
  READY22 ~~ READY32
  READY12 ~~ READY32
  
  READY13 ~~ READY23
  READY23 ~~ READY33
  READY13 ~~ READY33
  
  COMMIT11 ~~ COMMIT21
  COMMIT21 ~~ COMMIT31
  COMMIT11 ~~ COMMIT31
  
  COMMIT12 ~~ COMMIT22
  COMMIT22 ~~ COMMIT32
  COMMIT12 ~~ COMMIT32
  
  COMMIT13 ~~ COMMIT23
  COMMIT23 ~~ COMMIT33
  COMMIT13 ~~ COMMIT33
'
```

### Fit Model

```{r}
metricInvariance_fit <- cfa(
  metricInvariance_syntax,
  data = longitudinalMI,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  #std.lv = TRUE,
  fixed.x = FALSE)
```

### Model Summary

```{r}
summary(
  metricInvariance_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Modification Indices

```{r}
modificationindices(
  metricInvariance_fit,
  sort. = TRUE)
```

### Path Diagram

```{r}
#| fig-cap: "Path Diagram"

lavaanPlot::lavaanPlot2(
  metricInvariance_fit,
  stand = TRUE,
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(metricInvariance_fit)
```

### Compare Models

```{r}
anova(
  configuralInvarianceCorrelatedResiduals_fit,
  metricInvariance_fit
)
```

## Scalar ("Strong Factorial") Invariance

Evaluates whether the items' intercepts are the same across time.

1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. Allow within-indicator residuals to covary across time
1. For each indicator, constrain its factor loading to be the same across time
1. **For each indicator, constrain its intercept to be the same across time**

### Model Syntax

```{r}
scalarInvariance_syntax <- '
  # Factor Loadings
  jobsat_1 =~ NA*loadj1*JOBSAT11 + loadj2*JOBSAT12 + loadj3*JOBSAT13
  jobsat_2 =~ NA*loadj1*JOBSAT21 + loadj2*JOBSAT22 + loadj3*JOBSAT23
  jobsat_3 =~ NA*loadj1*JOBSAT31 + loadj2*JOBSAT32 + loadj3*JOBSAT33
  
  ready_1 =~ NA*loadr1*READY11 + loadr2*READY12 + loadr3*READY13
  ready_2 =~ NA*loadr1*READY21 + loadr2*READY22 + loadr3*READY23
  ready_3 =~ NA*loadr1*READY31 + loadr2*READY32 + loadr3*READY33
  
  commit_1 =~ NA*loadc1*COMMIT11 + loadc2*COMMIT12 + loadc3*COMMIT13
  commit_2 =~ NA*loadc1*COMMIT21 + loadc2*COMMIT22 + loadc3*COMMIT23
  commit_3 =~ NA*loadc1*COMMIT31 + loadc2*COMMIT32 + loadc3*COMMIT33
  
  # Factor Identification: Standardize Factors at T1
  
  ## Fix Factor Means at T1 to Zero
  jobsat_1 ~ 0*1
  ready_1 ~ 0*1
  commit_1 ~ 0*1
  
  ## Fix Factor Variances at T1 to One
  jobsat_1 ~~ 1*jobsat_1
  ready_1 ~~ 1*ready_1
  commit_1 ~~ 1*commit_1
  
  # Freely Estimate Factor Means at T2 and T3 (relative to T1)
  jobsat_2 ~ 1
  jobsat_3 ~ 1
  
  ready_2 ~ 1
  ready_3 ~ 1
  
  commit_2 ~ 1
  commit_3 ~ 1
  
  # Freely Estimate Factor Variances at T2 and T3 (relative to T1)
  jobsat_2 ~~ jobsat_2
  jobsat_3 ~~ jobsat_3
  
  ready_2 ~~ ready_2
  ready_3 ~~ ready_3
  
  commit_2 ~~ commit_2
  commit_3 ~~ commit_3
  
  # Fix Intercepts of Indicators Across Time
  JOBSAT11 ~ intj1*1
  JOBSAT21 ~ intj1*1
  JOBSAT31 ~ intj1*1
  JOBSAT12 ~ intj2*1
  JOBSAT22 ~ intj2*1
  JOBSAT32 ~ intj2*1
  JOBSAT13 ~ intj3*1
  JOBSAT23 ~ intj3*1
  JOBSAT33 ~ intj3*1
  
  READY11 ~ intr1*1
  READY21 ~ intr1*1
  READY31 ~ intr1*1
  READY12 ~ intr2*1
  READY22 ~ intr2*1
  READY32 ~ intr2*1
  READY13 ~ intr3*1
  READY23 ~ intr3*1
  READY33 ~ intr3*1
  
  COMMIT11 ~ intc1*1
  COMMIT21 ~ intc1*1
  COMMIT31 ~ intc1*1
  COMMIT12 ~ intc2*1
  COMMIT22 ~ intc2*1
  COMMIT32 ~ intc2*1
  COMMIT13 ~ intc3*1
  COMMIT23 ~ intc3*1
  COMMIT33 ~ intc3*1
  
  # Estimate Residual Variances of Manifest Variables
  JOBSAT11 ~~ JOBSAT11
  JOBSAT12 ~~ JOBSAT12
  JOBSAT13 ~~ JOBSAT13
  JOBSAT21 ~~ JOBSAT21
  JOBSAT22 ~~ JOBSAT22
  JOBSAT23 ~~ JOBSAT23
  JOBSAT31 ~~ JOBSAT31
  JOBSAT32 ~~ JOBSAT32
  JOBSAT33 ~~ JOBSAT33
  
  READY11 ~~ READY11
  READY12 ~~ READY12
  READY13 ~~ READY13
  READY21 ~~ READY21
  READY22 ~~ READY22
  READY23 ~~ READY23
  READY31 ~~ READY31
  READY32 ~~ READY32
  READY33 ~~ READY33
  
  COMMIT11 ~~ COMMIT11
  COMMIT12 ~~ COMMIT12
  COMMIT13 ~~ COMMIT13
  COMMIT21 ~~ COMMIT21
  COMMIT22 ~~ COMMIT22
  COMMIT23 ~~ COMMIT23
  COMMIT31 ~~ COMMIT31
  COMMIT32 ~~ COMMIT32
  COMMIT33 ~~ COMMIT33
  
  # Residual Covariances Within Indicator Across Time
  JOBSAT11 ~~ JOBSAT21
  JOBSAT21 ~~ JOBSAT31
  JOBSAT11 ~~ JOBSAT31
  
  JOBSAT12 ~~ JOBSAT22
  JOBSAT22 ~~ JOBSAT32
  JOBSAT12 ~~ JOBSAT32
  
  JOBSAT13 ~~ JOBSAT23
  JOBSAT23 ~~ JOBSAT33
  JOBSAT13 ~~ JOBSAT33
  
  READY11 ~~ READY21
  READY21 ~~ READY31
  READY11 ~~ READY31
  
  READY12 ~~ READY22
  READY22 ~~ READY32
  READY12 ~~ READY32
  
  READY13 ~~ READY23
  READY23 ~~ READY33
  READY13 ~~ READY33
  
  COMMIT11 ~~ COMMIT21
  COMMIT21 ~~ COMMIT31
  COMMIT11 ~~ COMMIT31
  
  COMMIT12 ~~ COMMIT22
  COMMIT22 ~~ COMMIT32
  COMMIT12 ~~ COMMIT32
  
  COMMIT13 ~~ COMMIT23
  COMMIT23 ~~ COMMIT33
  COMMIT13 ~~ COMMIT33
'
```

### Fit Model

```{r}
scalarInvariance_fit <- cfa(
  scalarInvariance_syntax,
  data = longitudinalMI,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  #std.lv = TRUE,
  fixed.x = FALSE)
```

### Model Summary

```{r}
summary(
  scalarInvariance_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Modification Indices

```{r}
modificationindices(
  scalarInvariance_fit,
  sort. = TRUE)
```

### Path Diagram

```{r}
#| fig-cap: "Path Diagram"

lavaanPlot::lavaanPlot2(
  scalarInvariance_fit,
  stand = TRUE,
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(scalarInvariance_fit)
```

### Compare Models

```{r}
anova(
  metricInvariance_fit,
  scalarInvariance_fit
)
```

## Residual ("Strict Factorial") Invariance

Evaluates whether the items' residual variances are the same across time.

1. Standardize the latent factor(s) at T1 (i.e., fix the mean to zero and the variance to one)
1. Allow within-indicator residuals to covary across time
1. For each indicator, constrain its factor loading to be the same across time
1. For each indicator, constrain its intercept to be the same across time
1. **For each indicator, constrain its residual variance to be the same across time**

### Model Syntax

```{r}
residualInvariance_syntax <- '
  # Factor Loadings
  jobsat_1 =~ NA*loadj1*JOBSAT11 + loadj2*JOBSAT12 + loadj3*JOBSAT13
  jobsat_2 =~ NA*loadj1*JOBSAT21 + loadj2*JOBSAT22 + loadj3*JOBSAT23
  jobsat_3 =~ NA*loadj1*JOBSAT31 + loadj2*JOBSAT32 + loadj3*JOBSAT33
  
  ready_1 =~ NA*loadr1*READY11 + loadr2*READY12 + loadr3*READY13
  ready_2 =~ NA*loadr1*READY21 + loadr2*READY22 + loadr3*READY23
  ready_3 =~ NA*loadr1*READY31 + loadr2*READY32 + loadr3*READY33
  
  commit_1 =~ NA*loadc1*COMMIT11 + loadc2*COMMIT12 + loadc3*COMMIT13
  commit_2 =~ NA*loadc1*COMMIT21 + loadc2*COMMIT22 + loadc3*COMMIT23
  commit_3 =~ NA*loadc1*COMMIT31 + loadc2*COMMIT32 + loadc3*COMMIT33
  
  # Factor Identification: Standardize Factors at T1
  
  ## Fix Factor Means at T1 to Zero
  jobsat_1 ~ 0*1
  ready_1 ~ 0*1
  commit_1 ~ 0*1
  
  ## Fix Factor Variances at T1 to One
  jobsat_1 ~~ 1*jobsat_1
  ready_1 ~~ 1*ready_1
  commit_1 ~~ 1*commit_1
  
  # Freely Estimate Factor Means at T2 and T3 (relative to T1)
  jobsat_2 ~ 1
  jobsat_3 ~ 1
  
  ready_2 ~ 1
  ready_3 ~ 1
  
  commit_2 ~ 1
  commit_3 ~ 1
  
  # Freely Estimate Factor Variances at T2 and T3 (relative to T1)
  jobsat_2 ~~ jobsat_2
  jobsat_3 ~~ jobsat_3
  
  ready_2 ~~ ready_2
  ready_3 ~~ ready_3
  
  commit_2 ~~ commit_2
  commit_3 ~~ commit_3
  
  # Constrain Intercepts of Indicators Across Time
  JOBSAT11 ~ intj1*1
  JOBSAT21 ~ intj1*1
  JOBSAT31 ~ intj1*1
  JOBSAT12 ~ intj2*1
  JOBSAT22 ~ intj2*1
  JOBSAT32 ~ intj2*1
  JOBSAT13 ~ intj3*1
  JOBSAT23 ~ intj3*1
  JOBSAT33 ~ intj3*1
  
  READY11 ~ intr1*1
  READY21 ~ intr1*1
  READY31 ~ intr1*1
  READY12 ~ intr2*1
  READY22 ~ intr2*1
  READY32 ~ intr2*1
  READY13 ~ intr3*1
  READY23 ~ intr3*1
  READY33 ~ intr3*1
  
  COMMIT11 ~ intc1*1
  COMMIT21 ~ intc1*1
  COMMIT31 ~ intc1*1
  COMMIT12 ~ intc2*1
  COMMIT22 ~ intc2*1
  COMMIT32 ~ intc2*1
  COMMIT13 ~ intc3*1
  COMMIT23 ~ intc3*1
  COMMIT33 ~ intc3*1
  
  # Constrain Residual Variances of Indicators Across Time
  JOBSAT11 ~~ resj1*JOBSAT11
  JOBSAT21 ~~ resj1*JOBSAT21
  JOBSAT31 ~~ resj1*JOBSAT31
  JOBSAT12 ~~ resj2*JOBSAT12
  JOBSAT22 ~~ resj2*JOBSAT22
  JOBSAT32 ~~ resj2*JOBSAT32
  JOBSAT13 ~~ resj3*JOBSAT13
  JOBSAT23 ~~ resj3*JOBSAT23
  JOBSAT33 ~~ resj3*JOBSAT33
  
  READY11 ~~ resr1*READY11
  READY21 ~~ resr1*READY21
  READY31 ~~ resr1*READY31
  READY12 ~~ resr2*READY12
  READY22 ~~ resr2*READY22
  READY32 ~~ resr2*READY32
  READY13 ~~ resr3*READY13
  READY23 ~~ resr3*READY23
  READY33 ~~ resr3*READY33
  
  COMMIT11 ~~ resc1*COMMIT11
  COMMIT21 ~~ resc1*COMMIT21
  COMMIT31 ~~ resc1*COMMIT31
  COMMIT12 ~~ resc2*COMMIT12
  COMMIT22 ~~ resc2*COMMIT22
  COMMIT32 ~~ resc2*COMMIT32
  COMMIT13 ~~ resc3*COMMIT13
  COMMIT23 ~~ resc3*COMMIT23
  COMMIT33 ~~ resc3*COMMIT33
  
  # Residual Covariances Within Indicator Across Time
  JOBSAT11 ~~ JOBSAT21
  JOBSAT21 ~~ JOBSAT31
  JOBSAT11 ~~ JOBSAT31
  
  JOBSAT12 ~~ JOBSAT22
  JOBSAT22 ~~ JOBSAT32
  JOBSAT12 ~~ JOBSAT32
  
  JOBSAT13 ~~ JOBSAT23
  JOBSAT23 ~~ JOBSAT33
  JOBSAT13 ~~ JOBSAT33
  
  READY11 ~~ READY21
  READY21 ~~ READY31
  READY11 ~~ READY31
  
  READY12 ~~ READY22
  READY22 ~~ READY32
  READY12 ~~ READY32
  
  READY13 ~~ READY23
  READY23 ~~ READY33
  READY13 ~~ READY33
  
  COMMIT11 ~~ COMMIT21
  COMMIT21 ~~ COMMIT31
  COMMIT11 ~~ COMMIT31
  
  COMMIT12 ~~ COMMIT22
  COMMIT22 ~~ COMMIT32
  COMMIT12 ~~ COMMIT32
  
  COMMIT13 ~~ COMMIT23
  COMMIT23 ~~ COMMIT33
  COMMIT13 ~~ COMMIT33
'
```

### Fit Model

```{r}
residualInvariance_fit <- cfa(
  residualInvariance_syntax,
  data = longitudinalMI,
  missing = "ML",
  estimator = "MLR",
  meanstructure = TRUE,
  #std.lv = TRUE,
  fixed.x = FALSE)
```

### Model Summary

```{r}
summary(
  residualInvariance_fit,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE)
```

### Modification Indices

```{r}
modificationindices(
  residualInvariance_fit,
  sort. = TRUE)
```

### Path Diagram

```{r}
#| fig-cap: "Path Diagram"

lavaanPlot::lavaanPlot2(
  residualInvariance_fit,
  stand = TRUE,
  coef_labels = TRUE)
```

To generate an interactive/modifiable path diagram, you can use the following syntax:

```{r}
#| eval: false

lavaangui::plot_lavaan(residualInvariance_fit)
```

### Compare Models

```{r}
anova(
  scalarInvariance_fit,
  residualInvariance_fit
)

round(cbind(
  configural = inspect(configuralInvarianceCorrelatedResiduals_fit, "fit.measures"), 
  weak = inspect(metricInvariance_fit, "fit.measures"),
  strong = inspect(scalarInvariance_fit, "fit.measures"),
  strict = inspect(residualInvariance_fit, "fit.measures")), 3)
```

# Bounded Estimation with Random Starts {#sec-boundedEstimationRandomStarts}

For more info, see De Jonckere and Rosseel (2022, 2025).

```{r}
HS.model <- '
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
'
```

```{r}
#| output: false

fit1 <- cfa(
  HS.model,
  data = HolzingerSwineford1939,
  missing = "ML",
  estimator = "MLR",
  bounds = "pos.var", # forces all variances of both observed and latent variables to be strictly nonnegative
  rstarts = 10, # random starts
  verbose = TRUE) # print all output

fit2 <- cfa(
  HS.model,
  data = HolzingerSwineford1939,
  missing = "ML",
  estimator = "MLR",
  bounds = "standard", # uses bounds for observed and latent variances, and for factor loadings
  rstarts = 10, # random starts
  verbose = TRUE) # print all output
```

# Power Analysis {#sec-powerAnalysis}

Overview of power analysis in SEM: <https://web.pdx.edu/~newsomj/semclass/ho_power.pdf> (archived at <https://perma.cc/NBB5-C7FL>)

Example of Monte Carlo simulation  using `lavaan`: <https://isaactpetersen.github.io/Principles-Psychological-Assessment/structural-equation-modeling.html#sec-monteCarloPowerAnalysis>

Example of Monte Carlo simulation using `Mplus`: <https://devpsylab.github.io/DataAnalysis/mplus.html#sec-monteCarlo>

Toolboxes and webapps:

- [`simsem` `R` package](https://cran.r-project.org/web/packages/simsem/index.html) (see vignettes [here](https://github.com/simsem/simsem/wiki/Vignette))
- [`semPower` `R` package](https://moshagen.github.io/semPower/) (based on [Moshagen & Bader, 2024](https://doi.org/10.3758/s13428-023-02254-7)), with basic functionality in a Shiny webapp: <https://sempower.shinyapps.io/sempower/>
- <https://yilinandrewang.shinyapps.io/pwrSEM/> (based on [Wang & Rhemtulla, 2001](https://doi.org/10.1177/2515245920918253))
- <https://schoemanna.shinyapps.io/mc_power_med/>
- <https://sjak.shinyapps.io/power4SEM/> (based on [Jak et al., 2020](https://doi.org/10.3758/s13428-020-01479-0))
- Basic webapps: <https://webpower.psychstat.org/wiki/models/index> and <https://www.danielsoper.com/statcalc/calculator.aspx?id=89>

# Path Diagrams {#sec-pathDiagrams}

For a list of tools to create path diagrams, see [here](figures.qmd#sec-pathDiagrams).

# Partial Regression Plot {#sec-partialRegressionPlot}

```{r}
# simulate extra observed covariates
N <- nrow(PoliticalDemocracy)
PoliticalDemocracy$z1 <- rnorm(N)
PoliticalDemocracy$z2 <- 0.3*PoliticalDemocracy$x1 + 0.3*PoliticalDemocracy$x2 + 0.3*PoliticalDemocracy$x3 + rnorm(N)

# model syntax
myModel <- ' 
 # latent variables 
   ind60 =~ x1 + x2 + x3 
   dem60 =~ y1 + y2 + y3 + y4 
   dem65 =~ y5 + y6 + y7 + y8 
 # regressions
   dem60 ~ ind60 
   dem65 ~ ind60 + dem60 + z1 + z2
 # residual covariances 
   y1 ~~ y5
   y2 ~~ y4 + y6 
   y3 ~~ y7 
   y4 ~~ y8
   y6 ~~ y8
'

# fit the model
fit <- sem(
  model = myModel,
  data = PoliticalDemocracy)

# extract the factor scores
factorScores <- as.data.frame(lavaan::lavPredict(fit))

# merge the factor scores with the original data
mydata <- cbind(PoliticalDemocracy, factorScores)

# residualize latent Y (dem65) on covariates
ry <- resid(lm(dem65 ~ dem60 + z1 + z2, data = mydata))

# residualize latent X (ind60) on same covariates
rx <- resid(lm(ind60 ~ dem60 + z1 + z2, data = mydata))

# create the partial regression plot
ggplot(
  data = data.frame(rx, ry),
  mapping = aes(rx, ry)) +
  geom_point() +
  geom_smooth(
    method = "lm") +
  labs(
    x = "Residualized ind60",
    y = "Residualized dem65",
    title = "Partial Regression Plot from SEM"
  )
```

# Session Info

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

sessionInfo()
```