Angelos Markos bio photo

Angelos Markos

Data Scientist, Runner, Assistant Professor

Email Twitter Facebook LinkedIn

What is lavaan?

lavaan is a free, open source R package for latent variable analysis. You can use lavaan to estimate a large variety of multivariate statistical models, including path analysis, confirmatory factor analysis, structural equation modeling and growth curve models.

What is “Confirmatory Factor Analysis for Applied Research”?

It is a very comprehensive book by Timothy A. Brown emphasizing practical and theoretical aspects of Confirmatory Factor Analysis (CFA), a well-established multivariate statistics method. CFA/SEM is mostly performed via commercial software, such as Mplus, EQS, Lisrel and Amos. The R package lavaan offers a very promising open source alternative. The companion website of Brown’s book on CFA currently contains the syntax code for running the book examples in most (if not all) of the above commercial packages. So what’s missing? A lavaan implementation. Here it is, following the chapters of the book.

Chapter 4 - Specification and Interpretation of CFA Models

Contains two examples a) Two-Factor CFA (Neuroticism, Extraversion) and b) CFA with Single Indicators: Health Status

#############################################################
# Run this code at your own risk! 
############################################################# 
# Chapter 4: Specification and Interpretation of CFA Models

install.packages("lavaan")
library(lavaan)

# Two Factor CFA Model of Neuroticism and Extraversion (Tables 4.2, 4.3, 4.4)

# input: the sample variance-covariance matrix of page 116
covmat <- '32.49
 24.4826 31.36
 26.6669 25.4106 40.96
 25.2772 23.557 27.7978 32.49
 -12.0042 -10.1472 -13.6704 -10.8756 36
 -11.1674 -9.7216 -11.904 -9.4358 25.11 38.44
 -9.617 -9.2249 -10.8346 -9.617 21.6828 23.0063 32.49
 -9.0014 -7.9654 -10.4653 -7.8204 17.9424 20.589 18.0667 31.36'

neodat.cov <- getCov(covmat, names=c("N1","N2","N3","N4","E1","E2","E3","E4"))

#set up model
new.model <- ' neuro =~ N1 + N2 + N3 + N4
extrav =~ E1 + E2 + E3 + E4'

# fit model
fit <- cfa(new.model, sample.cov=neodat.cov, sample.nobs=250, mimic="Mplus")

#summary statistics
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

parameterEstimates(fit)   #coefficients with se, 95%CI and z values etc. 
fitted(fit)               #the model implied variance-covariance matrix, slightly different values from Table 4.2 p. 117
resid(fit) # unstandardized residuals matrix
resid(fit,type="standardized") #standardized residuals matrix, there are differences with the values found in Table 4.2 p. 117

#############################################################
#CFA with Single Indicators: Health Status
#############################################################

# read Health Status data from the book's companion website
health <- read.table("http://people.bu.edu/tabrown/Ch4/fig4.3.dat")
# remove id variable
health <- health[,-1]
# assign variable names
names(health) <- c("ACTIV","SOMA","PAIN","MENTH","SOCF","VITAL","GENHLTH","AGE")
head(health)

#set up model
model <- '
# measurement model
psysfun =~ ACTIV + SOMA + PAIN
mentfun =~ MENTH + SOCF + VITAL
AGEF =~ AGE
GWB =~ GENHLTH
#fix variances of single indicators
AGE ~~ 0*AGE
GENHLTH ~~ 7.88*GENHLTH
# residual correlations
ACTIV ~~ SOMA
'

# fit model
fit <- cfa(model, data=health,mimic="Mplus")
# summary statistics
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

Chapter 5 - CFA Model Revision and Comparison

Contains two examples a) a Three-Factor CFA Model (Alcohol Motives) and b) E/CFA (EFA in the CFA framework)

################################################
# Chapter 5: CFA Model Revision and Comparison #
################################################

library(lavaan)

####################################################################
# Three-Factor CFA (Alcohol Motives QN) model (Figure 1, page 168) #
####################################################################

