5 Using the OpenMx Package for CFA

The OpenMx package in R is a port of the well-respected MX analytical software. It handles SEM and can easily be used for CFA. Here, the same two models that were run in lavaan will be run again, but an additional model will be run first.

library(OpenMx)

5.1 First OpenMx Model - Single Factor

First, fit a bad model that posits only one underlying factor

5.1.1 set new dataframe and set up basics

Some initial setups

# grab data
#wisc2 <- wisc1[2:12]   # remove ID from df
# set up basics
manifests <- names(wisc2)
observedCov <- cov(wisc2)
numSubjects <- nrow(wisc2)

5.1.2 Create the model using the mxModel function

This code initializes the model and sets vars/paths

latents="F1"
cfa2a <- OpenMx::mxModel("Common Factor Model",type="RAM", 
                 manifestVars = manifests, latentVars = latents,
                 # Now set the residual variance for manifest variables
                 mxPath(from=manifests, arrows=2,free=T,values=1,labels=paste("error",1:11,sep="")),
                 # set latent factor variance to 1
                 mxPath(from="F1",arrows=2,free=F,values=1,labels="varF1"),
                 # specify factor loadings
                 mxPath(from="F1",to=manifests, arrows=1,
                        free=T,values=1,labels=paste("i", 1:11,sep="")),
                 # specify the covariance matrix
                 mxData(observed=observedCov,type="cov",numObs=numSubjects)
) # close the model

5.1.3 Now use mxRun to fit the model

mxRun uses the model defined above.

cfa2a <- OpenMx::mxRun(cfa2a)

Running Common Factor Model with 22 parameters

Summarize the fit:

summary(cfa2a)

Summary of Common Factor Model 
 
free parameters:
      name matrix      row      col  Estimate Std.Error A
1       i1      A     info       F1 2.1084282 0.2049277  
2       i2      A     comp       F1 2.1025863 0.2088113  
3       i3      A    arith       F1 1.2725315 0.1731234  
4       i4      A    simil       F1 2.2597773 0.2232870  
5       i5      A    vocab       F1 2.1749799 0.2034502  
6       i6      A    digit       F1 1.0089237 0.2130399  
7       i7      A pictcomp       F1 1.3505514 0.2286743  
8       i8      A   parang       F1 0.8983356 0.2122489  
9       i9      A    block       F1 1.2315969 0.2112892  
10     i10      A   object       F1 1.0178377 0.2271327  
11     i11      A   coding       F1 0.2012793 0.2347101  
12  error1      S     info     info 3.9873781 0.5453011  
13  error2      S     comp     comp 4.3219872 0.5694074  
14  error3      S    arith    arith 3.6720914 0.4269853  
15  error4      S    simil    simil 4.9709933 0.6497852  
16  error5      S    vocab    vocab 3.8211750 0.5301597  
17  error6      S    digit    digit 6.2528075 0.6886600  
18  error7      S pictcomp pictcomp 6.7364634 0.7623141  
19  error8      S   parang   parang 6.2264645 0.6827525  
20  error9      S    block    block 5.7843994 0.6529180  
21 error10      S   object   object 7.0059957 0.7723127  
22 error11      S   coding   coding 8.1614176 0.8732111  

Model Statistics: 
               |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
       Model:             22                     44              5491.972
   Saturated:             66                      0              5375.122
Independence:             11                     55              5894.326
Number of observations/statistics: 175/66

chi-square:  χ² ( df=44 ) = 116.8505,  p = 1.599332e-08
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       28.85054               160.8505                 167.5084
BIC:     -110.40004               230.4758                 160.8087
CFI: 0.8430635 
TLI: 0.8038294   (also known as NNFI) 
RMSEA:  0.09726823  [95% CI (0.07176259, 0.122825)]
Prob(RMSEA <= 0.05): 0.0002660291
timestamp: 2025-06-25 19:06:04 
Wall clock time: 0.06910896 secs 
optimizer:  SLSQP 
OpenMx version number: 2.21.11 
Need help?  See help(mxSummary) 

What is in the model?

