Overview

This guide provides examples of how to use the functions in the MplusAutomation package. The package is designed to automate three major aspects of latent variable modeling in Mplus:

1. creating related groups of models
2. running batches of models
3. extracting and comparing model parameters and fit statistics

The guide tries to make few assumptions about the user’s familiarity with the R environment, as many Mplus user’s may be unfamiliar with R. Note that the MplusAutomation package was written to be compatible with Mplus version 5.21 or later. Although many functions are likely to work with earlier versions, there may be incompatibilities.

Note: some examples herein reference examples from the Mplus User’s Guide. The input and output files for these examples are assumed to reside in the directory:

C:\Program Files\Mplus\Mplus Examples\User's Guide Examples

If you have installed Mplus to a different location, please use the correct directory in the code below.

Installing and loading the package

Note: You only need to install the MplusAutomation package once on a given machine. Once you follow these instructions, you can just use library(MplusAutomation) to load the package in future sessions.

The MplusAutomation package is compatible with R and Mplus on Windows, Mac OS X, and Linux platforms.

The package was initially built using R 2.15.1 and has been tested with modern versions up through R 4.0.2. To obtain R follow this link.

After installing R and launching a new session, type the following at the command line:

install.packages("MplusAutomation")

If you’re running on Windows, you may need to run R as an administrator (right-click, Run as administrator).

To load the MplusAutomation package so that its functions are available to use, type:

library(MplusAutomation)

## Version:  1.1.1
## We work hard to write this free software. Please help us get credit by citing: 
## 
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 25, 621-638. doi: 10.1080/10705511.2017.1402334.
## 
## -- see citation("MplusAutomation").

Note: Mac users using versions of R prior to 4.0.0 may experiences crashes if they have not setup Tcl/Tk properly within R. The default distribution of Tcl/Tk shipped with Mac OS X is not compatible with R. Instead, users should download and install the Tcl/Tk distribution here.

sessionInfo()

## R version 4.3.2 (2023-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.2.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] MplusAutomation_1.1.1
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4       jsonlite_1.8.8     dplyr_1.1.4        compiler_4.3.2    
##  [5] tidyselect_1.2.0   Rcpp_1.0.11        parallel_4.3.2     jquerylib_0.1.4   
##  [9] scales_1.3.0       boot_1.3-28.1      yaml_2.3.8         fastmap_1.1.1     
## [13] lattice_0.22-5     coda_0.19-4        ggplot2_3.4.4      R6_2.5.1          
## [17] plyr_1.8.9         generics_0.1.3     knitr_1.45         backports_1.4.1   
## [21] checkmate_2.3.1    tibble_3.2.1       pander_0.6.5       gsubfn_0.7        
## [25] proto_1.0.0        munsell_0.5.0      bslib_0.6.1        pillar_1.9.0      
## [29] rlang_1.1.2        utf8_1.2.4         cachem_1.0.8       xfun_0.41         
## [33] sass_0.4.8         cli_3.6.2          magrittr_2.0.3     fastDummies_1.7.3 
## [37] digest_0.6.33      grid_4.3.2         xtable_1.8-4       texreg_1.39.3     
## [41] lifecycle_1.0.4    vctrs_0.6.5        data.table_1.14.10 evaluate_0.23     
## [45] glue_1.6.2         fansi_1.0.6        colorspace_2.1-0   rmarkdown_2.25    
## [49] httr_1.4.7         tools_4.3.2        pkgconfig_2.0.3    htmltools_0.5.7

update.packages(ask=FALSE, checkBuilt=TRUE)

Exporting data from R to Mplus

The helper function, prepareMplusData(), eases the task of transferring data from R to Mplus. Mplus requires that files not have a header row and that the variable names be specified within the Mplus input syntax. The prepareMplusData() function converts an R data.frame object (the typical way to represent two-dimensional data in R) to a tab-delimited file and it prints the corresponding Mplus syntax to the console. This syntax can then be pasted into the head of a new Mplus input file. Here are two basic examples of using the command:

prepareMplusData(
  my_data,
  filename = "C:/Data_Analysis/Prepare Mplus.dat",
  keepCols=c("id", "item1", "item3", "item6"))

prepareMplusData(
  my_other_data,
  filename = "C:/Data_Analysis/Prepare Dropped Mplus.dat",
  dropCols=c("baditem1", "baditem2", "baditem7"))

The first call requests that only the variables specified by keepCols be included in the resultant Mplus data file. In like manner, the second call requests that all variables except those specified by dropCols by included in the data file. If neither keepCols nor dropCols is specified, then the entire data.frame is output.

For factors and character variables (categorical), prepareMplusData converts each level into an integer value (first level = 1, second level = 2, etc.) and prints out information about how the levels were converted so that you can use them properly inside the Mplus syntax. Alternatively, if you use the dummyCode argument for prepareMplusData, you will get a set of dummy codes that can be used in the Mplus syntax.