#input: the lower part of a correlation matrix, see p.169, Table 5.2 (the lisrel syntax)
cormat <- '
1.000
0.300 1.000
0.229 0.261 1.000
0.411 0.406 0.429 1.000
0.172 0.252 0.218 0.481 1.000
0.214 0.268 0.267 0.579 0.484 1.000
0.200 0.214 0.241 0.543 0.426 0.492 1.000
0.185 0.230 0.185 0.545 0.463 0.548 0.522 1.000
0.134 0.146 0.108 0.186 0.122 0.131 0.108 0.151 1.000
0.134 0.099 0.061 0.223 0.133 0.188 0.105 0.170 0.448 1.000
0.160 0.131 0.158 0.161 0.044 0.124 0.066 0.061 0.370 0.350 1.000
0.087 0.088 0.101 0.198 0.077 0.177 0.128 0.112 0.356 0.359 0.507 1.000'

#convert it to a square matrix
Cmat <- getCov(cormat)

# the standard deviations of the variables (Table 5.2)
neodat.sdev = c(2.06, 1.52, 1.92, 1.41, 1.73, 1.77, 2.49, 2.27, 2.68, 1.75, 2.57, 2.66)

# convert the correlation matrix to a covariance matrix (lavaan needs a covariance matrix as input)
Dmat = diag(neodat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10","X11","X12")
colnames(covmat)=rownames(covmat)               

#model set up - this is the model presented in Figure 1, p.168 
# you can play around with this model following the text of Ch.5
new.model <- '#X4 loads on both coping and social factors
              coping =~ X1 + X2 + X3 + X4   
              social =~ X5 + X6 + X7 + X8 + X4
              enhance =~ X9 + X10 + X11 + X12
              # residual correlation
              X11 ~~ X12'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=500)

#results are almost identical due to rounding errors in the matrix conversion step
fit

summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


########################################################
#E/CFA (EFA in the CFA framework) page 193 - Table 5.8 #
########################################################

# read Alcohol Motives data from the book's companion website
alcohol <- read.table("http://people.bu.edu/tabrown/Ch5/efa.dat")
#remove last column
alcohol <- alcohol[,-13]
# assign variable names
names(alcohol) <- c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10","X11","X12")
head(alcohol)

#set up model
new.model <- '#Force the loadings of the first indicators to be free 
              coping =~ NA*X1 + X2+ X3 + X4 + X5 + X6 + X7 + X9 + X10 + X11  
              social =~ NA*X8 + X2 + X3 + X4 + X5 + X6 + X7 + X9 + X10 + X11 
              enhance =~ NA*X12 + X2 + X3 + X4 + X5 + X6 + X7 + X9 + X10 + X11 
              #factor variances fixed to 1.0
              coping ~~ 1*coping
              social ~~ 1*social
              enhance ~~ 1*enhance'

#fit model
fit <- sem(new.model, data=alcohol)

summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

Chapter 6 - CFA of Multitrait-Multimethod Matrices

Contains two examples a) a Correlated Methods Model of a Multitrait-Multimethod Matrix (MTMM) and b) a Correlated Uniqueness Model of a MTMM matrix

#####################################################
# Chapter 6: CFA of Multitrait-Multimethod Matrices #
#####################################################

library(lavaan)

######################################################
# Correlated Methods Model of the MTMM Matrix, p.218 #
######################################################

#Input matrix: the correlation matrix of Table 6.2 (Lisrel syntax)
cormat <- '
1.000
0.290  1.000
0.372  0.478  1.000
0.587  0.238  0.209  1.000
0.201  0.586  0.126  0.213  1.000
0.218  0.281  0.681  0.195  0.096  1.000
0.557  0.228  0.195  0.664  0.242  0.232  1.000
0.196  0.644  0.146  0.261  0.641  0.248  0.383  1.000
0.219  0.241  0.676  0.290  0.168  0.749  0.361  0.342  1.000'


Cmat <- getCov(cormat)

