sem.Rmd

---
title: "Structural Equation Modeling"
output:
  html_document:
    code_folding: show
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  error = TRUE,
  comment = "",
  class.source = "fold-show")

options(scipen = 999)
```

# Preamble

## Install Libraries

```{r, class.source = "fold-hide"}
#install.packages("remotes")
#remotes::install_github("DevPsyLab/petersenlab")
```

## Load Libraries

```{r, message = FALSE, warning = FALSE, class.source = "fold-hide"}
library("lavaan")
library("semTools")
library("semPlot")
library("lcsm")
library("MBESS")
library("tidyverse")
```

# Simulate Data

```{r, class.source = "fold-hide"}
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

# Overview

https://isaactpetersen.github.io/Principles-Psychological-Assessment/sem.html

# Analysis examples

https://isaactpetersen.github.io/Principles-Psychological-Assessment/sem.html#semModelExample-sem

# Plot Observed Growth Curve

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() +
  scale_x_continuous(
    breaks = 1:4,
    name = "Timepoint") +
  scale_y_continuous(
    name = "Score")
```

# Latent Growth Curve Model {#lgcm}

## 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
  t2 ~ 0
  t3 ~ 0
  t4 ~ 0
  
  # 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}
semPaths(
  lgcm1_fit,
  what = "Std.all",
  layout = "tree2",
  edge.label.cex = 1.5)
```

## Plot Trajectories

### Protoypical Growth Curve

Calculated from intercept and slope parameters:

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

ggplot() +
  xlab("Timepoint") +
  ylab("Score") +
  scale_x_continuous(
    limits = c(0, 3),
    labels = 1:4) +
  scale_y_continuous(
    limits = c(0, 5)) +
  geom_abline(
    mapping = aes(
      slope = lgcm1_slope,
      intercept = lgcm1_intercept))
```

Calculated manually:

```{r}
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)) +
  xlab("Timepoint") +
  ylab("Score") +
  scale_y_continuous(
    limits = c(0, 5)) +
  geom_line()
```

### Individuals' Growth Curves

Calculated from intercept and slope parameters:

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

ggplot(
  data = newData) +
  xlab("Timepoint") +
  ylab("Score") +
  scale_x_continuous(
    limits = c(0, 3),
    labels = 1:4) +
  scale_y_continuous(
    limits = c(-10, 20)) +
  geom_abline(
    mapping = aes(
      slope = slope,
      intercept = intercept))
```

Calculated manually:

```{r}
newData$t1 <- newData$intercept
newData$t4 <- newData$intercept + (timepoints - 1)*newData$slope

newData2 <- pivot_longer(
  newData,
  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))) +
  xlab("Timepoint") +
  ylab("Score") +
  scale_y_continuous(
    limits = c(-10, 20)) +
  geom_line()
```

### Individuals' Trajectories Overlaid with Prototypical Trajectory

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

ggplot(
  data = newData) +
  xlab("Timepoint") +
  ylab("Score") +
  scale_x_continuous(
    limits = c(0, 3),
    labels = 1:4) +
  scale_y_continuous(
    limits = c(-10, 20)) +
  geom_abline(
    mapping = aes(
      slope = slope,
      intercept = intercept)) +
  geom_abline(
    mapping = aes(
      slope = lgcm1_slope,
      intercept = lgcm1_intercept),
    color = "blue",
    linewidth = 2)
```

# Latent Change Score Model {#lcsm}

## Model Syntax

```{r}
bivariateLCSM_syntax <- specify_bi_lcsm(
  timepoints = 10,
  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 = FALSE), # 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)
```

## Plot Trajectories

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

# Random Intercept Cross-Lagged Panel Model {#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)
```

# Mediation {#mediation}

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

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

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

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

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

mediationFit_indirect
```

#### Proportion Mediated (*P*<sub>*M*</sub>) {#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 == "total") %>% 
  select(std.all) %>% 
  as.numeric

mediationFit_pm <- mediationFit_indirect / 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 renders the estimate of proportion mediated not a meaningful estimate of effect size (which explains why it the estimate exceeds 1.0; Fairchild & McDaniel, 2017).

#### Proportion of Variance in Y That is Explained by the Indirect Effect (*R*<sup>2</sup><sub>mediated</sub>) {#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$) {#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$) {#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

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

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

## Interaction Plot {#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 {#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 {#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(b1, b3)} \\
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(b1, b3) \\
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()
```

# Bounded Estimation with Random Starts {#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, message = FALSE, warning = FALSE, results = FALSE}
fit <- 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
```

# Power Analysis {#powerAnalysis}

https://isaactpetersen.github.io/Principles-Psychological-Assessment/sem.html#monteCarloPowerAnalysis

- https://yilinandrewang.shinyapps.io/pwrSEM/
- https://schoemanna.shinyapps.io/mc_power_med/
- https://sjak.shinyapps.io/power4SEM/
- https://sempower.shinyapps.io/sempower/
- https://webpower.psychstat.org/wiki/models/index

# Path Diagrams {#pathDiagrams}

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

# Session Info

```{r, class.source = "fold-hide"}
sessionInfo()
```