#slotNames(cfa2a@output)
names(cfa2a@output)

 [1] "matrices"                "algebras"               
 [3] "data"                    "SaturatedLikelihood"    
 [5] "IndependenceLikelihood"  "calculatedHessian"      
 [7] "vcov"                    "standardErrors"         
 [9] "gradient"                "hessian"                
[11] "expectations"            "fit"                    
[13] "fitUnits"                "Minus2LogLikelihood"    
[15] "maxRelativeOrdinalError" "minimum"                
[17] "estimate"                "infoDefinite"           
[19] "conditionNumber"         "status"                 
[21] "iterations"              "evaluations"            
[23] "mxVersion"               "frontendTime"           
[25] "backendTime"             "independentTime"        
[27] "wallTime"                "timestamp"              
[29] "cpuTime"                

5.2 Second OpenMx Model - the two factor model

Now fit a model that is the same as the initial model fit with lavaan in chapter 3. Two factors, verbal and performance are established with each manifest variable uniquely specified by only one of the two latent factors.

5.2.1 set new dataframe and set up basics

Same initial setups as above.

# grab data
wisc2 <- wisc1[2:12]   # remove ID from df
# set up basics
manifests <- names(wisc2)
verbalVars <- names(wisc2[1:6])
perfVars <- names(wisc2[7:11])
latents2 <- c("verbal","perf")
#manifests <- c("verbalVars","perfVars")
observedCov <- cov(wisc2[,1:11])
numSubjects <- nrow(wisc2)

5.2.2 Create the model using the mxModel function

Essentially initializes the model and sets vars/paths

cfa2b <- OpenMx::mxModel("Two Factor Model",type="RAM", 
                 manifestVars = manifests, latentVars = latents2,
                 # Now set the residual variance for manifest variables
                 mxPath(from = manifests, arrows=2,free=T,values=1, labels=paste("e",1:11,sep="")),
                 # set latent factor variances
                 mxPath(from=latents2,arrows=2,free=F,values=1,labels=c("p1","p2")),
                 # speciffy factor loadings
                 # Allow the factors to covary as per the lavaan model
                 mxPath(from="verbal", to="perf", arrows=2, 
                        free=T, values=1, labels="latcov1"),
                 # Specify the latent factor paths to manifests
                 mxPath(from="verbal", to = verbalVars,
                        arrows=1,
                        free=T,values=1),
                 mxPath(from="perf", to = perfVars,
                        arrows=1,
                        free=T,values=1),
                 mxData(observed=observedCov,type="cov",numObs=numSubjects)
) # close the model

5.2.3 Now use ’mxRun’to fit the model

mxRun uses the model defined above

cfa2b <- OpenMx::mxRun(cfa2b)

Running Two Factor Model with 23 parameters

Summarize the fit

summary(cfa2b)

Summary of Two Factor Model 
 
free parameters:
                        name matrix      row      col  Estimate  Std.Error A
1   Two Factor Model.A[1,12]      A     info   verbal 2.2061543 0.20134270  
2   Two Factor Model.A[2,12]      A     comp   verbal 2.0421965 0.21197225  
3   Two Factor Model.A[3,12]      A    arith   verbal 1.2998415 0.17269068  
4   Two Factor Model.A[4,12]      A    simil   verbal 2.2320429 0.22525465  
5   Two Factor Model.A[5,12]      A    vocab   verbal 2.2504518 0.20064799  
6   Two Factor Model.A[6,12]      A    digit   verbal 1.0527695 0.21230880 !
7   Two Factor Model.A[7,13]      A pictcomp     perf 1.7420017 0.24615917  
8   Two Factor Model.A[8,13]      A   parang     perf 1.2532318 0.22470470  
9   Two Factor Model.A[9,13]      A    block     perf 1.8459345 0.22453720  
10 Two Factor Model.A[10,13]      A   object     perf 1.6045702 0.23603617  
11 Two Factor Model.A[11,13]      A   coding     perf 0.2070222 0.25706998 !
12                        e1      S     info     info 3.5657068 0.51671427  
13                        e2      S     comp     comp 4.5722769 0.59475403  
14                        e3      S    arith    arith 3.6018519 0.42197020 !
15                        e4      S    simil    simil 5.0955466 0.66622728  
16                        e5      S    vocab    vocab 3.4871862 0.50752196 !
17                        e6      S    digit    digit 6.1624203 0.68092737  
18                        e7      S pictcomp pictcomp 5.5258964 0.77199802  
19                        e8      S   parang   parang 5.4628807 0.65874795  
20                        e9      S    block    block 3.8937546 0.65151342  
21                       e10      S   object   object 5.4673459 0.71961473  
22                       e11      S   coding   coding 8.1590901 0.87416703  
23                   latcov1      S   verbal     perf 0.5888370 0.07718452  