neodat.sdev = c(3.61, 3.66, 3.59, 2.94, 3.03, 2.85, 2.22, 2.42, 2.04)
Dmat = diag(neodat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("PARI","SZTI","SZDI","PARC","SZTC","SZDC","PARO","SZTO","SZDO")
colnames(covmat)= rownames(covmat)

#Correlated Methods Model p.218
#model set up - this is the model presented in Figure 6.1, p.218 
new.model <- 'paranoid =~ PARI +  PARC + PARO   
              schizotypal =~ SZTI + SZTC + SZTO
              schizoid =~ SZDI + SZDC + SZDO
              inventory =~ SZTI + PARI + SZDI
              clininter =~ PARC + SZTC + SZDC
              obsrating =~ PARO + SZTO + SZDO'
            
#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=500)

#The model has NOT converged! And this is exactly what we expected. 
# Brown does not present the results of this model for the same reason.
fit


#########################################################
# Correlated Uniqueness Model of the MTMM Matrix, p.220 #
#########################################################


#Input matrix: the correlation matrix of Table 6.3 (Lisrel syntax)
cormat <- '
1.000
0.290  1.000
0.372  0.478  1.000
0.587  0.238  0.209  1.000
0.201  0.586  0.126  0.213  1.000
0.218  0.281  0.681  0.195  0.096  1.000
0.557  0.228  0.195  0.664  0.242  0.232  1.000
0.196  0.644  0.146  0.261  0.641  0.248  0.383  1.000
0.219  0.241  0.676  0.290  0.168  0.749  0.361  0.342  1.000'


Cmat <- getCov(cormat)

neodat.sdev = c(3.61, 3.66, 3.59, 2.94, 3.03, 2.85, 2.22, 2.42, 2.04)
Dmat = diag(neodat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("PARI","SZTI","SZDI","PARC","SZTC","SZDC","PARO","SZTO","SZDO")
colnames(covmat)= rownames(covmat)

#Correlated Uniqueness Model of the MTMM Matrix, p.228
#model set up - this is the model presented in Figure 6.2, p.220 
new.model <- 'paranoid =~ PARI +  PARC + PARO   
              schizotypal =~ SZTI + SZTC + SZTO
              schizoid =~ SZDI + SZDC + SZDO
              # residual correlation
              PARI ~~ SZTI
              SZTI ~~ SZDI
              PARI ~~ SZDI
              PARC ~~ SZTC
              SZTC ~~ SZDC
              PARC ~~ SZDC
              PARO ~~ SZTO
              SZTO ~~ SZDO
              PARO ~~ SZDO'
#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=500)

#results same as in Table 6.4, page 226
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

Chapter 7 - CFA with Equality Constraints, Multiple Groups, and Mean Structures

Contains examples for a) Tau Equivalent and Parallel Indicators b) Longitudinal Measurement Invariance c) Multiple Groups CFA (Depression example) d) MIMIC Models (Phobia example)

#####################################################
# Chapter 7: CFA with Equality Constraints,         #
#            Multiple Groups, and Mean Structures   #
#####################################################

library(lavaan)

#####################################
# Equality Constraints p.236 - 252  #
#####################################

#Input matrix: the correlation matrix of Figure 7.1
cormat <- '
1.000
0.661  1.000
0.630  0.643  1.000
0.270  0.300  0.268  1.000
0.297  0.265  0.225  0.805  1.000
0.290  0.287  0.248  0.796  0.779  1.000'

Cmat <- getCov(cormat)

neodat.sdev = c(2.610,  2.660,  2.590,  1.940,  2.030,  2.050)
Dmat = diag(neodat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("x1","x2","x3","x4","x5","x6")
colnames(covmat)= rownames(covmat)

############ Congeneric Model #############
#model set up - this is the model presented in Figure 7.1, p.238 
new.model <- 'auditory =~ x1 +  x2 + x3   
              visual =~ x4 + x5 + x6'            


#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=200)

#results same as in Table 7.1, page 241
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


############# Tau-equivalent Model - Auditory Memory only ############
#model set up - this is the model presented in Figure 7.1, p.238 
new.model <- 'auditory =~ v1*x1 +  v1*x2 + v1*x3
              visual =~ x4 + x5 + x6'            

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=200)

#results same as in Table 7.3, page 246 (only model fit)
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


############# Tau-equivalent Model - Auditory & Visual memory ############
#model set up - this is the model presented in Figure 7.1, p.238 
new.model <- 'auditory =~ v1*x1 +  v1*x2 + v1*x3
              visual =~ v2*x4 + v2*x5 + v2*x6'
              
#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=200)

#results same as in Table 7.3, page 246 (only model fit)
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


