This vignette will demonstrate latent variable modeling and confirmatory analysis via the common factor model. Factor analysis is a rudimentary application of OpenMx, and hopefully will serve as an easily understandable introduction to the data structures and workflow patterns necessary to define and fit a model using OpenMx.

OpenMX handles model definition through the MxModel class, instantiated via the mxModel() function. Model definition requires input data for analysis and the model specification itself. Input data can be either in the form of a set of variable observations or by having a covariance matrix and associated variable means. The model specification can either be produced via specifying a path diagram or, for advanced users, the matrices to define the desired RAM or LISREL model.

We’ll walk through the construction of a single and multiple factor model first through specifying the path diagram, and then an alternate form of model specification manually defining the matrices in a RAM-type model. Additionally, we will demonstrate how to specify data for analysis using either raw data or a covariance matrix of manifest variables. Scripts for a variety of implementations are located within the /demo/ directory of the OpenMx GitHub repository and are named in accordance with their specific implementation (eg. “OneFactorModel_PathRaw.R”, “TwoFactorModel_MatrixCov.R”, etc.)

Common Factor Model

The common factor model supposes that observed variables measure or indicate the same latent variable. While there are a number of exploratory approaches that do not assume the number of latent factor(s), this example uses a confirmatory approach that assumes a single factor. A path diagram verions of the common factor model for a set of variables \(x_{1}-x_{6}\) are given below.

\begin{eqnarray} x{ij} = \mu{j} + \lambda{j} * \eta{i} + \epsilon_{ij} \end{eqnarray}

One Factor Model

Multiple Factor Model

The common factor model can be extended to include multiple latent variables. The model for any person and path diagram of the common factor model for a set of variables \(x_{1}-x_{3}\) and \(y_{1}-y_{3}\) are given below.

\begin{eqnarray} x{ij} = \mu{j} + \lambda{j} * \eta{1i} + \epsilon{ij}\ y{ij} = \mu{j} + \lambda{j} * \eta{2i} + \epsilon{ij} \end{eqnarray}

Two Factor Model

While 19 parameters are displayed in the equation and path diagram above (six manifest variances, six manifest means, six factor loadings and one factor variance), we must constrain either the factor variance or one factor loading to a constant to identify the model and scale the latent variable. As such, this model contains 18 parameters. The means and covariance matrix for six observed variables contain 27 degrees of freedom, and thus our model contains 9 degrees of freedom.

Input Data

Our first step to running this model is to include the data to be analyzed. The MxData class instantiated through mxData() is used for representing data within the MxModel.

As mentioned, a raw set of observations on variables can be provided as data. The data for this example is available from the myFADataRaw dataset, which contains 500 observations on nine observed variables. In accordance with the path diagram above, we specify the first six variables for analysis in the model.


dataRawOneFactor <- mxData( observed=myFADataRaw[,c("x1","x2","x3","x4","x5","x6")], type="raw" )

Alternatively, our data could be specificed via a covariance matrix with variable means and number of observations. A type argument of "cov" is specified in the below example with a covariance matrix (other inputs being the same).

    c(0.997, 0.642, 0.611, 0.672, 0.637, 0.677,
      0.642, 1.025, 0.608, 0.668, 0.643, 0.676,
      0.611, 0.608, 0.984, 0.633, 0.657, 0.626,
      0.672, 0.668, 0.633, 1.003, 0.676, 0.665,
      0.637, 0.643, 0.657, 0.676, 1.028, 0.654,
      0.677, 0.676, 0.626, 0.665, 0.654, 1.020),

myFADataMeans <- c(2.988, 3.011, 2.986, 3.053, 3.016, 3.010)
names(myFADataMeans) <- c("x1","x2","x3","x4","x5","x6")

dataCovOneFactor  <- mxData( observed=myFADataCov, type="cov", numObs=500,
                        mean=myFADataMeans )

Input Data, Multiple Factors

The only difference in setting up data for multiple factor analysis is to label the manifest variables appropriately. The below code reads in the first six variables from the myFADataRaw dataset once again, but this time the labels for \(y_n\) observed variables for the second factor, \(F_2\) are changed to match the two factor path diagram above. Variable naming is arbitrary and does not affect the model structure.

dataRawTwoFactor <- mxData( observed=myFADataRaw[,c("x1","x2","x3","y1","y2","y3")], type="raw" )

Setting up a covariance data input for multiple factors involves similar changes from the single factor case. The below code block sets up a covariance matrix of 9 observed variables, 6 of which are labeled as \(x_n\), and 3 of which are \(y_n\). For the MxData selection, \(x_{1}-x_{3}\) and \(y_{1}-y_{3}\) are used in line with our desired path diagram.

