bifactor model using lavaan

bifactor_lavaan.R

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


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)

cor(bifac)


# 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)

If the packages aren't installed you install them at the top of the script but never load them once they're installed.

In this sample, you have to insert 'librasry(xxx) inside every if, to permit package loading