############# Parallel Model - Auditory Memory only ############
#model set up - this is the model presented in Figure 7.1, p.238 
new.model <- 'auditory =~ v1*x1 +  v1*x2 + v1*x3
              visual =~ v2*x4 + v2*x5 + v2*x6
              #equality constraints of error variances (Auditory)
              f1*x1 ~~ f1*x1
              f1*x2 ~~ f1*x2
              f1*x3 ~~ f1*x3'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=200)

#results same as in Table 7.3, page 246 (only model fit)
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


############# Parallel Model - Auditory & Visual memory ############
#model set up - this is the model presented in Figure 7.1, p.238 
new.model <- 'auditory =~ v1*x1 +  v1*x2 + v1*x3
              visual =~ v2*x4 + v2*x5 + v2*x6
              #equality constraints of error variances (Auditory & Visual)
              f1*x1 ~~ f1*x1
              f1*x2 ~~ f1*x2
              f1*x3 ~~ f1*x3
              f2*x4 ~~ f2*x4
              f2*x5 ~~ f2*x5
              f2*x6 ~~ f2*x6'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=200)

#results same as in Table 7.4, page 248
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

#to see the internal representation of the model (this helps understand the model definition)
inspect(fit,what="list")

########################################
# Longitudinal Invariance p.252 - 268  #
########################################

#Input matrix: the correlation matrix of Figure 7.2
cormat <- '
1.000
0.736  1.000
0.731  0.648  1.000
0.771  0.694  0.700  1.000
0.685  0.512  0.496  0.508  1.000
0.481  0.638  0.431  0.449  0.726  1.000
0.485  0.442  0.635  0.456  0.743  0.672  1.000
0.508  0.469  0.453  0.627  0.759  0.689  0.695  1.000'

Cmat <- getCov(cormat)

