---
title: "Tabulating results from structural equation models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tabulating results from structural equation models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
run_mplus <- FALSE
```

```{r setup, message=FALSE}
library(tidySEM)
library(lavaan)
library(MplusAutomation)
```

`tidySEM` tabulates the results of different types of models in the same uniform way. This facilitates parsing the output into Tables and Figures for publication. The function to tabulate output is `table_results()`. The function to tabulate fit indices is `table_fit()`.

Let's use a classic `lavaan` tutorial example for a multiple group model, using the `HolzingerSwineford1939` data. The `tidySEM` package has a function `measurement()` to generate measurement models automatically. It guesses which latent variable an observed variable belongs to by splitting the names (by default, at the last `_` symbol), so it helps to rename the variables:

```{r, echo = TRUE, eval = FALSE}
df <- HolzingerSwineford1939
names(df)[7:15] <- paste0(rep(c("vis", "tex", "spe"), each = 3), "_", rep(1:3, 3))
df |>
  subset(select = c("school", "vis_1", "vis_2", "vis_3", "tex_1", "tex_2", "tex_3", "spe_1", 
"spe_2", "spe_3")) -> df
```
```{r, echo = FALSE, eval = TRUE}
df <- HolzingerSwineford1939
names(df)[7:15] <- paste0(rep(c("vis", "tex", "spe"), each = 3), "_", rep(1:3, 3))
subset(df, select = c("school", "vis_1", "vis_2", "vis_3", "tex_1", "tex_2", "tex_3", "spe_1", 
"spe_2", "spe_3")) -> df
```

Now, let's construct the model.

```{r, echo = TRUE, eval = FALSE}
df |>
  tidy_sem() |>
  measurement() -> model
```
```{r, eval = TRUE, echo = FALSE}
model <- measurement(tidy_sem(df))
```

## Output from lavaan

Now, let's run the model in `lavaan`. You can either use `lavaan` to run it, 

```{r, eval = FALSE, echo = TRUE}
model |>
  estimate_lavaan() -> fit_lav
```
```{r, echo = FALSE, eval = TRUE}
estimate_lavaan(model) -> fit_lav
```

The results can be tabulated using `table_results()`:

```{r}
table_results(fit_lav)
```

```{r}
table_fit(fit_lav)
```

## Output from OpenMx

Now, we'll reproduce the same analysis in 'OpenMx'. First, we run the model:

```{r, echo = TRUE, eval = FALSE}
model |>
  estimate_mx() -> fit_mx
table_results(fit_mx)
table_fit(fit_mx)
```

```{r, eval = TRUE, echo = FALSE}
estimate_mx(model) -> fit_mx
table_results(fit_mx)
table_fit(fit_mx)
```

## Output from Mplus

Now, we'll reproduce the same analysis in 'Mplus'. To illustrate the fact that `tidySEM` is compatible with existing solutions, we will specify the syntax for this example manually, using the package `MplusAutomation`. This code will only work on your machine if you have Mplus installed and R can find it. First, we run the model:

```{r, eval = FALSE, echo = TRUE}
fit_mplus <- mplusModeler(mplusObject(VARIABLE = "grouping IS school (1 = GW 2 = Pas);",
                                MODEL = c("visual BY vis_1 vis_2 vis_3;",
                                          "textual BY tex_1 tex_2 tex_3;",
                                          "speed BY spe_1 spe_2 spe_3;"),
                                usevariables = names(df),
                                rdata = df),
                    modelout = "example.inp",
                    run = 1L)
table_results(fit_mplus)
```
```{r, eval = run_mplus, echo = FALSE}
# fit <- mplusModeler(mplusObject(VARIABLE = "grouping IS school (1 = GW 2 = Pas);",
#                                 MODEL = c("visual BY x1 x2 x3;",
#                                           "textual BY x4 x5 x6;",
#                                           "speed BY x7 x8 x9;"),
#                                 usevariables = c(paste0("x", 1:9), "school"),
#                                 rdata = HolzingerSwineford1939),
#                     modelout = "example.inp",
#                     run = 1L)
# file.remove(list.files(pattern = "^example.+(inp|out|dat)$"))
#dput(fit$results$parameters)
#dput(fit, file = "mplusfit.R")
```
```{r eval = TRUE, echo = FALSE}
# Read the results
fit_mplus <- source("mplusfit.R")
fit_mplus <- fit_mplus$value
```
```{r}
table_results(fit_mplus)
table_fit(fit_mplus)
```