#install required libraries
library(lavaan)
library(semPlot)

# hierachical CFA
SCoA.hier <-  '
BD =~ bd1 + bd2 + bd3 + bd4 + bd5
IG =~ ig1 + ig2 + ig3
SI =~ si1 + si2 + si3 + si4 + si5
TI =~ ti1 + ti2 + ti3 + ti4 + ti5 + ti6 
CE =~ ce1 + ce2 + ce3 + ce4 + ce5 + ce6
PE =~ pe1 + pe2
SQ =~ sq1 + sq2
SF =~ sf1 + sf2 + sf3
# higher order factors
Irrelevance =~ BD + IG
Improve  =~ SI + TI
Affective =~ CE + PE
External =~ SF + SQ
'
#do CFA hierarchical analysis of model
SCoAHR_fit <- cfa (SCoA.hier, data=NZSCoAVI)
#get results of cfa
summary (SCoAHR_fit, standardized=TRUE)
#detailed fit indices
fitmeasures(SCoAHR_fit)

#to get gamma hat fit index
library(semTools)
moreFitIndices(SCoAHR_fit)
#to get x2/df ratio and p-value
fm <- fitMeasures(SCoAHR_fit)
fm[["chisq"]] / fm[["df"]]
fm["pvalue"] 
#to get conventional correlated factors diagram
semPaths(SCoAHR_fit, intercept = FALSE, whatLabel = "std",  residuals = FALSE, exoCov = TRUE) 

#to get lavaanplot style chart
library(lavaanPlot)
lavaanPlot(model = SCoAHR_fit, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"), coefs = TRUE, covs = TRUE, stars = TRUE)

# hierachical CFA adjusted to cope with negative error variance in IG
SCoA.hier2 <-  '
BD =~ bd1 + bd2 + bd3 + bd4 + bd5
SI =~ si1 + si2 + si3 + si4 + si5
TI =~ ti1 + ti2 + ti3 + ti4 + ti5 + ti6 
CE =~ ce1 + ce2 + ce3 + ce4 + ce5 + ce6
PE =~ pe1 + pe2
SQ =~ sq1 + sq2
SF =~ sf1 + sf2 + sf3
# higher order factors
Irrelevance =~ BD + ig1 + ig2 + ig3
Improve  =~ SI + TI
Affective =~ CE + PE
External =~ SF + SQ
'  
#do CFA hierarchical analysis of model
SCoAHR2_fit <- cfa (SCoA.hier2, data=NZSCoAVI)
#get results of cfa
summary (SCoAHR2_fit, standardized=TRUE)
#detailed fit indices
fitmeasures(SCoAHR2_fit)

#to get conventional correlated factors diagram
semPaths(SCoAHR2_fit, intercept = FALSE, whatLabel = "std",  residuals = FALSE, exoCov = TRUE) 

#to get lavaanplot style chart
library(lavaanPlot)
lavaanPlot(model = SCoAHR2_fit, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"), coefs = TRUE, covs = TRUE, stars = TRUE)