Written by: Josh Rosenberg

Primary Source: Joshua M. Rosenberg – September 1, 2017

According to the MPlus website, the R package MPlusAutomation serves three purposes:

- Creating related groups of models
- Running batches
- Extracting and tabulating model parameters and test statistics.

Because modeling involves comparing related models, (partially) automating these is compelling. It can make it easier to use model results in subsequent analyses and can cut down on copy and pasting output or results between programs.

So I tried it and liked it; it’s well-designed. Here’s a little example exploring a set of models for which one aspect of its specification changes between models. This example uses built-in data, so you should be able to do everything here in this example, **as long as you have purchased and installed MPlus**.

First, let’s load the package, which we can install from CRAN using `install.packages("MPlusAutomation")`

:

```
# install.packages("MPlusAutomation")
library(MplusAutomation)
```

`## Loading required package: methods`

I’m going to do something in what is maybe a bit of a clunky way: Hide the directory I am saving all of the input and output files in by evaluating a line of code but not displaying it here. The key thing tto know is that the variable `the_dir`

represents the working directory I’m saving these in; you should replace it with the directory, whether it’s a Dropbox folder, on your Desktop, or anywhere else, you’re saving these in. You can navigate to the directory and in the R Studio Files window, and then click “More” and then “Set as Working Directory” to easily set and use the specific working directory.

We will use the `iris`

dataset. Here is how it looks:

`head(iris)`

```
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
```

Next, we’ll use the super handy `prepareMPlusData()`

function to prepare a data file in the format MPlus uses, namely, a tab-separated `.dat`

file with no column names.

`prepareMplusData(iris[, -5], paste0(the_dir, "iris.dat"))`

Let’s unpack what we’re doing here.

- If you recall the contents of
`iris`

, you’ll note that the fifth column is the name of the species. Our goal in this analysis is to use a general mixture modeling or Latent Profile Analysis (APA) approach in MPlus, and so we’l only use the continous variables as predictor variables. - the somewhat-inelegant bit of code
`paste0(the_dir, "iris.dat")`

is basically pasting together are not-so-secretly-hid working directory with a name we chose for this specific`.dat`

file, namely,`iris.dat`

. In short, this says**save the prepared data file in this particular folder with this particular name**. Again, you’ll hvae to change`the_dir`

to be the folder (or working directory) that you chose, set, and are using.

Now, we have to specify the models using MPlus syntax.

We’ll specify three mixture models. The trick is that we’ll **save each of the following models (either in a .txt file or in MPlus style using a .inp.file type) with different names in the working directory we saved the data file to**.

Here we go:

- One for which we estimate different means between two profiles

```
TITLE: iris LPA
DATA:
File is iris.dat
VARIABLE:
Names are x1, x2, x3, x4;
Classes = c(2) ;
MODEL:
%overall%
x1 x2 x3 x4;
%c#1%
[x1-x4];
%c#2%
[x1-x4];
ANALYSIS:
Type is mixture;
OUTPUT:
Tech11;
```

- One for which we estimate different means between the two profiles
**and**the model is specified to estimate the correlation (or covariance) for the variables

```
TITLE: iris LPA
DATA:
File is iris.dat
VARIABLE:
Names are x1, x2, x3, x4;
Classes = c(2) ;
MODEL:
%overall%
x1 x2 x3 x4;
x1 WITH x2-x4;
x2 WITH x3-x4;
x3 WITH x4;
%c#1%
[x1-x4];
%c#2%
[x1-x4];
ANALYSIS:
Type is mixture;
OUTPUT:
Tech11;
```

- One for which we estimate different means
**and**the model is specified to*different*covariances (and variable variances) between the two profiles

```
TITLE: iris LPA
DATA:
File is iris.dat
VARIABLE:
Names are x1, x2, x3, x4;
Classes = c(2) ;
MODEL:
%c#1%
x1 x2 x3 x4;
x1 WITH x2-x4;
x2 WITH x3-x4;
x3 WITH x4;
[x1-x4];
%c#2%
x1 x2 x3 x4;
x1 WITH x2-x4;
x2 WITH x3-x4;
x3 WITH x4;
[x1-x4];
ANALYSIS:
Type is mixture;
OUTPUT:
Tech11;
```

