#
#
#
# Estimation of a Bi-Factor Model in packages 'sem' & 'lavaan' using simulated data.
#
#
# This script assumes you have worked through all the previous notes from
# the web page and you have downloaded, installed, and updated all available
# R packages.
# Load the following libraries if you have not already; when an additional
# library is needed, it will be listed in the script and loaded.
library(Rcmdr)
# Generate the data. This data is designed to produce one General factor (g) and
# three Sub-factors (f1, f2, & f3). All four latent factors are orthogonal.
# To see a graphical representation of this bi-factor model, paste the following
# URL into your browser and hit the enter key:
# http://researchsupport.unt.edu/class/Jon/R_SC/Module8/BiFactor_Figure.pdf
n <- 5000
library(psych)
library(MASS)
library(GPArotation)
group.factor.loads <- matrix(c(.6,.8,.4),nrow=3)
item.factor.loads <- matrix(c(
.8,0,0,
.7,0,0,
.6,0,0,
.6,0,0,
0,.7,0,
0,.6,0,
0,.5,0,
0,.5,0,
0,0,.6,
0,0,.5,
0,0,.5,
0,0,.4),
ncol=3,byrow=TRUE)
data.df <- sim.hierarchical(gload = group.factor.loads,
fload = item.factor.loads,
n = 1000, raw = TRUE, mu = c(rep(10,12)))
bifac <- data.frame(data.df$observed); rm(group.factor.loads, item.factor.loads, data.df)
rm(n)
names(bifac) <- c("x1","x2","x3","x4","x5","x6",
"x7","x8","x9","x10","x11","x12")
detach("package:psych")
detach("package:GPArotation")
detach("package:MASS")
summary(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)
####### Using the 'sem' package to estimate the bi-factor model.
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)
detach("package:sem")
####### Using the 'lavaan' package to estimate the bi-factor model.
library(lavaan)
##### Fit the Confirmatory Factor Analysis (CFA).
### Specify the Measurement Model.
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
'
### Fit the measurement model.
# * The "std.lv = TRUE" command sets the variance of all the latent factors to 1 (by default, the first loading
# of each latent variable is set to one).
# ** The "fit.measures = TRUE" command provides additional fit indices (only chi-square is given by default).
# *** The "standardize = TRUE" command provides standardized loadings (non-standardized loadings are given
# by default).
m.modelfit.2 <- cfa(m.model.2, data = bifac, std.lv = TRUE, information="observed")
summary(m.modelfit.2, fit.measures = TRUE, standardize = FALSE)
inspect(m.modelfit.2)
ls()
# Updated, Feb. 2015