twoFactorCov <- matrix(
      c(0.997, 0.642, 0.611, 0.342, 0.299, 0.337,
        0.642, 1.025, 0.608, 0.273, 0.282, 0.287,
        0.611, 0.608, 0.984, 0.286, 0.287, 0.264,
        0.342, 0.273, 0.286, 0.993, 0.472, 0.467,
        0.299, 0.282, 0.287, 0.472, 0.978, 0.507,
        0.337, 0.287, 0.264, 0.467, 0.507, 1.059),
          c("x1", "x2", "x3", "y1", "y2", "y3"),
          c("x1", "x2", "x3", "y1", "y2", "y3")),

twoFactorMeans <- c(2.988, 3.011, 2.986, 2.955, 2.956, 2.967)
names(twoFactorMeans) <- c("x1", "x2", "x3", "y1", "y2", "y3")
# Prepare Data
# -----------------------------------------------------------------------------

dataCovTwoFactor  <- mxData( observed=twoFactorCov, type="cov", numObs=500, means=twoFactorMeans )

Model Specification

Path Diagram Specification

One method of model specification is through defining the model via the existing path diagram. In this instance, our model instantiation will take as arguments the MxData defined previously, as well as lists of manifest variables & latent variables, and an arbitrary number of MxPath objects defining the paths between these variables. First, we will start by defining our manifest and latent variables in accordance to the path diagram.

latentVars <- "F1"

Next, we define our various paths between variables. The mxPath() function has a number of parameters for instantiating a similar set of paths in a single call, and all of these can be grouped together within the mxModel() call later.

First off, we define the latent variable variance path. The from field without a specified to value evaluates as if the paths go back to their origin, as we want with our variances. The arrows field of 2 indicates a double-headed path in this instance. This mxPath() call effectively represents one path, the variance on the latent variable \(F_1\), which we default pre-fitting to 1 and label "varF1" The free field of TRUE indicates that this value will be optimized once we fit the model.

# latent variance
latVar       <- mxPath( from="F1", arrows=2,
                        free=TRUE, values=1, labels ="varF1" )

Additionally, we need residual variance paths for the observed variables. These can all be created in one mxPath() call by specifying a list in from, as well as initial values and labels for these variances.

# residual variances
resVars      <- mxPath( from=c("x1","x2","x3","x4","x5","x6"), arrows=2,
                        free=TRUE, values=c(1,1,1,1,1,1),
                        labels=c("e1","e2","e3","e4","e5","e6") )

Next come the factor loadings. These are specified as directed paths (regressions) of the manifest variables on the latent variable "F1". As we have to scale the latent variable, the first factor loading has been given a fixed value of one by setting the first elements of the free and values arguments to FALSE and 1, respectively. Alternatively, the latent variable could have been scaled by fixing the factor variance to 1 in the previous mxPath() function and freely estimating all factor loadings. The five factor loadings that are freely estimated are all given starting values of 1 and labels "l2" through "l6".

# factor loadings
facLoads     <- mxPath( from="F1", to=c("x1","x2","x3","x4","x5","x6"), arrows=1,
                        free=c(FALSE,TRUE,TRUE,TRUE,TRUE,TRUE), values=c(1,1,1,1,1,1),
                        labels =c("l1","l2","l3","l4","l5","l6") )

Lastly, we must specify the mean structure for this model. As there are a total of seven variables in this model (six manifest and one latent), we have the potential for seven means. However, we must constrain at least one mean to a constant value, as there is not sufficient information to yield seven mean and intercept estimates from the six observed means. The six observed variables receive freely estimated intercepts, while the factor mean is fixed to a value of zero in the code below.

# means
means        <- mxPath( from="one", to=c("x1","x2","x3","x4","x5","x6","F1"), arrows=1,
                        free=c(T,T,T,T,T,T,FALSE), values=c(1,1,1,1,1,1,0),
                        labels =c("meanx1","meanx2","meanx3","meanx4","meanx5","meanx6",NA) )

Finally, we insantiate our unfitted model using mxModel(). This model begins with a name (“Common Factor Model Path Specification”) for the model and a type="RAM" argument. The name for the model may be omitted, or may be specified in any other place in the model using the name argument. Including type="RAM" allows the mxModel function to interpret the mxPath functions that follow and turn those paths into the corresponding RAM specification matrices. The raw data specification is included as the model data.

oneFactorPathModel <- mxModel("Common Factor Model Path Specification", type="RAM",
                        manifestVars=c("x1","x2","x3","x4","x5","x6"), latentVars="F1",
                        dataRawOneFactor, resVars, latVar, facLoads, means)