neodat.sdev = c(1.940,  2.030,  2.050,  1.990,  2.610,  2.660, 2.590,  2.550)
Dmat = diag(neodat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("A1","B1","C1","D1","A2","B2","C2","D2")
colnames(covmat)= rownames(covmat)

#Means: 1.500  1.320  1.450  1.410  6.600  6.420  6.560  6.310
#I couldn't find a way to bring in the indicator means, 
#so I used the meanstructure = TRUE argument in cfa call (which is not the same thing). 

##################################################################
#Results of the Equal Form Longitudinal Model of Job Satisfaction
##################################################################
#model set up - this is the model presented in Figure 7.2, p.254 
new.model <- 'SATIS1 =~ A1 + B1 + C1 + D1
              SATIS2 =~ A2 + B2 + C2 + D2
              A1 ~~ A2
              B1 ~~ B2
              C1 ~~ C2
              D1 ~~ D2'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=250, meanstructure = TRUE)

#results same as in Table 7.6, page 260 and fit same as in Table 7.7
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

##############################################################################
#Results of the Equal Factor Loadings Longitudinal Model of Job Satisfaction #
##############################################################################
#model set up - this is the model presented in Figure 7.2, p.254 
# loadings constrained to be equal using custom labels
new.model <- 'SATIS1 =~ v1*A1 +  v2*B1 + v3*C1 + v4*D1
              SATIS2 =~ v1*A2 + v2*B2 + v3*C2 + v4*D2
              A1 ~~ A2
              B1 ~~ B2
              C1 ~~ C2
              D1 ~~ D2'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=250,meanstructure = TRUE)

#Fit same as in Table 7.7
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


##############################################################################
#Equal indicator intercepts Longitudinal Model of Job Satisfaction #
##############################################################################
#model set up - this is the model presented in Figure 7.2, p.254 
# loadings constrained to be equal using custom labels
new.model <- 'SATIS1 =~ v1*A1 +  v2*B1 + v3*C1 + v4*D1
              SATIS2 =~ v1*A2 + v2*B2 + v3*C2 + v4*D2
              A1 ~~ A2
              B1 ~~ B2
              C1 ~~ C2
              D1 ~~ D2
              # intercepts constrained to be equal
              # using custom labels
              A1 ~ int2 * 1
              A2 ~ int2 * 1
              B1 ~ int3 * 1
              B2 ~ int3 * 1
              C1 ~ int4 * 1
              C2 ~ int4 * 1
              D1 ~ int5 * 1
              D2 ~ int5 * 1'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=250,meanstructure=TRUE)

#results are not the same as in Table 7.7 - no way to bring in means in lavaan (?)
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


##############################################################################
#Equal indicator error variances Longitudinal Model of Job Satisfaction #
##############################################################################
#model set up - this is the model presented in Figure 7.2, p.254
new.model <- 'SATIS1 =~ v1*A1 +  v2*B1 + v3*C1 + v4*D1
              SATIS2 =~ v1*A2 + v2*B2 + v3*C2 + v4*D2
              A1 ~~ A2
              B1 ~~ B2
              C1 ~~ C2
              D1 ~~ D2
              # intercepts constrained to be equal
              # using custom labels
              A1 ~ int2 * 1
              A2 ~ int2 * 1
              B1 ~ int3 * 1
              B2 ~ int3 * 1
              C1 ~ int4 * 1
              C2 ~ int4 * 1
              D1 ~ int5 * 1
              D2 ~ int5 * 1
              #error variances constrained to be equal
              #using custom labels
              f1*A1 ~~ f1*A1
              f1*A2 ~~ f1*A2
              f2*B1 ~~ f2*B1
              f2*B2 ~~ f2*B2
              f3*C1 ~~ f3*C1
              f3*C2 ~~ f3*C2
              f4*D1 ~~ f4*D1
              f4*D2 ~~ f4*D2'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=250,meanstructure=TRUE)

#results are not the same as in Table 7.7 - no way to bring in means in lavaan (?)
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

#to see the internal representation of the model (this helps understand the model definition)
inspect(fit,what="list")

#############################################
#### Multiple Group Analysis pp. 268-299 ####
#############################################

# read Depression example from the book's companion website
depres <- read.table("http://people.bu.edu/tabrown/Ch7/MDDALL.dat")

# assign variable names
names(depres) <- c("Sex","M1","M2","M3","M4","M5","M6","M7","M8","M9")
head(depres)

#model set up - this is the model presented in Figure 7.2, p.238 
new.model <- 'majdep =~ M1 + M2 + M3 + M4 + M5 + M6 + M7 + M8 + M9
              M1 ~~ M2'

### Separate analysis of males and females to inspect model fit

#model fit for Female
fit <- cfa(new.model, data=depres[Sex==0,-1])
#results are not the same as in Table xxxx
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

#model fit for Male
fit2 <- cfa(new.model, data=depres[Sex==1,-1])
#results are not the same as in Table xxxx
summary(fit2,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

#Tests of Measurement Invariance and Population Heterogeneity

#Measurement Invariance - same model fit as in Table 7.9 
# model 1: configural invariance
fit1 <- cfa(new.model, data=depres, group="Sex")
# model 2: weak invariance
fit2 <- cfa(new.model, data=depres, group="Sex",
            group.equal="loadings")
# model 3: strong invariance
fit3 <- cfa(new.model, data=depres, group="Sex",
            group.equal=c("loadings", "intercepts"))
# model 4: equal loadings + intercepts + residuals
fit4 <- cfa(new.model, data=depres, group="Sex",
            group.equal=c("loadings", "intercepts", "residuals"))
#Population heterogeneity- # same model fit as in Table 7.9 
# model 5: equal loadings + intercepts + residuals + factor variances
fit5 <- cfa(new.model, data=depres, group="Sex",
            group.equal=c("loadings", "intercepts", "residuals","lv.variances"))
# model 6: equal loadings + intercepts + residuals + factor variances + factor means
fit6 <- cfa(new.model, data=depres, group="Sex",
            group.equal=c("loadings", "intercepts", "residuals","lv.variances","means"))

############## Alternatively for Models 1-3 and equal loadings + intercepts + means ##########
############## with chi-square tests of differences ##########
install.packages("semTools")
library(semTools)
measurementInvariance(new.model, data=depres, group="Sex", strict=FALSE)
###########################################

###################################
##### MIMIC Models pp.304-316 #####
###################################

cormat <- '
1.000
0.705   1.000
0.724   0.646   1.000
0.213   0.195   0.190   1.000
0.149   0.142   0.128   0.521   1.000
0.155   0.162   0.135   0.557   0.479   1.000'

Cmat <- getCov(cormat)

neodat.sdev = c(2.26, 2.73, 2.11, 2.32, 2.61, 2.44)
Dmat = diag(neodat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("S1","S2","S3","A1","A2","A3")
colnames(covmat)= rownames(covmat)

# Step 1. Ensure that the two-factor model of Social Phobia
# and Agoraphobia is reasonable and good fitting in the full sample (N = 730).

#model set up - this is the CFA model of the correlation matrix presented in Figure 7.5, p.308 
new.model <- 'socialphobia =~ S1 + S2 + S3   
              agoraphobia =~ A1 + A2 + A3
              socialphobia ~~ agoraphobia'            

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=730)

#results same as p. 308
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


# Step 2. MIMIC Model of Fig. 7.5, p. 308.

#Introduce sex
cormat <- '
1.000
0.705   1.000
0.724   0.646   1.000
0.213   0.195   0.190   1.000
0.149   0.142   0.128   0.521   1.000
0.155   0.162   0.135   0.557   0.479   1.000
-0.019  -0.024  -0.029  -0.110  -0.074  -0.291   1.000'


Cmat <- getCov(cormat)

neodat.sdev = c(2.26, 2.73, 2.11, 2.32, 2.61, 2.44, 0.50)
Dmat = diag(neodat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("S1","S2","S3","A1","A2","A3","sex")
colnames(covmat)= rownames(covmat)

new.model <- 'socialphobia =~ S1 + S2 + S3   
              agoraphobia =~ A1 + A2 + A3
              socialphobia ~~ agoraphobia
              #regressions
              socialphobia ~ sex
              agoraphobia ~ sex'            

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=730)

#results same as p. 315 - TABLE 7.15. Results of MIMIC Model of Social Phobia and Agoraphobia
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)


### Extra: Evaluation of measurement invariance by fixing all direct effects 
### between the covariate and the indicators to zero and then inspecting modification indices.
### pp.314-316

new.model <- 'socialphobia =~ S1 + S2 + S3   
              agoraphobia =~ A1 + A2 + A3
              socialphobia ~~ agoraphobia
              #regressions
              socialphobia ~ sex
              agoraphobia ~ sex
              #fix direct effects between sex and the indicators to zero
              A3 ~ 0*sex
              A2 ~ 0*sex
              A1 ~ 0*sex'            

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=730)