Model Statistics: 
               |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
       Model:             23                     43              5445.730
   Saturated:             66                      0              5375.122
Independence:             11                     55              5894.326
Number of observations/statistics: 175/66

chi-square:  χ² ( df=43 ) = 70.60803,  p = 0.005008876
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:      -15.39197               116.6080                 123.9193
BIC:     -151.47777               189.3981                 116.5642
CFI: 0.9405261 
TLI: 0.9239287   (also known as NNFI) 
RMSEA:  0.06057096  [95% CI (0.02657469, 0.08959509)]
Prob(RMSEA <= 0.05): 0.2334782
timestamp: 2025-06-25 19:06:05 
Wall clock time: 0.09197903 secs 
optimizer:  SLSQP 
OpenMx version number: 2.21.11 
Need help?  See help(mxSummary) 

What is in the model?

slotNames(cfa2b@output)

NULL

names(cfa2b@output)

 [1] "matrices"                "algebras"               
 [3] "data"                    "SaturatedLikelihood"    
 [5] "IndependenceLikelihood"  "calculatedHessian"      
 [7] "vcov"                    "standardErrors"         
 [9] "gradient"                "hessian"                
[11] "expectations"            "fit"                    
[13] "fitUnits"                "Minus2LogLikelihood"    
[15] "maxRelativeOrdinalError" "minimum"                
[17] "estimate"                "infoDefinite"           
[19] "conditionNumber"         "status"                 
[21] "iterations"              "evaluations"            
[23] "mxVersion"               "frontendTime"           
[25] "backendTime"             "independentTime"        
[27] "wallTime"                "timestamp"              
[29] "cpuTime"                

5.2.4 Can we use semPlot to draw the OpenMx model fit?

The semPlot package is very powerful and can recognize many lm and sem model objects. We can use the identical code that we used in chapter 2 for the lavaan model.

# Note that the base plot, including standardized path coefficients plots positive coefficients green
# and negative coefficients red.  Red-green colorblindness issues anyone?
# I redrew it here to choose a blue and red.  But all the coefficients in this example are
# positive,so they are shown with the skyblue.
# more challenging to use colors other than red and green.  not in this doc
semPlot::semPaths(cfa2b, residuals=F,sizeMan=7,"std",
         posCol=c("skyblue4", "red"),
         #edge.color="skyblue4",
         edge.label.cex=1.2,layout="circle2")

# or we could draw the paths in such  a way to include the residuals:
#semPaths(fit1, sizeMan=7,"std",edge.color="skyblue4",edge.label.cex=1,layout="circle2")
# the base path diagram can be drawn much more simply:
#semPaths(fit1)
# or
semPlot::semPaths(cfa2b,"std")

5.3 Third OpenMx Model

Now fit a model that is the same as the second model fit with lavaan above. Two factors, verbal and performance, plus paths from both latents to “comp”.

5.3.1 set new dataframe and set up basics

Same initial setups as above.

# grab data
wisc2 <- wisc1[2:12]   # remove ID from df
# set up basics
manifests <- names(wisc2)
verbalVars <- names(wisc2[1:6])
perfVars <- c(names(wisc2[7:11]),"comp")
latents2 <- c("verbal","perf")
#manifests <- c("verbalVars","perfVars")
observedCov <- cov(wisc2[,1:11])
numSubjects <- nrow(wisc2)

5.3.2 Create the model using the mxModel function

Essentially initializes the model and sets vars/paths