Path Specification for Multiple Factors

We will focus on the changes this model makes to the mxPath functions and their relevant inputs. To start, the last three variables of our manifestVars argument have changed from \(x_4, x_5, x_6\) to \(y_1, y_2, y_3\), and latentVars now includes a second variable, F_2. The mxPath functions which detail residual variances and intercepts accommodate the changes in manifest and latent variables but carry out identical functions to the common factor model.

manifestVars <- c("x1","x2","x3","y1","y2","y3")
latentVars <- c("F1","F2")

# residual variances
resVars      <- mxPath( from=c("x1", "x2", "x3", "y1", "y2", "y3"), arrows=2,
                        free=TRUE, values=c(1,1,1,1,1,1),
                        labels=c("e1","e2","e3","e4","e5","e6") )
# means
means        <- mxPath( from="one", to=c("x1","x2","x3","y1","y2","y3","F1","F2"),
                        free=c(T,T,T,T,T,T,F,F), values=c(1,1,1,1,1,1,0,0),
                                 "meany1","meany2","meany3",NA,NA) )

The remaining mxPath functions provide some changes to the model. The latVars mxPath function specifies the variances and covariance of the two latent variables. Like previous examples, we’ve omitted the to argument for this set of two-headed paths. Unlike previous examples, we’ve set the connect argument to unique.pairs, which creates all unique paths between the variables. As omitting the to argument is identical to putting identical variables in the from and to arguments, we are creating all unique paths from and to our two latent variables. This results in three paths: from \(F_1\) to \(F_1\) (the variance of \(F_1\)), from \(F_1\) to \(F_2\) (the covariance of the latent variables), and from \(F_2\) to \(F_2\) (the variance of \(F_2\)).

# latent variances and covariance
latVars      <- mxPath( from=c("F1","F2"), arrows=2, connect="unique.pairs",
                        free=TRUE, values=c(1,.5,1), labels=c("varF1","cov","varF2") )

You can check your understanding by printing out the paths to the console.

## mxPath
## F1 <-> F1 [value=1, free=TRUE, label='varF1']
## F1 <-> F2 [value=0.5, free=TRUE, label='cov']
## F2 <-> F2 [value=1, free=TRUE, label='varF2']

The final two mxPath functions define the factor loadings for each of the latent variables. We’ve split these loadings into two functions, one for each latent variable. The first loading for each latent variable is fixed to a value of one, just as in the previous example.

# factor loadings for x variables
facLoadsX    <- mxPath( from="F1", to=c("x1","x2","x3"), arrows=1,
                        free=c(F,T,T), values=c(1,1,1), labels=c("l1","l2","l3") )
# factor loadings for y variables
facLoadsY    <- mxPath( from="F2", to=c("y1","y2","y3"), arrows=1,
                        free=c(F,T,T), values=c(1,1,1), labels=c("l4","l5","l6") )

Finally, we instantiate the mxModel object. Specifying the two factor model is virtually identical to the single factor case, with the only change being our manifest and latent variable labels.

twoFactorPathModel <- mxModel("Two Factor Model Path Specification", type="RAM",
                        manifestVars=c("x1", "x2", "x3", "y1", "y2", "y3"),
                        dataRawTwoFactor, resVars, latVars, facLoadsX, facLoadsY, means)

Matrix Specification

Alternatively to specifying the paths for an MxModel, we have the option of manually specifying the matrices in a RAM model. In this case, MxModel specification is carried out using mxMatrix functions to create matrices. In the following section, we will go through the necessary setup to do manual model setup for the same one and two factor examples as path specification.

The \(A\) matrix specifies all of the asymmetric paths or regressions in our model. In the common factor model, these parameters are the factor loadings. This matrix is square, and contains as many rows and columns as variables in the model (manifest and latent, typically in that order). Regressions are specified in the A matrix by placing a free parameter in the row of the dependent variable and the column of independent variable.

The common factor model requires that one parameter (typically either a factor loading or factor variance) be constrained to a constant value. In our model, we will constrain the first factor loading to a value of 1, and let all other loadings be freely estimated. All factor loadings have a starting value of one and labels of "l1" - "l6".

