---
title: "Full CFA/SEM workflow"
author: "Rémi Thériault"
date: "August 29, 2022"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Full CFA/SEM workflow}
  %\VignetteEngine{knitr::rmarkdown}
---

## CFA example

```{r cfaplot1, message = FALSE}
# Load library
library(lavaan)
library(lavaanExtra)
library(tibble)

# Define latent variables
latent <- list(visual = paste0("x", 1:3),
               textual = paste0("x", 4:6),
               speed = paste0("x", 7:9))

# Write the model, and check it
cfa.model <- write_lavaan(latent = latent)
cat(cfa.model)

# Fit the model fit and plot with `lavaanExtra::cfa_fit_plot`
# to get the factor loadings visually (optionally as PDF)
fit.cfa <- cfa_fit_plot(cfa.model, HolzingerSwineford1939)
```

<img src="https://lavaanExtra.remi-theriault.com/articles/images/fullworkflow/cfaplot.png" width="30%" />

```{r cfa2}
# Get fit indices
nice_fit(fit.cfa)

# We can get it prettier with the `rempsyc::nice_table` integration
nice_fit(fit.cfa, nice_table = TRUE)

```

But let's say you had a bad fit and wanted to remove the three items with the lowest loadings, you can do so without have to respecify the model, only what items you wish to remove:

```{r cfaplot4}
# Fit the model fit and plot with `lavaanExtra::cfa_fit_plot`
# to get the factor loadings visually (as PDF)
fit.cfa2 <- cfa_fit_plot(cfa.model, HolzingerSwineford1939,
                         remove.items = paste0("x", c(2:3, 7)))
```

<img src="https://lavaanExtra.remi-theriault.com/articles/images/fullworkflow/cfaplot2.png" width="30%" />

Let's compare the fit to see if it's better now:

```{r cfaplot5}
nice_fit(lst(fit.cfa, fit.cfa2), nice_table = TRUE)
```

It is! If you like this table, you may also wish to save it to Word. Also easy:

```{r cfaplot6}
# Save fit table as an object
fit_table <- nice_fit(lst(fit.cfa, fit.cfa2), nice_table = TRUE)

# Save fit table to Word!
flextable::save_as_docx(fit_table, path = "fit_table.docx")
```

```{r, echo=FALSE}
unlink("fit_table.docx")
```

## SEM example

Here is a structural equation model example. We start with a path analysis first.

### Saturated model

One might decide to look at the saturated `lavaan` model first.

```{r saturated}
# Calculate scale averages
data <- HolzingerSwineford1939
data$visual <- rowMeans(data[paste0("x", 1:3)])
data$textual <- rowMeans(data[paste0("x", 4:6)])
data$speed <- rowMeans(data[paste0("x", 7:9)])

# Define our variables
M <- "visual"
IV <- c("ageyr", "grade")
DV <- c("speed", "textual")

# Define our lavaan lists
mediation <- list(speed = M, textual = M, visual = IV)
regression <- list(speed = IV, textual = IV)
covariance <- list(speed = "textual", ageyr = "grade")

# Write the model, and check it
model.saturated <- write_lavaan(mediation, regression, covariance)
cat(model.saturated)

```

This looks good so far, but we might also want to check our indirect effects (mediations). For this, we have to obtain the path names by setting `label = TRUE`. This will allow us to define our indirect effects and feed them back to `write_lavaan`.

```{r saturated2}
# We can run the model again. 
# However, we set `label = TRUE` to get the path names
model.saturated <- write_lavaan(mediation, regression, covariance, label = TRUE)
cat(model.saturated)

```

Here, if we check the mediation section of the model, we see that it has been "augmented" with the path names. Those are `visual_speed`, `visual_textual`, `ageyr_visual`, and `grade_visual`. The logic for the determination of the path names is predictable: it is always the predictor variable, on the left, followed by the predicted variable, on the right. So if we were to test all possible indirect effects, we would define our `indirect` object as such:

```{r saturated3}
# Define indirect object
indirect <- list(ageyr_visual_speed = c("ageyr_visual", "visual_speed"),
                 ageyr_visual_textual = c("ageyr_visual", "visual_textual"),
                 grade_visual_speed = c("grade_visual", "visual_speed"),
                 grade_visual_textual = c("grade_visual", "visual_textual"))

# Write the model, and check it
model.saturated <- write_lavaan(mediation, regression, covariance, 
                                indirect, label = TRUE)
cat(model.saturated)

```

If preferred (e.g., when dealing with long variable names), one can choose to use letters for the predictor variables. Note however that this tends to be somewhat more confusing and ambiguous.

```{r letters}
# Write the model, and check it
model.saturated <- write_lavaan(mediation, regression, covariance, 
                                label = TRUE, use.letters = TRUE)
cat(model.saturated)
```