Again, each of those models has to be saved in the working directory we saved the data in.

Now we can run the models using `runModels()`

, which runs the model in MPlus. You can see here what I named each of the three files:

```
# Model 1
runModels(paste0(the_dir, "2-iris-LPA_means.inp"))
# Model 2
runModels(paste0(the_dir, "2-iris-LPA_means_correlated.inp"))
# Model 3
runModels(paste0(the_dir, "2-iris-LPA_means_correlated_free_variances.inp"))
```

Now we’re in business and can look at the output using `readModels()`

:

```
m1 <- readModels(paste0(the_dir, "2-iris-LPA_means.out"))
m2 <- readModels(paste0(the_dir, "2-iris-LPA_means_correlated.out"))
m3 <- readModels(paste0(the_dir, "2-iris-LPA_means_correlated_free_variances.out"))
```

We can now inspect the fit statistics and other summary information for the three models:

`m1$summaries$BIC`

`## [1] 1042.968`

`m2$summaries$BIC`

`## [1] 688.097`

`m3$summaries$BIC`

`## [1] 574.018`

And can examine parameters, as well (ignore the `nrow(...)`

for now; the last row is not necessary, and this – clunkily, but for now, easily – does not print it):

`m1$parameters[[1]][-nrow(m1$parameters[[1]]), ]`

```
## paramHeader param est se est_se pval LatentClass
## 1 Means X1 5.006 0.049 102.032 0 1
## 2 Means X2 3.423 0.055 61.909 0 1
## 3 Means X3 1.471 0.026 55.788 0 1
## 4 Means X4 0.250 0.016 15.938 0 1
## 5 Variances X1 0.328 0.042 7.853 0 1
## 6 Variances X2 0.121 0.017 7.347 0 1
## 7 Variances X3 0.459 0.063 7.340 0 1
## 8 Variances X4 0.123 0.013 9.126 0 1
## 9 Means X1 6.265 0.068 92.358 0 2
## 10 Means X2 2.873 0.034 85.125 0 2
## 11 Means X3 4.911 0.085 57.798 0 2
## 12 Means X4 1.678 0.043 38.643 0 2
## 13 Variances X1 0.328 0.042 7.853 0 2
## 14 Variances X2 0.121 0.017 7.347 0 2
## 15 Variances X3 0.459 0.063 7.340 0 2
## 16 Variances X4 0.123 0.013 9.126 0 2
```

`m2$parameters[[1]][-nrow(m2$parameters[[1]]), ]`

```
## paramHeader param est se est_se pval LatentClass
## 1 X1.WITH X2 0.113 0.019 5.805 0 1
## 2 X1.WITH X3 0.305 0.050 6.104 0 1
## 3 X1.WITH X4 0.114 0.019 6.112 0 1
## 4 X2.WITH X3 0.098 0.022 4.359 0 1
## 5 X2.WITH X4 0.056 0.010 5.330 0 1
## 6 X3.WITH X4 0.193 0.024 8.175 0 1
## 7 Means X1 5.006 0.049 101.442 0 1
## 8 Means X2 3.428 0.053 64.589 0 1
## 9 Means X3 1.462 0.024 60.137 0 1
## 10 Means X4 0.246 0.015 16.674 0 1
## 11 Variances X1 0.331 0.042 7.870 0 1
## 12 Variances X2 0.120 0.016 7.574 0 1
## 13 Variances X3 0.460 0.063 7.333 0 1
## 14 Variances X4 0.123 0.013 9.137 0 1
## 15 X1.WITH X2 0.113 0.019 5.805 0 2
## 16 X1.WITH X3 0.305 0.050 6.104 0 2
## 17 X1.WITH X4 0.114 0.019 6.112 0 2
## 18 X2.WITH X3 0.098 0.022 4.359 0 2
## 19 X2.WITH X4 0.056 0.010 5.330 0 2
## 20 X3.WITH X4 0.193 0.024 8.175 0 2
## 21 Means X1 6.262 0.066 94.947 0 2
## 22 Means X2 2.872 0.033 86.743 0 2
## 23 Means X3 4.906 0.082 59.719 0 2
## 24 Means X4 1.676 0.042 39.652 0 2
## 25 Variances X1 0.331 0.042 7.870 0 2
## 26 Variances X2 0.120 0.016 7.574 0 2
## 27 Variances X3 0.460 0.063 7.333 0 2
## 28 Variances X4 0.123 0.013 9.137 0 2
```

