#
#
###### Familiarization with SEM using the 'lavaan' package. ######
#
#
# 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 library which is necessary to import an SPSS data file.

library(foreign)

### Upload the SAS SEM Example data ("SEMData.sav") which is SIMULATION DATA.
# This data can be imported directly from the web.

exsem <- read.spss("http://www.unt.edu/rss/class/Jon/R_SC/Module8/SAS_Ex/SEMData.sav", use.value.labels=TRUE, 
  max.value.labels=Inf, to.data.frame=TRUE)
detach("package:foreign")
head(exsem)

attach(exsem)

####### Lavaan #######
# Lavaan is a relatively new package (May 2010) which was created to make it easier for new R users 
# to conduct latent variable modeling (e.g. confirmatory factor analysis, structural equation modeling, 
# and latent growth curve models). The lavaan package requires at least 5 other packages: MASS, boot,
# mnormt, pbivnorm, & quadprog

library(lavaan)

# The primary benefit of the Lavaan package is an intuitive set of operators and structure for specifying elements
# of a model. As an example, consider the following generic (not used here) example from the Lavaan documentation.
# NOT RUN; simply shown as generic examples.

myModel <- '
            # Regression or path equations
              y  ~ F1 + F2 + x1 + x2
              F1 ~ F2 + F3
              F2 ~ F3 + x1 + x2
            # Latent variable definitions
              F1 =~ y1 + y2 + y3
              F2 =~ y4 + y5 + y6
              F3 =~ y7 + y8 + y9 + y10
            # Variances and covariances
              y1 ~~ y1
              y1 ~~ y2
              F1 ~~ F2
            # Intercepts
              y1 ~ 1
              F1 ~ 1
           '

# The developers of the Lavaan package wanted to provide an easy to use and free alternative to 
# commercially available modeling software. The developers state that "results of the lavaan 
# package are typically very close, if not identical, to the results of the comercial package 
# Mplus" (Rosseel, 2010) Lavaan contains an arguement for comparing its results directly to 
# other modeling software. The optional arguement is 'mimic.Mplus=FALSE'. If 'TRUE', an 
# attempt is made to mimic Mplus behavior as close as possible. If 'FALSE', results are more 
# in line with EQS or LISREL.

# Rosseel, Y. (2010). lavaan: An R package for structural equation modeling and more. 
# Available at: http://cran.r-project.org/web/packages/lavaan/index.html

# See also: http://lavaan.ugent.be/

#################################################################################################
############### CURRENT EXAMPLES BELOW ###############

##### Fit the Confirmatory Factor Analysis (CFA).

### Specify the Measurement Model.

m.model.1 <- '
Personality  =~ extro + open + agree
Engagement   =~ social + cognitive + physical + cultural
Crystallized =~ vocab + abstruse + block
Fluid        =~ common + sets + series
'

### 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.model.fit <- cfa(m.model.1, data = exsem, std.lv = TRUE)
summary(m.model.fit, fit.measures = TRUE, standardize = TRUE)

### Specify the Structural Model. 

# Recall, by default, the lavaan package model specification functions (cfa & sem) constrain the first loading 
# of each latent variable to 1. 

s.model.1 <- '
Personality  =~ extro + open + agree
Engagement   =~ social + cognitive + physical + cultural
Crystallized =~ vocab + abstruse + block
Fluid        =~ common + sets + series
Engagement   ~ Personality
Crystallized ~ Engagement + Personality
Fluid        ~ Engagement + Personality
Crystallized ~~ Fluid
'

### Fit the Structural Model.

s.model.fit <- sem(s.model.1, data = exsem)
summary(s.model.fit, fit.measures = TRUE, standardize = TRUE)

# To review the residuals, use the 'residuals' function.

resids <- residuals(s.model.fit)
resids

hist(resids$cov)

# As a comparison, you can see that these results are virtually identical to 
# those furnished by SAS:
# http://www.unt.edu/rss/class/Jon/SAS_SC/SEM/SEM2/SAS_Module8_SEM2.htm


# Minor updates: Feb. 2011.