In this case, the path names are `a_speed`, `a_textual`, `a_visual`, and `b_visual`. So we would define our `indirect` object as such:

```{r letters2}
# Define indirect object
indirect <- list(ageyr_visual_speed = c("a_visual", "a_speed"),
                 ageyr_visual_textual = c("a_visual", "a_textual"),
                 grade_visual_speed = c("b_visual", "a_speed"),
                 grade_visual_textual = c("b_visual", "a_textual"))

# Write the model, and check it
model.saturated <- write_lavaan(mediation, regression, covariance, 
                                indirect, label = TRUE, use.letters = TRUE)
cat(model.saturated)

```

There is also an experimental feature that attempts to produce the indirect effects automatically. This feature requires specifying your independent, dependent, and mediator variables as "IV", "M", and "DV", respectively, in the `indirect` object. In our case, we have already defined those earlier, so we can just feed the proper objects.

```{r saturated4}
# Define indirect object
indirect <- list(IV = IV, M = M, DV = DV)

# Write the model, and check it
model.saturated <- write_lavaan(mediation, regression, covariance, 
                                indirect, label = TRUE)
cat(model.saturated)

```

We are now satisfied with our model, so we can finally fit it!

```{r saturated5}
# Fit the model with `lavaan`
fit.saturated <- sem(model.saturated, data = data)

# Get regression parameters only 
# And make it pretty with the `rempsyc::nice_table` integration
lavaan_reg(fit.saturated, nice_table = TRUE, highlight = TRUE)

```

So `speed` as predicted by `ageyr` isn't significant. We could remove that path from the model it if we are trying to make a more parsimonious model. Let's make the non-saturated path analysis model next.

### Path analysis model

Because we use `lavaanExtra`, we don't have to redefine the entire model: simply what we want to update. In this case, the regressions and the indirect effects.

```{r path}
regression <- list(speed = "grade", textual = IV)

# We can run the model again, setting `label = TRUE` to get the path names
model.path <- write_lavaan(mediation, regression, covariance, label = TRUE)
cat(model.path)
# We check that we have removed "ageyr" correctly from "speed" in the 
# regression section. OK.

# Define just our indirect effects of interest
indirect <- list(age_visual_speed = c("ageyr_visual", "visual_speed"),
                 grade_visual_textual = c("grade_visual", "visual_textual"))

# We run the model again, with the indirect effects
model.path <- write_lavaan(mediation, regression, covariance, 
                           indirect, label = TRUE)
cat(model.path)

# Fit the model with `lavaan`
fit.path <- sem(model.path, data = data)

# Get regression parameters only
lavaan_reg(fit.path)

# We can get it prettier with the `rempsyc::nice_table` integration
lavaan_reg(fit.path, nice_table = TRUE, highlight = TRUE)

# We only kept significant regressions. Good (for this demo).

# Get covariance indices
lavaan_cov(fit.path)

# We can get it prettier with the `rempsyc::nice_table` integration
lavaan_cov(fit.path, nice_table = TRUE)

# Get nice fit indices with the `rempsyc::nice_table` integration
nice_fit(lst(fit.cfa, fit.saturated, fit.path), nice_table = TRUE)

# Let's get the indirect effects only
lavaan_ind(fit.path)

# We can get it prettier with the `rempsyc::nice_table` integration
lavaan_ind(fit.path, nice_table = TRUE)

# Get modification indices only
modindices(fit.path, sort = TRUE, maximum.number = 5)
```

For reference, this is our model, visually speaking

![](https://lavaanExtra.remi-theriault.com/articles/images/fullworkflow/holzinger_model.png)

We could also attempt to draw it with `lavaanExtra::nice_tidySEM`, a convenience wrapper around the amazing `tidySEM` package.

```{r nice_tidySEM}
labels <- list(ageyr = "Age (year)",
               grade = "Grade",
               visual = "Visual",
               speed = "Speed",
               textual = "Textual")
layout <- list(IV = IV, M = M, DV = DV)

nice_tidySEM(fit.path, layout = layout, label = labels,
             hide_nonsig_edges = TRUE, label_location = .60)

```


### Latent model

Finally, perhaps we change our mind and decide to run a full SEM instead, with latent variables. Fear not: we don't have to redo everything again. We can simply define our latent variables and proceed. In this example, we have *already* defined our latent variable for our CFA earlier, so we don't even need to write that again!

```{r latent}
model.latent <- write_lavaan(mediation, regression, covariance, 
                             indirect, latent, label = TRUE)
cat(model.latent)

# Run model
fit.latent <- sem(model.latent, data = HolzingerSwineford1939)

# Get nice fit indices with the `rempsyc::nice_table` integration
nice_fit(lst(fit.cfa, fit.saturated, fit.path, fit.latent), nice_table = TRUE)

```