Here, in the mtcars dataset, we convert the number of gears (in the car) to a factor so that prepareMplusData will know to convert its unique values (3, 4, 5) into an integer-valued variable (1, 2, 3) that can be used in Mplus syntax. For the number of cylinders (cyl) and transmission type (am), we ask prepareMplusData to generate dummy codes for each level of these factors. Note that one dummy code is returned for each level, so no reference category is implied. Thus, it is up to you to drop one of the dummy codes in the Mplus syntax to denote the reference category and avoid multicollinearity.

data(mtcars)
mtcars$gear <- factor(mtcars$gear)
prepareMplusData(mtcars, "mtcars.dat", dummyCode = c("cyl", "am"))

## Converting the following variables to dummy codes: cyl, am

## 
## -----
## Factor: gear 
## Conversion:
##  level number
##      3      1
##      4      2
##      5      3
## -----

## The file(s)
##  'mtcars.dat' 
## currently exist(s) and will be overwritten

## TITLE: Your title goes here
## DATA: FILE = "mtcars.dat";
## VARIABLE: 
## NAMES = mpg disp hp drat wt qsec vs gear carb cyl_4 cyl_6 cyl_8 am_0 am_1; 
## MISSING=.;

Running batches of Mplus models

A major purpose of the MplusAutomation package is to allow for easy runs of batches/groups of Mplus models. Oftentimes, one wants to compare a group of related models, such as testing for different forms of measurement invariance. Depending on the complexity of the models, Mplus can take several minutes to many hours to run each model. The routine is designed to run a group of related models located within a directory (or nested subdirectories).

runModels("C:/Program Files/Mplus/Mplus Examples/Addendum Examples")

ComparingLCAvCFA/LCA/1-class LCA.inp
ComparingLCAvCFA/LCA/2-class LCA.inp
ComparingLCAvCFA/LCA/3-class LCA.inp

ComparingLCAvCFA/CFA/1-factor CFA.inp
ComparingLCAvCFA/CFA/2-factor CFA.inp
ComparingLCAvCFA/CFA/3-factor CFA.inp

runModels(
  "C:/Data_Analysis/ComparingLCAvCFA",
  recursive=TRUE)

runModels(
  "C:/Data_Analysis/ComparingLCAvCFA",
  recursive=TRUE,
  logFile="C:/CFALCA-Comparison-Log.txt")|

runModels(
  "C:/Data_Analysis/ComparingLCAvCFA",
  recursive=TRUE,
  logFile=NULL)

runModels(
  "C:/Data_Analysis/ComparingLCAvCFA",
  recursive=TRUE,
  replaceOutfile="never")

runModels(
  "C:/Data_Analysis/ComparingLCAvCFA",
  recursive=TRUE,
  replaceOutfile="modifiedDate")

runModels(
  "C:/Data_Analysis/ComparingLCAvCFA",
  recursive=TRUE,
  showOutput=TRUE)

User-friendly interface to runModels()

A wrapper function, runModels_Interactive(), is included in the MplusAutomation package, which provides a simple dialog box for specifying which models to run. To start the interface, type the following:

runModels_Interactive()

The picture below documents the appearance of this interface:

Figure. Example of using runModels() in an interactive graphical interface.

Although one can provide parameters to the function to set the initial states of the interface, it is rarely necessary to do so, so the syntax above should suffice.

readModels(): Extracting all supported data from Mplus output

The readModels() function provides a single method for extracting all available information (that MplusAutomation currently handles) from one or more output files. This is the preferred method for information from Mplus output files. More specifically, this function provides a wrapper around four other user-accessible functions: extractModelSummaries, extractModelParameters, extractModIndices, and getSavedata_Data. The readModels() function also uses a number of additional functions not accessible to the user to read output such as TECH1, TECH3, TECH4, TECH9, and RESIDUALS. The readModels() function extracts model summaries, parameters, modification indices, and saved data (using the Mplus SAVEDATA syntax), into a single list object (stored as mplus.model). The top-level elements of the list represent distinct output files and are named by the corresponding Mplus output filename. Each mplus.model element of the list is composed of several sub-elements, depending on the information in the output file. Elements of an mplus.model object include:

1. summaries: model summary statistics
2. parameters: model parameter estimates
3. class_counts: class counts and proportions for models that include a categorical latent variable.
4. mod_indices: model modification indices
5. savedata_info: information about SAVEDATA output corresponding to this model
6. savedata: data.frame corresponding to SAVEDATA: FILE = …
7. bparameters: an mcmc.list object containing the draws from the MCMC chains for a Bayesian model that uses the SAVEDATA: BPARAMETERS command
8. residuals: a list containing relevant information from OUTPUT: RESIDUALS.
9. tech1: a list containing parameter specification and starting values from OUTPUT: TECH1.
10. tech3: a list containing parameter covariance and correlation matrices from OUTPUT: TECH3.
11. tech4: a list containing means, covariances, and correlations for latent variables from OUTPUT: TECH4.
12. tech9: a list containing error messages for each replication of a Monte Carlo simulation study from OUTPUT: TECH9.
13. lcCondMeans: conditional latent class means, obtained using auxiliary(e) syntax in latent class models
14. gh5: a list containing data from the gh5 (graphics) file corresponding to this output. (Requires rhdf5 package).