#results same as p. 315 - TABLE 7.15. Results of MIMIC Model of Social Phobia and Agoraphobia
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

Chapter 8 - Higher-Order CFA, Scale Reliability Evaluation, Formative Indicators

Contains examples for a) a Higher Order model CFA b) estimating a scale’s reliability using the semtools package c) a Formative Indicators model

################################################
# Chapter 8: Other Types of CFA Models,        #
#            Higher-Order Factor Analysis      #
#            Scale Reliability and             #
#            Formative Indicators              #
################################################

library(lavaan)

#############################################
# Higher-Order Factor Analysis p.320 - 337  #
#############################################

#Input matrix: the correlation matrix of Figure 8.1
cormat <- '
1.00
0.78  1.00
0.80  0.77  1.00
0.56  0.51  0.48  1.00
0.52  0.51  0.46  0.78  1.00
0.59  0.51  0.51  0.80  0.79  1.00
0.16  0.15  0.17  0.14  0.18  0.16  1.00
0.19  0.13  0.18  0.14  0.16  0.16  0.81  1.00
0.12  0.17  0.17  0.17  0.20  0.16  0.75  0.80  1.00
0.16  0.13  0.17  0.15  0.16  0.18  0.56  0.52  0.50  1.00
0.16  0.14  0.18  0.15  0.16  0.18  0.51  0.58  0.51  0.81  1.00
0.16  0.15  0.14  0.16  0.16  0.14  0.52  0.57  0.52  0.80  0.79  1.00'

Cmat <- getCov(cormat)

solvedat.sdev = c(1.40, 2.10, 1.30, 2.30, 2.40, 1.90, 2.00, 2.90, 2.30, 3.10, 2.20, 1.20)

