bifactor_lavaan.R

# script stolen and adapted from http://goo.gl/kUzaov


if (!require('sem')) install.packages('sem')
if (!require('lavaan')) install.packages('lavaan')
if (!require('Rcmdr')) install.packages('Rcmdr')


# creating a dataset for g + 3 factors, all orthogonal

n  <- 5000
g  <- rnorm(n)
f1 <- rnorm(n)
f2 <- rnorm(n)
f3 <- rnorm(n)

x1 <- .9*f1 + .7*g + rnorm(n,0,sqrt(1 - (.9*.7)^2))
x2 <- .8*f1 + .7*g + rnorm(n,0,sqrt(1 - (.8*.7)^2))
x3 <- .7*f1 + .7*g + rnorm(n,0,sqrt(1 - (.7*.7)^2))
x4 <- .6*f1 + .7*g + rnorm(n,0,sqrt(1 - (.6*.7)^2))
x5 <- .6*f2 + .7*g + rnorm(n,0,sqrt(1 - (.6*.7)^2))
x6 <- .7*f2 + .7*g + rnorm(n,0,sqrt(1 - (.7*.7)^2))
x7 <- .8*f2 + .7*g + rnorm(n,0,sqrt(1 - (.8*.7)^2))
x8 <- .9*f2 + .7*g + rnorm(n,0,sqrt(1 - (.9*.7)^2))
x9 <- .6*f3 + .7*g + rnorm(n,0,sqrt(1 - (.6*.7)^2))
x10 <- .7*f3 + .7*g + rnorm(n,0,sqrt(1 - (.7*.7)^2))
x11 <- .8*f3 + .7*g + rnorm(n,0,sqrt(1 - (.8*.7)^2))
x12 <- .9*f3 + .7*g + rnorm(n,0,sqrt(1 - (.9*.7)^2))

bifac <- data.frame(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12)
rm(n,f1,f2,f3,g,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12)
summary(bifac)
nrow(bifac)

numSummary(bifac, statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.25,.5,.75,1))

# Variance of each variable listed on the diagonal of the covariance matrix.

cov(bifac)

############
#Bifactor with SEM package
############

library(sem)

# Put the covariance matrix into an object. 

cov.sem <- cov(bifac)

## Specify the Measurement Model.

m.model.1 <- specifyModel()
general.factor -> x1, loadg1, NA
general.factor -> x2, loadg2, NA
general.factor -> x3, loadg3, NA
general.factor -> x4, loadg4, NA
general.factor -> x5, loadg5, NA
general.factor -> x6, loadg6, NA
general.factor -> x7, loadg7, NA
general.factor -> x8, loadg8, NA
general.factor -> x9, loadg9, NA
general.factor -> x10, loadg10, NA
general.factor -> x11, loadg11, NA
general.factor -> x12, loadg12, NA
factor1 -> x1, loadf11, NA
factor1 -> x2, loadf12, NA
factor1 -> x3, loadf13, NA
factor1 -> x4, loadf14, NA
factor2 -> x5, loadf25, NA
factor2 -> x6, loadf26, NA
factor2 -> x7, loadf27, NA
factor2 -> x8, loadf28, NA
factor3 -> x9, loadf39, NA
factor3 -> x10, loadf310, NA
factor3 -> x11, loadf311, NA
factor3 -> x12, loadf312, NA
general.factor <-> factor1, NA, 0
general.factor <-> factor2, NA, 0
general.factor <-> factor3, NA, 0
factor1 <-> factor2, NA, 0
factor1 <-> factor3, NA, 0
factor2 <-> factor3, NA, 0
x1 <-> x1, evar1, NA
x2 <-> x2, evar2, NA
x3 <-> x3, evar3, NA
x4 <-> x4, evar4, NA
x5 <-> x5, evar5, NA
x6 <-> x6, evar6, NA
x7 <-> x7, evar7, NA
x8 <-> x8, evar8, NA
x9 <-> x9, evar9, NA
x10 <-> x10, evar10, NA
x11 <-> x11, evar11, NA
x12 <-> x12, evar12, NA
general.factor <-> general.factor, NA, 1
factor1 <-> factor1, NA, 1
factor2 <-> factor2, NA, 1
factor3 <-> factor3, NA, 1

## Run the SEM (on the covariance matrix object). Notice in the output all 
# the loadings for the General Factor are very close to 0.7 (as specified 
# when generating the data). Likewise, all the sub-factor loadings very 
# closely correspond to the 0.9, 0.8, 0.7, 0.6 which was specified when 
# generating the data above.

m.modelfit.1 <- sem(m.model.1, cov.sem, 5000, maxiter = 10000)
summary(m.modelfit.1, conf.level=0.95)

#########
#Bifactor with lavaan Package
#########

# Fit the Confirmatory Factor Analysis (CFA).

m.model.2 <- '
general.factor  =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12
factor1         =~ x1 + x2 + x3 + x4
factor2         =~ x5 + x6 + x7 + x8
factor3         =~ x9 + x10 + x11 + x12
general.factor  ~~ 0*factor1
general.factor  ~~ 0*factor2
general.factor  ~~ 0*factor3
factor1         ~~ 0*factor2
factor1         ~~ 0*factor3
factor2         ~~ 0*factor3
'

 
#  

m.modelfit.2 <- cfa(m.model.2, data = bifac, std.lv = TRUE, information="observed") # std.lv - if TRUE, the metric of each latent variable is determined by fixing their variances to 1.0. If FALSE, the metric of each latent variable is determined by fixing the factor loading of the first indicator to 1.0.
summary(m.modelfit.2, fit.measures = FALSE, standardized = FALSE) fit.measures = TRUE command provides additional fit indices: CFI, TLI, AIC, BIC, sample-size adjusted BIC, RMSEA # standardized = TRUE command provides a summary with standardized loadings as well as non-standardized loadings
summary(m.modelfit.2, fit.measures = TRUE, standardized = FALSE)
inspect(m.modelfit.2)