Basic example:

allOutput <- readModels(
  "C:/Data_Files/CFANesting",
  recursive=TRUE)

## assuming there are multiple files in this directory
## just model summaries could retained as a data.frame as follows:

library(plyr)
justSummaries <- do.call("rbind.fill",
  sapply(allOutput,"[", "summaries"))

The rbind.fill function is provided by the plyr package and is used to combine data.frames where the columns do not align perfectly (as happens when Mplus output files differ in the form of their summary statistics). The sapply call is used to extract only the summaries element from each top-level element (i.e., file).

Extracting model summary statistics

Another major purpose of the package is to allow for easy extraction of model summary statistics from one or more models. Such summary statistics include items such as log-likelihood values, root mean squared error of approximation (RMSEA), and Akaike’s Information (AIC).

Please note: The preferred way to obtain summary statistics (and other Mplus output) is the readModels() function, which is described above.

The extractModelSummaries() function is designed to extract model summaries from a group of models located within a directory (or nested within subdirectories). This function returns a data.frame containing one row per model, with columns representing several fit statistics. Note that extractModelSummaries can also extract summaries from a single file by simply passing in a file, not a directory, as the target.

A basic call to the function includes the directory containing output files to be parsed:

mySummaries <- extractModelSummaries(
  "C:/Data_Analysis/ComparingLCAvCFA",
  recursive = TRUE)

Now, the variable mySummaries is a data.frame containing summary statistics about models contained in the ComparingLCAvCFA directory.

As with runModels(), extractModelSummaries() includes a recursive parameter that specifies whether to parse output files located in subdirectories beneath the target directory (defaults to FALSE).

In addition, extractModelSummaries() also includes a parameter, filefilter, that allows the user to parse only files matching certain search criteria. filefilter accepts a Perl-compatible regular expression string. If you’re unfamiliar with regular expressions in Perl, I suggest these two websites:

• http://www.pcre.org/pcre.txt
• https://www.regular-expressions.info/

Note that many regular expression in Perl rely on backslashes (\) for defining character classes, escaping certain characters, and so on. In R, backslashes contained in strings must be doubled (i.e., \\).

Here is an example of filtering only files that match “ex4” followed by any characters Note that the function automatically searches only files with the .out extension, so it isn’t necessary to include .out in the file filter.

Here is a more complex filter that matches filenames that begin with the digits 1, 2, or 3 (for 1-class, 2-class, or 3-class output files) and also contain the text “Threshold”:

summaryStats <- extractModelSummaries(
  "C:/Data_Analysis/Multiclass Models",
  filefilter="[123]{1}-class.*Threshold.*")

Summarizing model fit statistics in tabular form

Once summary statistics for a group of models have been extracted, it is often useful to display them in tabular form to compare fit among models, sorted by a particular criterion (e.g., AIC).

The MplusAutomation package provides three routines for tabulating model summary statistics. At this time, there are three table-generating functions, which are detailed below: showSummaryTable(), HTMLSummaryTable(), LatexSummaryTable().

As their names suggest, these functions can create tables for on-screen display (showSummaryTable()), as an HTML file containing the table (HTMLSummaryTable()), or as a LaTex-formatted table (LatexSummaryTable() ).

showSummaryTable(
  summaryStats,
  keepCols = c("Title", "LL", "AIC", "BIC", "CFI"),
  sortBy = "AIC")

showSummaryTable(
  summaryStats,
  dropCols = c("InputInstructions", "Observations", "Parameters"),
  sortBy = "CFI")

HTMLSummaryTable(
  summaryStats,
  filename = "C:/MyModelSummary.html",
  display = TRUE,
  keepCols = c("Title", "LL", "AIC", "BIC", "AICC"),
  sortBy = "AIC")

myLatexTable <- LatexSummaryTable(
  summaryStats,
  keepCols = c("Title", "BIC", "Parameters"),
  sortBy = "Parameters",
  caption = "Comparing CFA vs. LCA according to number of parameters",
  label="CFALCATab")

<<echo=TRUE, results=tex>>=

myLatexTable <- LatexSummaryTable(
  summaryStats, 
  keepCols = c("Title", "BIC", "Parameters"), 
  sortBy = "Parameters", 
  caption = "Comparing CFA vs. LCA according to number of parameters", 
  label = "CFALCATab")

print(myLatexTable)
@

Extracting model modification indices

