#
#
#
# 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://www.unt.edu/rss/class/Jon/R_SC/Module8/BiFactor_Figure.pdf

n  <- 5000
f1 <- rnorm(n)
f2 <- rnorm(n)
f3 <- rnorm(n)
g  <- 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)

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



# Minor Updates, Nov. 2011