# asymmetric paths
matrA        <- mxMatrix( type="Full", nrow=7, ncol=7,
                          free=  c(F,F,F,F,F,F,F,
                          byrow=TRUE, name="A" )

The second matrix in a RAM model is the \(S\) matrix, which specifies the symmetric or covariance paths in our model. This matrix is symmetric and square, and contains as many rows and columns as variables in the model (manifest and latent, typically in that order). The symmetric paths in our model consist of six residual variances and one factor variance. All of these variances are given starting values of one and labels "e1" - "e6" and "varF1".

# symmetric paths
matrS        <- mxMatrix( type="Symm", nrow=7, ncol=7,
                          free=  c(T,F,F,F,F,F,F,
                          labels=c("e1",NA,  NA,  NA,  NA,  NA,  NA,
                                   NA, "e2", NA,  NA,  NA,  NA,  NA,
                                   NA,  NA, "e3", NA,  NA,  NA,  NA,
                                   NA,  NA,  NA, "e4", NA,  NA,  NA,
                                   NA,  NA,  NA,  NA, "e5", NA,  NA,
                                   NA,  NA,  NA,  NA,  NA, "e6", NA,
                                   NA,  NA,  NA,  NA,  NA,  NA, "varF1"),
                          byrow=TRUE, name="S" )

The third matrix in our RAM model is the \(F\) or filter matrix. Our data contains six observed variables, but the \(A\) and \(S\) matrices contain seven rows and columns. For our model to define the covariances present in our data, we must have some way of projecting the relationships defined in the \(A\) and \(S\) matrices onto our data. The \(F\) matrix filters the latent variables out of the expected covariance matrix, and can also be used to reorder variables.

The \(F\) matrix will always contain the same number of rows as manifest variables and columns as total (manifest and latent) variables. If the manifest variables in the \(A\) and \(S\) matrices precede the latent variables and are in the same order as the data, then the \(F\) matrix will be the horizontal adhesion of an identity matrix and a zero matrix. This matrix never contains free parameters, and is made with the mxMatrix function below.

# filter matrix
matrF        <- mxMatrix( type="Full", nrow=6, ncol=7,
                          byrow=TRUE, name="F" )

The last matrix of our model is the \(M\) matrix, which defines the means and intercepts for our model. This matrix describes all of the regressions on the constant in a path model, or the means conditional on the means of exogenous variables. This matrix contains a single row, and one column for every manifest and latent variable in the model. In our model, the latent variable has a constrained mean of zero, while the manifest variables have freely estimated means, labeled "meanx1" through "meanx6".

# means
matrM        <- mxMatrix( type="Full", nrow=1, ncol=7,
                          name="M" )

The final parts of this model are the expectation function and the fit function. The expectation defines how the specified matrices combine to create the expected covariance matrix of the data. The fit defines how the expectation is compared to the data to create a single scalar number that is minimized. In a RAM specified model, the expected covariance matrix is defined as:

\begin{eqnarray} ExpCovariance = F * (I - A){-1} * S * ((I - A){-1})' * F' \end{eqnarray}

The expected means are defined as:

\begin{eqnarray} ExpMean = F * (I - A){-1} * M \end{eqnarray}

The free parameters in the model can then be estimated using maximum likelihood for covariance and means data, and full information maximum likelihood for raw data. Although users may define their own expected covariance matrices using mxExpectationNormal and other functions in OpenMx, the mxExpectationRAM function computes the expected covariance and means matrices when the \(A\), \(S\), \(F\) and \(M\) matrices are specified. The \(M\) matrix is required both for raw data and for covariance data that includes a means vector. The mxExpectationRAM function takes four arguments, which are the names of the \(A\), \(S\), \(F\) and \(M\) matrices in your model. The mxFitFunctionML yields maximum likelihood estimates of structural equation models. It uses full information maximum likelihood when the data are raw.

exp          <- mxExpectationRAM("A","S","F","M",
                                  dimnames=c(manifestVars, latentVars))
funML        <- mxFitFunctionML()

We now have all of the prerequisite structures for defining our model. As opposed to the path specification, the model now includes data of observed variables (i.e., raw data or manifest variable covariances in MxData form), model matrices, an expectation function, and a fit function. As such, the model has all the required elements to define the expected covariance matrix and estimate parameters.

oneFactorMatrixModel <- mxModel("Common Factor Model Matrix Specification",
                          dataRawOneFactor, matrA, matrS, matrF, matrM, exp, funML)

Matrix Specification for Multiple Factors

To account for the second latent variable \(F_2\), we must modify the constituent RAM matricies.

The \(A\) matrix, shown below, is used to specify the regressions of the manifest variables on the factors. The first three manifest variables ("x1"-"x3") are regressed on "F1", and the second three manifest variables ("y1"-"y3") are regressed on "F2". We must again constrain the model to identify and scale the latent variables, which we do by constraining the first loading for each latent variable to a value of one.

# asymmetric paths
matrA        <- mxMatrix( type="Full", nrow=8, ncol=8,
                          free=  c(F,F,F,F,F,F,F,F,
                          byrow=TRUE, name="A" )

The \(S\) matrix has an additional row and column, and two additional parameters. For the two factor model, we must add a variance term for the second latent variable and a covariance between the two latent variables.

# symmetric paths
matrS        <- mxMatrix( type="Symm", nrow=8, ncol=8,
                          free=  c(T,F,F,F,F,F,F,F,
                          labels=c("e1",NA,  NA,  NA,  NA,  NA,  NA,  NA,
                                   NA, "e2", NA,  NA,  NA,  NA,  NA,  NA,
                                   NA,  NA, "e3", NA,  NA,  NA,  NA,  NA,
                                   NA,  NA,  NA, "e4", NA,  NA,  NA,  NA,
                                   NA,  NA,  NA,  NA, "e5", NA,  NA,  NA,
                                   NA,  NA,  NA,  NA,  NA, "e6", NA,  NA,
                                   NA,  NA,  NA,  NA,  NA,  NA,"varF1","cov",
                                   NA,  NA,  NA,  NA,  NA,  NA,"cov","varF2"),
                          byrow=TRUE, name="S" )

The F and M matrices contain only minor changes. The F matrix is now of order 6x8, but the additional column is simply a column of zeros. The M matrix contains an additional column (with only a single row), which contains the mean of the second latent variable. As this model does not contain a parameter for that latent variable, this mean is constrained to zero.

matrF        <- mxMatrix( type="Full", nrow=6, ncol=8,
                          byrow=TRUE, name="F" )
matrM        <- mxMatrix( type="Full", nrow=1, ncol=8,
                          name="M" )

Finally, we specify our expectation and fit functions as well as instantiating the MxModel itself. These are completely identical to the one factor model, only with the data variable changed from dataRawOneFactor to dataRawTwoFactor.

exp          <- mxExpectationRAM("A","S","F","M",
                                  dimnames=c(manifestVars, latentVars))
funML        <- mxFitFunctionML()

twoFactorMatrixModel <- mxModel("Two Factor Model Matrix Specification",
                          dataRawTwoFactor, matrA, matrS, matrF, matrM, exp, funML)

Running Model & Results

Our model can now be run using the mxRun function. mxRun returns an MxModel structured on its input MxModel with the free parameters fitted. This fitting does not modify the input model, and as such mxRun must assign its output to a new object.

twoFactorMatrixFit <- mxRun(twoFactorMatrixModel)
## Running Two Factor Model Matrix Specification with 19 parameters

summary(twoFactorMatrixFit) will provide us with the estimation and standard errors on all free paramaters, along with several model statistics and criteria for assessing model fit.

## Summary of Two Factor Model Matrix Specification 
## free parameters:
##      name matrix row col  Estimate  Std.Error A
## 1      l2      A   2   7 0.9723169 0.05914081  
## 2      l3      A   3   7 0.9313485 0.05777022  
## 3      l5      A   5   8 1.0489157 0.09346097  
## 4      l6      A   6   8 1.0530553 0.09542648  
## 5      e1      S   1   1 0.3348589 0.03634049  
## 6      e2      S   2   2 0.3984917 0.03755572  
## 7      e3      S   3   3 0.4091477 0.03661980  
## 8      e4      S   4   4 0.5404280 0.04824857  
## 9      e5      S   5   5 0.4808798 0.04721738  
## 10     e6      S   6   6 0.5570633 0.05106953  
## 11  varF1      S   7   7 0.6604139 0.06622390  
## 12    cov      S   7   8 0.2951547 0.03933735  
## 13  varF2      S   8   8 0.4505102 0.06260275  
## 14 meanx1      M   1  x1 2.9879512 0.04461557  
## 15 meanx2      M   1  x2 3.0113433 0.04522939  
## 16 meanx3      M   1  x3 2.9861141 0.04431702  
## 17 meanx4      M   1  y1 2.9554046 0.04451830  
## 18 meanx5      M   1  y2 2.9561472 0.04419374  
## 19 meanx6      M   1  y3 2.9673326 0.04597058  
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:             19                   2981              7548.105
##    Saturated:             27                   2973                    NA
## Independence:             12                   2988                    NA
## Number of observations/statistics: 500/3000
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       1586.105               7586.105                 7587.688
## BIC:     -10977.642               7666.182                 7605.875
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2022-02-21 08:50:03 
## Wall clock time: 0.1143224 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.20.6 
## Need help?  See help(mxSummary)

If more detailed output is desired, accessing the $output field on an MxModel object will provide more detailed information including the fit matricies in raw form, the derived Hessian matrix, and data about execution environment and timing.