This vignette illustrated `tidySEM`

’s ability to perform
latent class growth analysis, or growth mixture modeling. The simulated
data used for this example are inspired by work in progress by Plas and
colleagues, on heterogeneity in depression trajectories among Dutch
military personnel who were deployed to Afghanistan. The original data
were collected as part of the *Prospection in Stress-related Military
Research (PRISMO)* study, which examined of psychological problems
after deployment in more than 1,000 Dutch military personnel from
2005-2019.

First, we load all required packages:

```
library(tidySEM)
library(ggplot2)
library(MASS)
```

We first examined the descriptive statistics for the sum score scales:

```
# Get descriptives
<- plas_depression
df <- descriptives(df)
desc <- desc[, c("name", "mean", "median", "sd", "min", "max",
desc "skew_2se", "kurt_2se")]
::kable(desc, caption = "Item descriptives") knitr
```

Note that all variables were extremely right-skewed due to censoring at the lower end of the scale.

We can examine these distributions visually as well:

```
<- reshape(df, direction = "long", varying = names(df))
df_plot ggplot(df_plot, aes(x = scl)) + geom_density() + facet_wrap(~time) +
theme_bw()
```

As this type of skew can result in convergence problems in LCGA, we compared several transformations to reduce skew: The square and cube root, log, inverse, and Box-Cox transformations.

```
<- df_plot
df_scores # Store original range of SCL
<- range(df_scores$scl)
rng_scl # Log-transform
$log <- scales::rescale(log(df_scores$scl), to = c(0,
df_scores1))
# Square root transform
$sqrt <- scales::rescale(sqrt(df_scores$scl), to = c(0,
df_scores1))
# Cube root transform
$qrt <- scales::rescale(df_scores$scl^0.33, to = c(0,
df_scores1))
# Reciprocal transform
$reciprocal <- scales::rescale(1/df_scores$scl, to = c(0,
df_scores1))
# Define function for Box-Cox transformation
<- function(x, lambda) {
bc ^lambda) - 1)/lambda)
(((x
}# Inverse Box-Cox transformation
<- function(x, lambda) {
invbc * lambda) + 1)^(1/lambda)
((x
}# Box-Cox transform
<- MASS::boxcox(lm(df_scores$scl ~ 1), plotit = FALSE)
b <- b$x[which.max(b$y)]
lambda $boxcox <- bc(df_scores$scl, lambda)
df_scores# Store range of Box-Cox transformed data
<- range(df_scores$boxcox)
rng_bc $boxcox <- scales::rescale(df_scores$boxcox, to = c(0,
df_scores1))
# Rescale SCL
$scl <- scales::rescale(df_scores$scl, to = c(0, 1)) df_scores
```

We can plot these transformations:

```
# Make plot data
<- do.call(rbind, lapply(c("scl", "log", "sqrt", "qrt",
df_plot "boxcox"), function(n) {
data.frame(df_scores[c("time", "id")], Value = df_scores[[n]],
Transformation = n)
}))# Plot
ggplot(df_plot, aes(x = Value, colour = Transformation)) + geom_density() +
facet_wrap(~time) + scale_y_sqrt() + xlab("scl (rescaled to 0-1)") +
theme_bw()
```

Evidently, the Box-Cox transformation reduced skew the most. Consequently, we proceeded with the Box-Cox transformed scores for analysis.

```
<- df_scores[, c("id", "time", "boxcox")]
dat <- reshape(dat, direction = "wide", v.names = "boxcox", timevar = "time",
dat idvar = "id")
names(dat) <- gsub("boxcox.", "scl", names(dat))
```

Next, we estimated a latent class growth model for SCL. The model
included an overall intercept, centered at T1, `i`

. To model
the potential effect of deployment on depresion, we also included a
dummy variable that was zero before deployment, and 1 after deployment,
`step`

. Finally, to model potential change (or recovery) in
depression post-deployment, we included a linear slope from T2-T6,
`s`

. All variances of growth parameters were fixed to zero
due to the sparse nature of the data. In this vignette, we do not
consider more than 5 classes, because the analyses are computationally
very intensive and the data were simulated from a 3-class model.

It is important to highlight that in LCGA, the subgroups will be limited by the specify growth structure, meaning that LCA will identify distinctive growth patterns within the intercept, step, and slope growth. For example, if there is a subgroup that follows a quadratic growth pattern this models will not be able to identify it.

**NOTE: The time scales in this model are not correct; it
currently assumes that all measurements are equidistant. Feel free to
experiment with adjusting this.**

```
set.seed(27796)
"id"]] <- NULL
dat[[<- mx_growth_mixture(model = "
res_step i =~ 1*scl1 + 1*scl2 + 1*scl3 +1*scl4 +1*scl5 +1*scl6
step =~ 0*scl1 + 1*scl2 + 1*scl3 +1*scl4 +1*scl5 +1*scl6
s =~ 0*scl1 + 0*scl2 + 1*scl3 +2*scl4 +3*scl5 +4*scl6
scl1 ~~ vscl1*scl1
scl2 ~~ vscl2*scl2
scl3 ~~ vscl3*scl3
scl4 ~~ vscl4*scl4
scl5 ~~ vscl5*scl5
scl6 ~~ vscl6*scl6
i ~~ 0*i
step ~~ 0*step
s ~~ 0*s
i ~~ 0*s
i ~~ 0*step
s ~~ 0*step",
classes = 1:5, data = dat)
# Additional iterations because of convergence problems for
# model 1:
1]] <- mxTryHardWideSearch(res_step[[1]], extraTries = 50) res_step[[
```