cfa2c <- OpenMx::mxModel("TwoFac, Comp Common",type="RAM", 
                 manifestVars = manifests, latentVars = latents2,
                 # Now set the residual variance for manifest variables
                 mxPath(from = manifests, arrows=2,free=T,values=1, labels=paste("e",1:11,sep="")),
                 # set latent factor variances
                 mxPath(from=latents2,arrows=2,free=F,values=1,labels=c("p1","p2")),
                 # speciffy factor loadings
                 # Allow the factors to covary as per the lavaan model
                 mxPath(from="verbal", to="perf", arrows=2, 
                        free=T, values=1, labels="latcov1"),
                 # Specify the latent factor paths to manifests
                 mxPath(from="verbal", to = verbalVars,
                        arrows=1,
                        free=T,values=1),
                 mxPath(from="perf", to = perfVars,
                        arrows=1,
                        free=T,values=1),
                 mxData(observed=observedCov,type="cov",numObs=numSubjects)
) # close the model

5.3.3 Now use ‘mxRun’ to fit the model

mxRun uses the model defined above

cfa2c <- OpenMx::mxRun(cfa2c)

Running TwoFac, Comp Common with 24 parameters

Summarize the fit

summary(cfa2c)

Summary of TwoFac, Comp Common 
 
free parameters:
                           name matrix      row      col  Estimate  Std.Error A
1   TwoFac, Comp Common.A[1,12]      A     info   verbal 2.2560643 0.20030711  
2   TwoFac, Comp Common.A[2,12]      A     comp   verbal 1.4912880 0.25643470 !
3   TwoFac, Comp Common.A[3,12]      A    arith   verbal 1.3067883 0.17304451  
4   TwoFac, Comp Common.A[4,12]      A    simil   verbal 2.2047514 0.22731986  
5   TwoFac, Comp Common.A[5,12]      A    vocab   verbal 2.2734881 0.20058708  
6   TwoFac, Comp Common.A[6,12]      A    digit   verbal 1.0747643 0.21233735  
7   TwoFac, Comp Common.A[2,13]      A     comp     perf 0.8841463 0.27029754  
8   TwoFac, Comp Common.A[7,13]      A pictcomp     perf 1.7896201 0.24239834 !
9   TwoFac, Comp Common.A[8,13]      A   parang     perf 1.1888316 0.22509928 !
10  TwoFac, Comp Common.A[9,13]      A    block     perf 1.8227666 0.22163101 !
11 TwoFac, Comp Common.A[10,13]      A   object     perf 1.6328343 0.23323255 !
12 TwoFac, Comp Common.A[11,13]      A   coding     perf 0.2004685 0.25520584  
13                           e1      S     info     info 3.3430411 0.50942368  
14                           e2      S     comp     comp 4.3310604 0.55728973 !
15                           e3      S    arith    arith 3.5837408 0.42201358  
16                           e4      S    simil    simil 5.2166719 0.68262345  
17                           e5      S    vocab    vocab 3.3829818 0.50743874  
18                           e6      S    digit    digit 6.1156232 0.67754036  
19                           e7      S pictcomp pictcomp 5.3576276 0.75572310 !
20                           e8      S   parang   parang 5.6202038 0.66563236 !
21                           e9      S    block    block 3.9787852 0.63669887  
22                          e10      S   object   object 5.3759025 0.70817556 !
23                          e11      S   coding   coding 8.1617269 0.87416234  
24                      latcov1      S   verbal     perf 0.5332247 0.08203393 !

Model Statistics: 
               |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
       Model:             24                     42              5435.732
   Saturated:             66                      0              5375.122
Independence:             11                     55              5894.326
Number of observations/statistics: 175/66

chi-square:  χ² ( df=42 ) = 60.61003,  p = 0.03139638
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:      -23.38997               108.6100                 116.6100
BIC:     -156.31098               184.5649                 108.5643
CFI: 0.9599098 
TLI: 0.9475009   (also known as NNFI) 
RMSEA:  0.05031876  [95% CI (0, 0.08138947)]
Prob(RMSEA <= 0.05): 0.4663969
timestamp: 2025-06-25 19:06:11 
Wall clock time: 0.07219601 secs 
optimizer:  SLSQP 
OpenMx version number: 2.21.11 
Need help?  See help(mxSummary) 

