Using MPlus from R with MPlusAutomation
Sep 1, 2017
9 minute read

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

  1. Creating related groups of models
  2. Running batches
  3. 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)

I’m going to do something in what is maybe a bit of a clunky way: Hide the directory I’m saving all of the input and output files in by evaluating a line of code but not displaying it here, because it’s not necessary to see and is a pain to copy and paste. The key thing to 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 - that you’re saving these in.

We will use the iris dataset built-in to R. It can be loaded just by typing iris. 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’ll only use the continuous 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 have 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:

  1. One for which we estimate different means between two profiles
TITLE: iris LPA

DATA:
    File is iris.dat
    
VARIABLE: 
    Names are Sepal_Length Sepal_Width Petal_Length Petal_Width;
    Classes = c(2) ;
            
MODEL:
    %overall%
    Sepal_Length Sepal_Width Petal_Length Petal_Width; 
    
    %c#1%
    [Sepal_Length Sepal_Width Petal_Length Petal_Width];

    %c#2%
    [Sepal_Length Sepal_Width Petal_Length Petal_Width];
    
ANALYSIS: 
    Type is mixture;
            
OUTPUT:
    Tech11;
  1. 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 Sepal_Length Sepal_Width Petal_Length Petal_Width;
    Classes = c(2) ;
            
MODEL:
    %overall%
    Sepal_Length Sepal_Width Petal_Length Petal_Width; 
    Sepal_Length WITH Sepal_Width Petal_Length Petal_Width;
    Sepal_Width WITH Petal_Length Petal_Width;
    Petal_Length WITH Petal_Width;

    %c#1%
    [Sepal_Length Sepal_Width Petal_Length Petal_Width];
    
    %c#2%
    [Sepal_Length Sepal_Width Petal_Length Petal_Width];
    
ANALYSIS: 
    Type is mixture;
            
OUTPUT:
    Tech11;
  1. 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 Sepal_Length Sepal_Width Petal_Length Petal_Width;
    Classes = c(2) ;
            