Note that the first model showed convergence problems, throwing the
error: *The model does not satisfy the first-order optimality
conditions to the required accuracy, and no improved point for the merit
function could be found during the final linesearch.* To address
this problem, we performed additional iterations to

find a better solution, using
`OpenMx::mxTryHardWideSearch()`

. This also illustrates that
`tidySEM`

mixture models inherit from `OpenMx`

’s
`MxModel`

, and thus, different `OpenMx`

functions
can be used to act on models specified via `tidySEM`

.

The fifth model also evidenced convergence problems, but this (as we will see) is because the solution is overfitted.

To determine the correct number of classes, we considered the following criteria:

- We do not consider classes with, on average, fewer than 5 participants per parameter in a class due to potential local underidentification
- Lower values for information criteria (AIC, BIC, saBIC) indicate better fit
- Significant Lo-Mendell-Rubin LRT test indicates better fit for \(k\) vs \(k-1\) classes
- We do not consider solutions with entropy < .90 because poor class separability compromises interpretability of the results
- We do not consider solutions with minimum posterior classification probability < .90 because poor class separability compromises interpretability of the results

```
# Get fit table fit
<- table_fit(res_step)
tab_fit # Select columns
c("Name", "Classes", "LL", "Parameters", "BIC", "Entropy",
tab_fit[, "prob_min", "n_min", "warning", "lmr_p")]
```

According to the Table, increasing the number of classes keeps increasing model fit according to all ICs except the BIC, which increased after 3 classes.

The first two LMR tests are significant, indicating that a 2- and 3-class solution were a significant improvement over a 1- and 2-class solution, respectively. However, solutions with >3 classes had entropy and minimum posterior classification probability below the pre-specified thresholds. Models with >3 solutions also had fewer than five observations per parameter. This suggests that the preferred model should be selected from 1-3 classes.

A scree plot indicates that the largest decrease in ICs occurs from 1-2 classes, and the inflection point for all ICs is at 3 classes. Moreover, the BIC increased after 3 classes. A three-class solution thus appears to be the most parsimonious solution with good fit.

`plot(tab_fit, statistics = c("AIC", "BIC", "saBIC"))`

Based on the aforementioned criteria, we selected a 3-class model for
further analyses. First, to prevent label switching, we re-order these
classes by the value of the intercept `i`

. Then, we report
the estimated parameters.

```
<- mx_switch_labels(res_step[[3]], param = "M[1,7]",
res_final decreasing = FALSE)
<- table_results(res_final, columns = NULL)
tab_res # Select rows and columns
<- tab_res[tab_res$Category %in% c("Means", "Variances"),
tab_res c("Category", "lhs", "est", "se", "pval", "confint", "name")]
tab_res
```

As evident from these results, Class 1 started at a relatively lower level of depressive symptoms, experienced a decrease after deployment, followed by increase over time. Class 2 started at a moderate level of depressive symptoms, experienced an increase after deployment, followed by significant increase over time from T2-T6. Class 3 started at a relatively higher level, experienced an increase after deployment, followed by stability.

To test whether parameters are significantly different between
classes, we can use Wald tests. Wald tests can be specified for all
parameters in the model, using the hypothesis syntax from the
`bain`

package for informative hypothesis testing.

To identify the names of parameters in the model, we can use the
`name`

column of the results table above. Alternatively, to
see all parameters in the model, run:

`names(coef(res_final))`

Next, specify equality constrained hypotheses. For example, a hypothesis that states that the mean intercept is equal across groups is specified as follows:

`"class1.M[1,7] = class2.M[1,7] & class1.M[1,7] = class3.M[1,7]`

It is also possible to consider comparisons between two classes at a time. When conducting many significance tests, consider correcting for multiple comparisons however.

```
<- wald_test(res_final, "
wald_tests class1.M[1,7] = class2.M[1,7]&
class1.M[1,7] = class3.M[1,7];
class1.M[1,8] = class2.M[1,8]&
class1.M[1,8] = class3.M[1,8];
class1.M[1,9] = class2.M[1,9]&
class1.M[1,9] = class3.M[1,9]")
# Rename the hypothesis
$Hypothesis <- c("Mean i", "Mean step", "Mean slope")
wald_tests::kable(wald_tests, digits = 2, caption = "Wald tests") knitr
```

All Wald tests are significant, indicating that there are significant differences between the intercepts, step function, and slopes of the three classes.

Finally, we can plot the growth trajectories. This can help interpret the results better, as well as the residual heterogeneity around class trajectories.

```
<- plot_growth(res_step[[3]], rawdata = TRUE, alpha_range = c(0,
p 0.05))
# Add Y-axis breaks in original scale
<- seq(0, 1, length.out = 5)
brks <- round(invbc(scales::rescale(brks, from = c(0, 1), to = rng_bc),
labs
lambda))<- p + scale_y_continuous(breaks = seq(0, 1, length.out = 5),
p labels = labs) + ylab("SCL (rescaled from Box-Cox)")
p
```

Note that the observed individual trajectories show very high variability within classes.