Dmat = diag(solvedat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("P1","P2","P3","C1","C2","C3","E1","E2","E3","S1","S2","S3")
colnames(covmat)= rownames(covmat)

#model set up - this is the model presented in Figure 8.1, p.322 
new.model <- 'PROBSLV =~ P1 + P2 + P3   
              COGRES =~ C1 + C2 + C3
              EXPREMOT =~ E1 + E2 + E3
              SOCSPT =~ S1 + S2 + S3
              PROBFOC =~ PROBSLV + COGRES
              EMOTFOC =~ EXPREMOT + SOCSPT'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=275)

#results very close to Table 8.3, page 333 
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

###########################################
# Scale Reliability Estimation p. 337-351 #
###########################################

#Input matrix: the correlation matrix of Figure 8.1
cormat <- '
1.000
0.594  1.000
0.607  0.613  1.000
0.736  0.765  0.717  1.000
0.378  0.321  0.360  0.414  1.000
0.314  0.301  0.345  0.363  0.732  1.000
0.310  0.262  0.323  0.337  0.665  0.583  1.000
0.317  0.235  0.276  0.302  0.632  0.557  0.796  1.000'

Cmat <- getCov(cormat)

solvedat.sdev = c(1.150,  1.200,  1.570,  2.820,  1.310,  1.240,  1.330,  1.290)

Dmat = diag(solvedat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("Y1","Y2","Y3","Y4","Y5","Y6","Y7","Y8")
colnames(covmat)= rownames(covmat)

#model set up - this is the model presented in Figure 8.2, p.340 
new.model <- 'INTRUS =~ Y1 + Y2 + Y3 + Y4
              AVOID =~ Y5 + Y6 + Y7 + Y8
              Y7 ~~ Y8'

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=500)

#model fit is close but not identical to p.339 
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

#EXTRA: For scale reliability estimation you can use the reliability function
#of the semTools package
install.packages("semTools")
library(semTools)
reliability(fit) 
#OUTPUT
#       INTRUS     AVOID
#alpha  0.8399990 0.8865353
#omega  0.9193497 0.8232744
#omega2 0.9193497 0.8232744
#omega3 0.9195232 0.8234095
# the omega values equal the point estimates of scale reliability
# shown in p.344, but standard errors or confidence intervals are not calculated.

#For handmade calculations of omega and omega_h you can visit 
#https://bearspace.baylor.edu/Alex_Beaujean/R/Omega.pdf
#for an example

#######################################################
##### Models with Formative Indicators pp.351-362 #####
#######################################################

#Input matrix: the correlation matrix of Figure 8.5
cormat <- '
1.000
0.700  1.000
0.713  0.636  1.000
0.079  0.066  0.076  1.000
0.088  0.058  0.070  0.681  1.000
0.084  0.056  0.074  0.712  0.633  1.000
0.279  0.248  0.240  0.177  0.155  0.170  1.000
0.250  0.214  0.222  0.157  0.143  0.152  0.373  1.000
0.280  0.236  0.251  0.173  0.178  0.171  0.448  0.344  1.000'

Cmat <- getCov(cormat)

solvedat.sdev = c(2.5, 2.1, 3.0, 4.1, 3.9, 4.4, 1.2, 1.0, 1.2)

Dmat = diag(solvedat.sdev)
covmat = Dmat%*%Cmat%*%Dmat

# assign row and column names to the covariance matrix
rownames(covmat) = c("Y1","Y2","Y3","Y4","Y5","Y6","X1","X2","X3")
colnames(covmat)= rownames(covmat)

#Life Stress example
#model set up - this is the model presented in Figure 8.5, p.359 
new.model <- 'SATIS =~ Y1 + Y2 + Y3
OPTIM =~ Y4 + Y5 + Y6
STRESS ~ 1*X1 + X2 +X3  
STRESS =~ NA*SATIS + OPTIM
X1 ~~ X2
X1 ~~ X3
X2 ~~ X3'            

#model fit
fit <- cfa(new.model, sample.cov=covmat, sample.nobs=500)

#results same as in Table 8.9, page 360-361 (only model fit)
summary(fit,standardized=TRUE, fit.measures=TRUE, rsq=TRUE, modindices=TRUE)

#### EXTRA: Visualize this model with semPlot ####
install.packages("semPlot")
library(semPlot)
semPaths(fit, "std", title = FALSE, curvePivot = TRUE)

Chapter 9 - Data Issues in CFA: Missing, Non-normal, and Categorical Data (coming soon)

Chapter 10 - Statistical Power and Sample Size (coming soon)