`m3$parameters[[1]][-nrow(m3$parameters[[1]]), ]`

```
## paramHeader param est se est_se pval LatentClass
## 1 X1.WITH X2 0.097 0.022 4.469 0.000 1
## 2 X1.WITH X3 0.016 0.010 1.655 0.098 1
## 3 X1.WITH X4 0.010 0.004 2.486 0.013 1
## 4 X2.WITH X3 0.011 0.008 1.418 0.156 1
## 5 X2.WITH X4 0.009 0.005 1.763 0.078 1
## 6 X3.WITH X4 0.006 0.003 2.316 0.021 1
## 7 Means X1 5.006 0.049 101.439 0.000 1
## 8 Means X2 3.428 0.053 64.591 0.000 1
## 9 Means X3 1.462 0.024 60.132 0.000 1
## 10 Means X4 0.246 0.015 16.673 0.000 1
## 11 Variances X1 0.122 0.022 5.498 0.000 1
## 12 Variances X2 0.141 0.033 4.267 0.000 1
## 13 Variances X3 0.030 0.007 4.222 0.000 1
## 14 Variances X4 0.011 0.003 3.816 0.000 1
## 15 X1.WITH X2 0.121 0.027 4.467 0.000 2
## 16 X1.WITH X3 0.449 0.070 6.377 0.000 2
## 17 X1.WITH X4 0.166 0.026 6.282 0.000 2
## 18 X2.WITH X3 0.141 0.033 4.330 0.000 2
## 19 X2.WITH X4 0.079 0.015 5.295 0.000 2
## 20 X3.WITH X4 0.286 0.031 9.107 0.000 2
## 21 Means X1 6.262 0.066 94.948 0.000 2
## 22 Means X2 2.872 0.033 86.743 0.000 2
## 23 Means X3 4.906 0.082 59.719 0.000 2
## 24 Means X4 1.676 0.042 39.652 0.000 2
## 25 Variances X1 0.435 0.059 7.332 0.000 2
## 26 Variances X2 0.110 0.017 6.442 0.000 2
## 27 Variances X3 0.675 0.086 7.822 0.000 2
## 28 Variances X4 0.179 0.018 10.148 0.000 2
```

Cool!

An especially powerful feature of `MPlusAutomation`

is the ability to use `runModels()`

and `readModels()`

in conjunction with the function `createModels()`

, because these functions allow you to specify an entire folder, in addition to a specific model or output file.

The `createModels()`

function creates a set of models using a template. Here is an example that would create models with different numbers of profiles, from `1`

to `9`

. Here is doing that for the model with different means between profiles specified (model, 1 above):

```
[[init]]
iterators = classes;
classes = 1:9;
filename = "[[classes]]-iris-LPA.inp";
outputDirectory = the_dir;
[[/init]]
TITLE: iris LPA
DATA:
File is iris.dat
VARIABLE:
Names are x1 x2 x3 x4;
Classes = c([[classes]]) ;
MODEL:
%overall%
x1 x2 x3 x4;
[x1-x4];
ANALYSIS:
Type is mixture;
OUTPUT:
Tech11;
```

We would then run the following series of functions after saving the above-specified model as “iris_lpa_template.inp” (note that we save the output of `readModels()`

to a list object):

```
createModels(paste0(the_dir, "iris_lpa_template.inp"))
runModels(the_dir)
models_list <- readModels(the_dir)
```

We can then inspect specific models using list-indexing: The first model is saved as `models_list[[1]]`

, for example. Or, we can print the output for all of the models by typing `models_list`

.

You can learn more about MPLusAutomation here.

#### Latest posts by Josh Rosenberg (see all)

- An R package for sensitivity analysis (konfound) - April 16, 2018
- Explorations in Markov Chain Monte Carlo – comparing results from MCMCglmm and lme4 - March 31, 2018
- Find the top rail-trails in each state using mixed effects models - February 28, 2018