MODEL:
    %c#1%
    Sepal_Length Sepal_Width Petal_Length Petal_Width; 

    Sepal_Length WITH Sepal_Width Petal_Length Petal_Width;
    Sepal_Width WITH Petal_Length Petal_Width;
    Petal_Length WITH Petal_Width;
    
    [Sepal_Length Sepal_Width Petal_Length Petal_Width];
    
    %c#2%
    Sepal_Length Sepal_Width Petal_Length Petal_Width; 
    
    Sepal_Length WITH Sepal_Width Petal_Length Petal_Width;
    Sepal_Width WITH Petal_Length Petal_Width;
    Petal_Length WITH Petal_Width;
    
    [Sepal_Length Sepal_Width Petal_Length Petal_Width];
    
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 SEPAL_LENG 5.006 0.049 102.032    0           1
## 2        Means SEPAL_WIDT 3.423 0.055  61.909    0           1
## 3        Means PETAL_LENG 1.471 0.026  55.788    0           1
## 4        Means PETAL_WIDT 0.250 0.016  15.938    0           1
## 5    Variances SEPAL_LENG 0.328 0.042   7.853    0           1
## 6    Variances SEPAL_WIDT 0.121 0.017   7.347    0           1
## 7    Variances PETAL_LENG 0.459 0.063   7.340    0           1
## 8    Variances PETAL_WIDT 0.123 0.013   9.126    0           1
## 9        Means SEPAL_LENG 6.265 0.068  92.358    0           2
## 10       Means SEPAL_WIDT 2.873 0.034  85.125    0           2
## 11       Means PETAL_LENG 4.911 0.085  57.798    0           2
## 12       Means PETAL_WIDT 1.678 0.043  38.643    0           2
## 13   Variances SEPAL_LENG 0.328 0.042   7.853    0           2
## 14   Variances SEPAL_WIDT 0.121 0.017   7.347    0           2
## 15   Variances PETAL_LENG 0.459 0.063   7.340    0           2
## 16   Variances PETAL_WIDT 0.123 0.013   9.126    0           2
m2$parameters[[1]][-nrow(m2$parameters[[1]]), ]
##      paramHeader      param   est    se  est_se pval LatentClass
## 1  SEPAL_LE.WITH SEPAL_WIDT 0.113 0.019   5.805    0           1
## 2  SEPAL_LE.WITH PETAL_LENG 0.305 0.050   6.104    0           1
## 3  SEPAL_LE.WITH PETAL_WIDT 0.114 0.019   6.112    0           1
## 4  SEPAL_WI.WITH PETAL_LENG 0.098 0.022   4.359    0           1
## 5  SEPAL_WI.WITH PETAL_WIDT 0.056 0.010   5.330    0           1
## 6  PETAL_LE.WITH PETAL_WIDT 0.193 0.024   8.175    0           1
## 7          Means SEPAL_LENG 5.006 0.049 101.442    0           1
## 8          Means SEPAL_WIDT 3.428 0.053  64.589    0           1
## 9          Means PETAL_LENG 1.462 0.024  60.137    0           1
## 10         Means PETAL_WIDT 0.246 0.015  16.674    0           1
## 11     Variances SEPAL_LENG 0.331 0.042   7.870    0           1
## 12     Variances SEPAL_WIDT 0.120 0.016   7.574    0           1
## 13     Variances PETAL_LENG 0.460 0.063   7.333    0           1
## 14     Variances PETAL_WIDT 0.123 0.013   9.137    0           1
## 15 SEPAL_LE.WITH SEPAL_WIDT 0.113 0.019   5.805    0           2
## 16 SEPAL_LE.WITH PETAL_LENG 0.305 0.050   6.104    0           2
## 17 SEPAL_LE.WITH PETAL_WIDT 0.114 0.019   6.112    0           2
## 18 SEPAL_WI.WITH PETAL_LENG 0.098 0.022   4.359    0           2
## 19 SEPAL_WI.WITH PETAL_WIDT 0.056 0.010   5.330    0           2
## 20 PETAL_LE.WITH PETAL_WIDT 0.193 0.024   8.175    0           2
## 21         Means SEPAL_LENG 6.262 0.066  94.947    0           2
## 22         Means SEPAL_WIDT 2.872 0.033  86.743    0           2
## 23         Means PETAL_LENG 4.906 0.082  59.719    0           2
## 24         Means PETAL_WIDT 1.676 0.042  39.652    0           2
## 25     Variances SEPAL_LENG 0.331 0.042   7.870    0           2
## 26     Variances SEPAL_WIDT 0.120 0.016   7.574    0           2
## 27     Variances PETAL_LENG 0.460 0.063   7.333    0           2
## 28     Variances PETAL_WIDT 0.123 0.013   9.137    0           2
m3$parameters[[1]][-nrow(m3$parameters[[1]]), ]
##      paramHeader      param   est    se  est_se  pval LatentClass
## 1  SEPAL_LE.WITH SEPAL_WIDT 0.097 0.022   4.469 0.000           1
## 2  SEPAL_LE.WITH PETAL_LENG 0.016 0.010   1.655 0.098           1
## 3  SEPAL_LE.WITH PETAL_WIDT 0.010 0.004   2.486 0.013           1
## 4  SEPAL_WI.WITH PETAL_LENG 0.011 0.008   1.418 0.156           1
## 5  SEPAL_WI.WITH PETAL_WIDT 0.009 0.005   1.763 0.078           1
## 6  PETAL_LE.WITH PETAL_WIDT 0.006 0.003   2.316 0.021           1
## 7          Means SEPAL_LENG 5.006 0.049 101.439 0.000           1
## 8          Means SEPAL_WIDT 3.428 0.053  64.591 0.000           1
## 9          Means PETAL_LENG 1.462 0.024  60.132 0.000           1
## 10         Means PETAL_WIDT 0.246 0.015  16.673 0.000           1
## 11     Variances SEPAL_LENG 0.122 0.022   5.498 0.000           1
## 12     Variances SEPAL_WIDT 0.141 0.033   4.267 0.000           1
## 13     Variances PETAL_LENG 0.030 0.007   4.222 0.000           1
## 14     Variances PETAL_WIDT 0.011 0.003   3.816 0.000           1
## 15 SEPAL_LE.WITH SEPAL_WIDT 0.121 0.027   4.467 0.000           2
## 16 SEPAL_LE.WITH PETAL_LENG 0.449 0.070   6.377 0.000           2
## 17 SEPAL_LE.WITH PETAL_WIDT 0.166 0.026   6.282 0.000           2
## 18 SEPAL_WI.WITH PETAL_LENG 0.141 0.033   4.330 0.000           2
## 19 SEPAL_WI.WITH PETAL_WIDT 0.079 0.015   5.295 0.000           2
## 20 PETAL_LE.WITH PETAL_WIDT 0.286 0.031   9.107 0.000           2
## 21         Means SEPAL_LENG 6.262 0.066  94.948 0.000           2
## 22         Means SEPAL_WIDT 2.872 0.033  86.743 0.000           2
## 23         Means PETAL_LENG 4.906 0.082  59.719 0.000           2
## 24         Means PETAL_WIDT 1.676 0.042  39.652 0.000           2
## 25     Variances SEPAL_LENG 0.435 0.059   7.332 0.000           2
## 26     Variances SEPAL_WIDT 0.110 0.017   6.442 0.000           2
## 27     Variances PETAL_LENG 0.675 0.086   7.822 0.000           2
## 28     Variances PETAL_WIDT 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.