# How to do Simple Confirmatory Factory Analysis/SEM in R?

I have a data set consisting of 166 observations on 24 variables (a1, a2, a3, a4, b1, ..., b4, ..., f4). Responses where made on a 6-point scale from 1 to 6:

d <- read.table("http://pastebin.com/raw.php?i=m1ZJuKLH")
str(d)
## 'data.frame':   166 obs. of  24 variables:
##  $a1: int 7 7 7 1 1 7 7 7 7 1 ... ##$ a2: int  7 7 4 7 5 5 1 7 7 4 ...
##  $a3: int 7 7 5 1 1 7 7 3 7 1 ... ##$ a4: int  7 7 6 7 1 7 1 7 7 1 ...
##  $b1: int 1 2 5 7 1 4 1 7 4 2 ... ## [...] ##$ f4: int  6 6 4 1 7 7 7 7 7 1 ...


Variables with a common letter (i.e., all as, all bs, ..., all fs) are of the same type.

I would like to test whether the latent factors behind the variable types are correlated with each other (or whether you can set the correlation to 0). An additional complication is that the latent variables c to f share a further latent variable (termed c.to.f) and I only want to test the correlation between a, b and c.to.f.

I would like to run a confirmatory factor analysis (which essentially is a structural equation model) in R testing this. There are at least two mature packages of doing so sem and openMX.
I am interested in opinions/code on which package would be the best or perhaps easiest to specify such a model.

Edit: I would like to accept an answer which includes code samples.

• Have you tried John Fox's sem package? It's capable of doing a fairly stripped down model. Alternately, there's the OpenMX package which is considerably more involved to learn. – AdamO Jan 3 '14 at 20:38
• @Henrik try lavaan github.com/yrosseel/lavaan which is currently the most comprehensive SEM package for R. It allows you to fix/equate coefficients so you can also incorporate the higher order factor. The syntax is also very intuitive. – Momo Jan 3 '14 at 23:56

A CFA is pretty easy to do in R with OpenMx, sem, or lavaan. Since a CFA is such a vanilla case of SEM, all three are pretty easy to implement and offer helpful walkthroughs within their respective documentations. I personally use OpenMx or lavaan. One thing to keep in mind if you use OpenMx is that it won't give you fit statistics by default, you have to specify a saturated model first (or use the semTools package to do this for you).

Because OpenMx hasn't been updated for R version 3 yet (unless you compile from source), here's an example taken from the lavaan walkthrough. It is a CFA with 3 latent variables with three indicators, with covariances among all three latents. More information on the dataset used can be found in the link above.

# load the lavaan package
require(lavaan)

# specify the model
HS.model <- " visual  =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9 "

# fit a full CFA model
fit <- cfa(HS.model, data = HolzingerSwineford1939)

# fit an orthogonal CFA model
fitOrth <- cfa(HS.model, data = HolzingerSwineford1939, orthogonal = TRUE)

# Likelihood ratio test between full and orthogonal model
anova(fit, fitOrth)

# display summary output for full model
summary(fit, fit.measures=TRUE)


Here we see that the orthogonal model (all three covariances set to zero) fits significantly worse than a full CFA. Two things to keep in mind with this code:

1) In this specification, loadings for x1, x4, x7 are fixed to 1 by default to set the scale of the CFA. This can be changed by moving the variables around.

2) Again by default, residual variances are added automatically. This can be changed by adding residual regression weights in the model syntax.