What is in the model?

slotNames(cfa2c@output)

NULL

names(cfa2c@output)

 [1] "matrices"                "algebras"               
 [3] "data"                    "SaturatedLikelihood"    
 [5] "IndependenceLikelihood"  "calculatedHessian"      
 [7] "vcov"                    "standardErrors"         
 [9] "gradient"                "hessian"                
[11] "expectations"            "fit"                    
[13] "fitUnits"                "Minus2LogLikelihood"    
[15] "maxRelativeOrdinalError" "minimum"                
[17] "estimate"                "infoDefinite"           
[19] "conditionNumber"         "status"                 
[21] "iterations"              "evaluations"            
[23] "mxVersion"               "frontendTime"           
[25] "backendTime"             "independentTime"        
[27] "wallTime"                "timestamp"              
[29] "cpuTime"                

5.3.4 Draw the OpenMx model number 3

The semPlot package is very powerful and can recognize many ‘lm’ and SEM model objects. We can use the identical code that we used in chapter 2 for the lavaan model.

# Note that the base plot, including standardized path coefficients plots positive coefficients green
# and negative coefficients red.  Red-green colorblindness issues anyone?
# I redrew it here to choose a blue and red.  But all the coefficients in this example are
# positive,so they are shown with the skyblue.
# more challenging to use colors other than red and green.  not in this doc
semPlot::semPaths(cfa2c, residuals=F,sizeMan=7,"std",
         posCol=c("skyblue4", "red"),
         #edge.color="skyblue4",
         edge.label.cex=1.2,layout="circle2")

# or we could draw the paths in such  a way to include the residuals:
#semPaths(fit1, sizeMan=7,"std",edge.color="skyblue4",edge.label.cex=1,layout="circle2")
# the base path diagram can be drawn much more simply:
#semPaths(fit1)
# or
semPlot::semPaths(cfa2c,"std")

5.4 Compare the OpenMx models

In OpenMx a convenient function exists for model comparisons. This code compares the same two models that we initially compared in the lavaan approach that used the anova function.

mxCompare(cfa2c,cfa2b)

                 base       comparison ep minus2LL df      AIC   diffLL diffdf
1 TwoFac, Comp Common             <NA> 24 5435.732 42 5483.732       NA     NA
2 TwoFac, Comp Common Two Factor Model 23 5445.730 43 5491.730 9.997993      1
            p
1          NA
2 0.001567109

Allaire, J., Xie, Y., McPherson, J., Luraschi, J., Ushey, K., Atkins, A., … Iannone, R. (2018). Rmarkdown: Dynamic documents for r. Retrieved from https://CRAN.R-project.org/package=rmarkdown

Boker, S. M., Neale, M. C., Maes, H. H., Spiegel, M., Brick, T. R., Estabrook, R., … Kirkpatrick, R. M. (2019). OpenMx: Extended structural equation modelling. Retrieved from https://CRAN.R-project.org/package=OpenMx

Epskamp, S., & Simon Stuber, with contributions from. (2017). semPlot: Path diagrams and visual analysis of various SEM packages’ output. Retrieved from https://CRAN.R-project.org/package=semPlot

Fox, J., Nie, Z., & Byrnes, J. (2017). Sem: Structural equation models. Retrieved from https://CRAN.R-project.org/package=sem

Fox, J., Weisberg, S., & Price, B. (2018). Car: Companion to applied regression. Retrieved from https://CRAN.R-project.org/package=car

Korkmaz, S., Goksuluk, D., & Zararsiz, G. (2018). MVN: Multivariate normality tests. Retrieved from https://CRAN.R-project.org/package=MVN

Revelle, W. (2019). Psych: Procedures for psychological, psychometric, and personality research. Retrieved from https://CRAN.R-project.org/package=psych

Rosseel, Y. (2018). Lavaan: Latent variable analysis. Retrieved from https://CRAN.R-project.org/package=lavaan

Xie, Y. (2018). Knitr: A general-purpose package for dynamic report generation in r. Retrieved from https://CRAN.R-project.org/package=knitr