The extractModIndices() function extracts the model modification indices from the MODEL MODIFICATION INDICES sections of one or more Mplus output files. It is up to the user to request modification indices using OUTPUT: MODINDICES, and for some models, such indices are not available (as noted in the WARNINGS section of the Mplus output.

Model modification indices are returned as a data.frame of the form:

modV1 operator modV2     MI    EPC Std_EPC StdYX_EPC
BORD1       ON BORD9 12.427 -0.284  -0.284    -0.231
BORD9     WITH BORD1 12.427 -0.053  -0.053    -0.222

These columns follow the conventions used in Mplus where EPC refers to the expected parameter change if the designated relationship is estimated.

Extracting model parameters

The extractModelParameters() function extracts the model parameters from the MODEL RESULTS and STANDARDIZED MODEL RESULTS sections of a given Mplus output file. Examples of such parameters include the parameter estimate, std. err, param/s.e., and two-tailed p-value.

Further, extractmodelParameters() supports extraction of results from many output files, with the results being returned as a list object, one element per output file. When available, unstandardized and standardized (StdYX, StdY, Std) parameters are extracted from each output file into a list object whose elements are data.frame objects. Relatedly, the resultType parameter has been deprecated and will be removed in a future version.

modelResults <- extractModelParameters(
  "C:/Data_Analysis/Mplus Output.out")

unstandardizedResults <- modelResults$unstandardized

#equivalently
standardizedResults <- modelResults[["stdyx.standardized"]]

ComparingLCAvCFA/LCA/1-class LCA.out
ComparingLCAvCFA/LCA/2-class LCA.out
ComparingLCAvCFA/LCA/3-class LCA.out

ComparingLCAvCFA/CFA/1-factor CFA.out
ComparingLCAvCFA/CFA/2-factor CFA.out
ComparingLCAvCFA/CFA/3-factor CFA.out

allModelParameters <- extractModelParameters(
  "C:/Data_Analysis/ComparingLCAvCFA",
  recursive = TRUE)

names(allModelParameters)

## ComparingLCAvCFA.LCA.1.class.LCA.out
## ComparingLCAvCFA.LCA.2.class.LCA.out
## ComparingLCAvCFA.LCA.3.class.LCA.out
## ComparingLCAvCFA.CFA.1.factor.CFA.out
## ComparingLCAvCFA.CFA.2.factor.CFA.out
## ComparingLCAvCFA.CFA.3.factor.CFA.out

TwoFacCFA.STDYX <- allModelParameters$ComparingLCAvCFA.CFA.2.factor.CFA.out$stdyx.standardized

unstandardizedOnly <- sapply(allModelParameters, "[", "unstandardized")

oldNames <- names(allModelParameters)
unstandardizedOnly <- sapply(allModelParameters, "[", "unstandardized")
names(unstandardizedOnly) <- oldNames

#add the filename as a field in the data.frame (so it's uniquely identified when combined)
lapply(names(unstandardizedOnly), function(element) {
    unstandardizedOnly[[element]]$filename <<- element
})

#this will only work if all data.frames have identical columns (i.e., same Mplus output fields)
combinedParameters <- do.call("rbind", unstandardizedOnly)

Capitalizing on the graphics strength of R to visualize results

One of the major strengths of R is its graphics functionality. Numerous functions are provided by base R, such as hist(), plot(), and curve(). Furthermore, impressively flexible and powerful graphics functionality is now provided in R by the grid graphics framework. An useful introduction to graphics in R is provided by Paul Murrell’s book, R Graphics, Second Edition published by CRC Press. In particular, the ggplot2 and lattice graphics package for R provide powerful functions for R graphics.

By contrast, Mplus has very basic graphics functionality that lacks the flexibility and robustness of R.

Now that we have illustrated how to import Mplus model parameter estimates into R, here are just a couple of examples of how useful graphs in can be developed from extractModelParameters() data.

Example: Plotting means and standard errors from a finite mixture model

The example model to be plotted comes from a finite mixture model with seven continuous indicators of a latent construct, each scaled to zero mean and unit variance. The intention of this plot is to visualize the means and standard error of each indicator across the latent classes.

  library(MplusAutomation)
  library(ggplot2)
  modelParams <- extractModelParameters("output_to_plot.out)$unstandardized
  modelParams <- subset(modelParams,
    paramHeader=="Means" &
      LatentClass != "Categorical.Latent.Variables",
    select=c("LatentClass", "param", "est", "se"))
  
  limits <- aes(ymax = est + se, ymin=est - se)
    
    fmmMeanPlot <- ggplot(modelParams, aes(x=param, y=est)) +
        geom_pointrange(limits) +
        scale_x_discrete("") +
        geom_hline(yintercept=0, color="grey50") +
        facet_grid(LatentClass ~ .) +
        theme_bw() +
        ylab("Mean Value") +
        coord_flip()
    print(fmmMeanPlot)

This relatively brief code snippet provides a relatively useful plot of this sort:

Figure. Example of graphing finite mixutre model results from Mplus using ggplot2.

Comparing summaries and parameters across models

An important part of interpreting results from latent variable models is the comparison of model fit indices and parameter estimates across related sets of models. For example, when an additional covariate is included in the model, how do the other parameter estimates change? The compareModels() function is designed to compare model fit indices and/or parameter estimates across two models.

This function also computes chi-square difference tests for nested models estimated with the ML, MLM, MLR, WLS, or WLSM estimators using the diffTest parameter.

To use compareModels(), I recommend using readModels() to extract various fit statistics and parameters from two or more models. You can pass in the results of extractModelSummaries() or extractModelParameters() to compareModels(), but the output will be limited to summaries or parameters, respectively.

Here is a brief example of how one might use compareModels().

parallelModels <- readModels("10_14_Harsh_SelfCon_Impul")

compareModels(parallelModels[["backport.from.grand.model.out"]],
  parallelModels[["backport.from.grand.model.slopesonw1.out"]],
  show = c("diff", "pdiff", "summaries", "unique"),
  equalityMargin = c(param = .05, pvalue = .02),
  sort = "type", diffTest = TRUE, showNS = FALSE)

Mplus model comparison
----------------------

------
Model 1:  10_14_Harsh_SelfCon_Impul/backport from grand model.out
Model 2:  10_14_Harsh_SelfCon_Impul/backport from grand model slopesonw1.out
------

Model Summary Comparison
------------------------

             m1                                    m2
Title        Harsh Impul Self-Control Trivariate   Harsh Impul Self-Control Trivariate
               LGCM with W1 Covariates - No Direct   LGCM with W1 Covariates - Direct
               Influence of W1 on Slopes and         Influence of W1 on Slopes and
               Intercepts                            Intercepts
Observations 2187                                  2187
Estimator    MLR                                   MLR
Parameters   87                                    103
LL           -69851.911                            -69820.77
AIC          139877.821                            139847.539
BIC          140372.876                            140433.639
ChiSqM_Value 280.781                               222.538
ChiSqM_DF    129                                   113
CFI          0.984                                 0.988
TLI          0.977                                 0.981
SRMR         0.025                                 0.022

  MLR Chi-Square Difference test for nested models based on loglikelihood
  -----------------------------------------------------------------------

  Difference Test Scaling Correction:  1.033437
  Chi-square difference:  60.2668
  Diff degrees of freedom:  16
  P-value:  0

  Note: The chi-square difference test assumes that these models are nested.
  It is up to you to verify this assumption.

  MLR Chi-Square Difference test for nested models
  --------------------------------------------

  Difference Test Scaling Correction: 1.03125
  Chi-square difference: 60.44
  Diff degrees of freedom: 16
  P-value: 0

Note: The chi-square difference test assumes that these models are nested.
  It is up to you to verify this assumption.

=========

Model parameter comparison
--------------------------


  Parameter estimates that differ between models (param. est. diff > 0.05)
  ----------------------------------------------
        paramHeader      param m1_est m2_est . m1_se m2_se . m1_est_se m2_est_se . m1_pval m2_pval
          SCC_S1.ON      IMP_I  0.155  0.073 | 0.036 0.050 |     4.255     1.478 |   0.000   0.139
         HPC_I.WITH      IMP_I  0.128  0.182 | 0.066 0.069 |     1.927     2.634 |   0.054   0.008
         SCC_I.WITH      HPC_I  1.391  1.455 | 0.189 0.183 |     7.372     7.949 |   0.000   0.000
         SCC_I.WITH      IMP_I  0.213  0.352 | 0.127 0.134 |     1.682     2.624 |   0.093   0.009
        SCC_S1.WITH     SCC_S2 -0.584 -0.649 | 0.160 0.157 |    -3.638    -4.141 |   0.000   0.000
         Intercepts A11XEMOTSP  4.523  4.592 | 0.568 0.568 |     7.964     8.084 |   0.000   0.000
         Intercepts      HPC_I  5.272  5.528 | 0.225 0.243 |    23.446    22.723 |   0.000   0.000
         Intercepts      HPC_S  0.725  0.583 | 0.092 0.092 |     7.865     6.358 |   0.000   0.000
         Intercepts      IMP_I -0.661 -0.544 | 0.155 0.167 |    -4.263    -3.252 |   0.000   0.001
         Intercepts      SCC_I  4.113  4.299 | 0.390 0.409 |    10.547    10.515 |   0.000   0.000
         Intercepts     SCC_S2  4.811  5.008 | 0.382 0.418 |    12.604    11.987 |   0.000   0.000
 Residual.Variances   A10CHARP  2.376  2.299 | 0.158 0.156 |    15.016    14.719 |   0.000   0.000
 Residual.Variances   A10CSELF  6.871  6.544 | 0.554 0.534 |    12.405    12.245 |   0.000   0.000
 Residual.Variances   A12CSELF  6.298  6.206 | 0.406 0.401 |    15.531    15.493 |   0.000   0.000
 Residual.Variances   A14CHARP  2.019  1.959 | 0.154 0.150 |    13.149    13.075 |   0.000   0.000
 Residual.Variances   A14CSELF  2.547  2.654 | 0.903 0.896 |     2.821     2.961 |   0.005   0.003
 Residual.Variances      HPC_I  2.634  2.764 | 0.166 0.162 |    15.846    17.079 |   0.000   0.000
 Residual.Variances      SCC_I  7.086  7.433 | 0.574 0.552 |    12.351    13.463 |   0.000   0.000
 Residual.Variances     SCC_S1  0.606  0.680 | 0.140 0.131 |     4.312     5.203 |   0.000   0.000


  P-values that differ between models (p-value diff > 0.02)
  -----------------------------------
 paramHeader param m1_est m2_est . m1_se m2_se . m1_est_se m2_est_se . m1_pval m2_pval
    HPC_S.ON IMP_I  0.049  0.028 | 0.011 0.015 |     4.243     1.894 |   0.000   0.058
   SCC_S1.ON IMP_I  0.155  0.073 | 0.036 0.050 |     4.255     1.478 |   0.000   0.139
  HPC_I.WITH IMP_I  0.128  0.182 | 0.066 0.069 |     1.927     2.634 |   0.054   0.008
  SCC_I.WITH IMP_I  0.213  0.352 | 0.127 0.134 |     1.682     2.624 |   0.093   0.009


  Parameters unique to model 1: 0
  -----------------------------

  None


  Parameters unique to model 2: 16
  -----------------------------

 paramHeader     param m2_est m2_se m2_est_se m2_pval
    HPC_S.ON W1PHARPUN  0.038 0.008     4.884   0.000
   SCC_S1.ON   W1PSELF  0.034 0.014     2.389   0.017


  14 filtered from output (fixed and/or n.s.)

    HPC_S.ON.IMPLW1, HPC_S.ON.W1PNEGEM, HPC_S.ON.W1PSELF, SCC_S1.ON.IMPLW1,
    SCC_S1.ON.W1PHARPUN, SCC_S1.ON.W1PNEGEM, SCC_S2.ON.IMPLW1,
    SCC_S2.ON.W1PHARPUN, SCC_S2.ON.W1PNEGEM, SCC_S2.ON.W1PSELF,
    IMP_S.ON.IMPLW1, IMP_S.ON.W1PHARPUN, IMP_S.ON.W1PNEGEM, IMP_S.ON.W1PSELF


==============

As can be discerned above, the example compareModels() call above compares two nested models from the list returned by readModels(). The show argument requests a comparison of parameter value difference ("diff"), p-value differences ("pdiff"), summary statistics ("summaries"), and parameters unique to each model ("unique").

Parameter estimate differences must exceed .05 to be displayed, and p-value differences must exceed .02 (this is specified by the equalityMargin argument).

The parameter comparisons are sorted by type (factor loadings, regressions, covariances, residual variances, etc.). A chi-square difference test is requested using diffTest=TRUE. And non-significant parameters are removed from the model comparison (showNS=FALSE).

Creating a group of models from an Mplus template file (createModels())

The third major focus of the MplusAutomation package is to provide tools that automate the process of creating input files for a related group of models. Perhaps the simplest example of a related group of models is latent class analysis, where one runs a certain model with different numbers of classes, but the input files are otherwise the same.

The createModels() function converts a single Mplus template file into a set of related Mplus input files. The template language is a simple extension of the Mplus language that allows dynamic values to be inserted into Mplus syntax files while reusing most of the code. The basic notion for template files is that multiple input files are created by looping over one or more variables (called “iterators”) and substituting specific values that change for each model, thereby allowing for the creation of related input files that share much of the code. As a basic example, one might iterate over a set of outcome variables (e.g., positive emotions, negative emotions, and conscientiousness) in a growth model where the dataset remains fixed, but the outcome variables change.

Mplus template syntax files are divided into two sections: the init section and the body section. The init section consists of definitions for variables to be inserted in the Mplus syntax, instructions for determining the filename and directory for created input files, and the variables to loop over to create multiple Mplus input files. The body section consists of Mplus syntax with template tags included where certain variables will change (e.g., the number of classes, names of outcome variables, etc.). Tags in Mplus syntax language are demarcated by double brackets. For example, the tag [[outcomeName]] requests that the value of provided in the init section be inserted into the Mplus syntax file.

Before we get into details, a particularly simple example may make the general concepts more tangible. This example is adapted from Mplus User’s Guide Example 7.3: LCA with binary latent class indicators using automatic starting values with random starts.

[[init]]
iterators = classes;
classes = 1:7;
filename = "[[classes]]-class LCA 7.3.inp";
outputDirectory = "C:/MplusAutomation/LCA Outputs";
[[/init]]

TITLE: this is an example of a LCA with binary
latent class indicators using automatic
starting values with random starts
DATA: FILE IS ex7.3.dat;
VARIABLE: NAMES ARE u1-u4 x1-x10;
USEVARIABLES = u1-u4;
CLASSES = c ([[classes]]);
CATEGORICAL = u1-u4;
AUXILIARY = x1-x10 (e);
ANALYSIS: TYPE = MIXTURE;
OUTPUT: TECH1 TECH8 TECH10;

The above template file instructs the createModels() function to loop over a variable called classes. The classes variable is defined as the integers from 1-7 (the colon in 1:7 indicates a sequence). Files should be named according to the number of classes, so when classes = 5, then filename = "[[classes]]-class LCA 7.3.inp"; will evaluate to 5-class LCA 7.3.inp. All files will be saved in the directory C:/MplusAutomation/LCA Outputs. Note that R`` uses forward slashes, not backslashes, to indicate directory paths.

In the body section of this simple example, the only thing that changes is the definition CLASSES = c([[classes]]). And because classes was defined as the integers from 1-7, 7 Mplus input files will be created by createModels() with the major difference being the number of classes. As createModels() loops over the classes variable, the current value of classes is inserted in the body section.

[[init]]
iterators = outcome;
outcome = 1:4;
outcomeNames#outcome = Conscien Extraver Agreeabl Openness;
filename = "CFA for [[outcomeNames#outcome]].inp";
outputDirectory = C:/;
[[/init]]

outcomeNames#outcome = "Conscientiousness Score" "Extraversion Score" "Agreeableness Score" "Openness Score";

[[init]]
iterators = outcome model classes;
outcome = 1:4;
model = 1:5;
classes = 2 3 4;
filename = "[[classes]]-class Model.inp";
outputDirectory = C:/Data/[[outcome]]/[[model]];
[[/init]]

[[init]]
iterators = outcome model;
outcome = 1:2;
model = 1:3;
outcomeName#outcome = Outcome1 Outcome2;
modelName#model = Model1 Model2 Model3;
outputDirectory = "C:/CFA/[[outcomeName#outcome]]/[[modelName#model]]";
filename="testfile.inp"
[[/init]]

CFA/
CFA/Outcome1/Model1
CFA/Outcome1/Model2
CFA/Outcome1/Model3
CFA/Outcome2/Model1
CFA/Outcome2/Model2
CFA/Outcome2/Model3

[[init]]
iterators = outcome model;
outcome = 1:5;
model = 1:3;
outcomeDirNames#outcome = Conscientiousness Extraversion Agreeableness Openness Neuroticism;
modelNames#model = Poisson "Negative Binomial" "Negative Binomial Hurdle";
filename = "[[modelNames#model]] Growth Model.inp";
outputDirectory = "Template Output/[[outcomeDirNames#outcome]];"
[[/init]]

[[init]]
iterators = model;
model = 1:3;
t1ModelNames#model = Outcome1_t1 Outcome2_t1 Outcome3_t1;
t2ModelNames#model = Outcome1_t2 Outcome2_t2 Outcome3_t2;
t3ModelNames#model = Outcome1_t3 Outcome2_t3 Outcome3_t3;
time1 = 0;
time2 = 1;
time3 = 3;
outputDirectory=output;
filename=[[modelOutcome#model]].inp;
[[/init]]
TITLE: [[modelOutcomes#model]] Latent Growth Curve Model
DATA:   FILE = "testData.dat";
VARIABLE:   NAMES ARE   id outcome1 outcome2 outcome3 covariate1;
MODEL: i s | [[t1ModelNames#model]]@[[time1]] [[t2ModelNames#model]]@[[time2]] [[t3ModelNames#model]]@[[time3]];
PLOT: TYPE=PLOT3;
    SERIES=[[t1ModelNames#model]] ([[time1]]) [[t2ModelNames#model]] ([[time2]]) [[t3ModelNames#model]] ([[time3]]);
    i s ON covariate1;

[[init]]
iterators = outcome;
outcome = 1:3;
modelName#outcome = "Outcome 1" "Outcome 2" "Outcome 3";
wave1Name#outcome = outcome1_w1 outcome2_w1 outcome3_w1;
wave2Name#outcome = outcome1_w2 outcome2_w2 outcome3_w2;
wave3Name#outcome = outcome1_w3 outcome2_w3 outcome3_w3;
outputDirectory=test;
filename="[[modelName#outcome]] Growth Model.inp";
[[/init]]
TITLE: [[modelName#outcome]] Growth Model
DATA:   FILE = "testData.dat";
VARIABLE:   NAMES ARE   id outcome1_w1 outcome1_w2 outcome1_w3
    outcome2_w1 outcome2_w2 outcome2_w3
    outcome3_w1 outcome3_w2 outcome3_w3;
    
    USEVARIABLES ARE [[wave1Name#outcome]] [[wave2Name#outcome]] [[wave3Name#outcome]];
MODEL:
    i s | [[wave1Name#outcome]]@0 [[wave2Name#outcome]]@1 [[wave3Name#outcome]]@2;

CLASSES = c ([[classes]]);

[[init]]
iterators=model;
model=1:5;
var1=test1;
var2=test2;
var3=test3;
[[/init]]
[[model > 2]]
COUNT ARE [[var1]] [[var2]] [[var3]];
[[/model > 2]]

[[modelNames#model != "Poisson GMM"]]
!Mplus code specific to non-Poisson GMM models here
[[/modelNames#model != "Poisson GMM"]]

[[init]]
iterators = outcome group model classes zeroclass;
outcome = 1:10;
group = 2 5;
model = 1:3;
classes = 1:4;
zeroclass = 1:2;

outcomenames#outcome = Paranoid Schizoid Schizotypal Antisocial
    Borderline Histrionic Narcissistic Avoidant Dependent OCPD;

groupnames#group = "Low Risk" "High Risk";

modelnames#model = "Normal LGCM" "Poisson GMM" "Poisson LCGA";

zeroclassnames#zeroclass = "" " with zero class";

#wave names are with respect to the outcome iterator
w1name#outcome = Paran1 Szoid1 Sztyp1 Anti1 Border1 Hist1
    Narc1 Avoid1 Depend1 OCPD1;
w2name#outcome = Paran2 Szoid2 Sztyp2 Anti2 Border2 Hist2
    Narc2 Avoid2 Depend2 OCPD2;
w3name#outcome = Paran3 Szoid3 Sztyp3 Anti3 Border3 Hist3
    Narc3 Avoid3 Depend3 OCPD3;

filename = "[[classes]]-class [[groupnames#group]] [[outcomenames#outcome]] [[modelnames#model]][[zeroclassnames#zeroclass]].inp";
outputDirectory = "PD GMM/[[outcomenames#outcome]]/[[groupnames#group]]/Unconditional_Models/[[modelnames#model]]";
[[/init]]

TITLE: [[classes]]-class [[outcomenames#outcome]] [[groupnames#group]] [[modelnames#model]] Unconditional Model[[zeroclassnames#zeroclass]]
DATA:   FILE = "personality_mplus.dat";
VARIABLE:   NAMES ARE
    id group sex age Paran1 Szoid1 Sztyp1 Anti1 Border1 Hist1 Narc1 Avoid1 Depend1
        OCPD1 PaAg1 Sadist1 SelfDef1 Paran2 Szoid2 Sztyp2 Anti2 Border2 Hist2 Narc2
        Avoid2 Depend2 OCPD2 Paran3 Szoid3 Sztyp3 Anti3 Border3 Hist3 Narc3 Avoid3
        Depend3 OCPD3;

    MISSING ARE .;
    USEVARIABLES ARE  [[w1name#outcome]] [[w2name#outcome]] [[w3name#outcome]];
    USEOBSERVATIONS ARE group EQ [[group]]; ![[groupnames#group]] Only
    
    [[model > 1]]
    COUNT ARE [[w1name#outcome]] [[w2name#outcome]] [[w3name#outcome]];
    [[/model > 1]]
    CLASSES = c ([[classes]]);

ANALYSIS:
    TYPE = MIXTURE;
    STARTS = 1000 10;
    K-1STARTS = 750 6;
    PROCESSORS = 4;
    [[model = 2]]
    ALGORITHM = INTEGRATION;
    [[/model = 2]]
MODEL:

%OVERALL%
Int Slope | [[w1name#outcome]]@0 [[w2name#outcome]]@0.97 [[w3name#outcome]]@2.77;

[[classes > 1]]
[[zeroclass = 2]]
!creates a class with all zeros at all time points
[[model=1]]
%c#2%
[Int@0 Slope@0];
Int@0 Slope@0;
[[/model=1]]
[[model>1]]
%c#2%
[Int@-15 Slope@-15];
Int@0 Slope@0;
[[/model>1]]
[[/zeroclass = 2]]
[[/classes > 1]]

PLOT:
    Type = PLOT3;
    Series =  [[w1name#outcome]] (0) [[w2name#outcome]] (0.97) [[w3name#outcome]] (2.77);
OUTPUT: TECH1 TECH4 [[model != 1]]TECH10 [[/model != 1]][[classes > 1]]TECH11 TECH14 [[/classes > 1]]STANDARDIZED RESIDUAL;

iterators = model;
model = 1:3;
A = a_test1 a_test2 "a_test3[[B#model]]";
B = "b_test1 [[A#model]]" "b_test2 [[A#model]]" "b_test3 [[A#model]]";

[[model > 1]]
[[classes != 1]]
var3@0; !Mplus code here
[[/classes]]
[[/model > 1]]