#SEM model mimic the path model

#to create bifactor 4 main factors SCoA model with Assessment checklist
#create latent factors
SCoA6Checkfull.model <- '
BD =~ bd1 + bd2 + bd3 + bd4 + bd5 + ig1 + ig2 + ig3
IMP =~ si1 + si2 + si3 + si4 + si5 + ti1 + ti2 + ti3 + ti4 + ti5 + ti6 
AFFECT =~ ce1 + ce2 + ce3 + ce4 + ce5 + ce6 + pe1 + pe2
EXTRNL =~ sq1 + sq2 + sf1 + sf2 + sf3 + sf4
general =~ bd1 + bd2 + bd3 + bd4 + bd5 + ce1 + ce2 + ce3 + ce4 + ce5 + ce6 + ig1 + ig2 + ig3 + pe1 + pe2 + sf1 + sf2 + sf3 + si1 + si2 + si3 + sf4 + si4 + si5 + sq1 + sq2 + ti1 + ti2 + ti3 + ti4 + ti5 + ti6
Check1 =~ A_Check1 + A_Check5 + A_Check6 + A_Check7 + A_Check10
Check2 =~ A_Check2 + A_Check3 + A_Check4 + A_Check8 + A_Check9 + A_Check11 + A_Check12 

# specify correlations among factors
     general ~~ 0*BD + 0*IMP + 0*AFFECT + 0*EXTRNL
	 Check1 ~~ Check2
	 BD ~~ IMP + AFFECT + EXTRNL
	 IMP ~~ AFFECT + EXTRNL
	 AFFECT ~~ EXTRNL
#specify regressions from path model
	Check1 + Check2 ~ general
	Check1 + Check2 ~ BD
	Check1 + Check2 ~ AFFECT
	Check1 + Check2 ~ IMP
	Check1 + Check2 ~ EXTRNL
'

# Fit the model to the data using sem instead of cfa. The difference  lies in the default settings.
fit_sem <- sem(model = SCoA6Checkfull.model, data = NZSCoAVI, ordered=c("A_Check1","A_Check5", "A_Check6","A_Check7", "A_Check10", "A_Check2", "A_Check3", "A_Check4", "A_Check8", "A_Check9", "A_Check11", "A_Check12", bd1", "bd2", "bd3", "bd4", "bd5", "ce1", "ce2", "ce3", "ce4", "ce5", "ce6", "ig1", "ig2", "ig3", "pe1", "pe2", "sf1", "sf2", "sf3", "sf4", "si1", "si2", "si3", "si4", "si5", "sq1", "sq2", "ti1", "ti2", "ti3", "ti4", "ti5", "ti6"), estimator="WLS")

# Show the (minimal) output
summary(fit_sem, standardized = TRUE)
#full fit info
fitmeasures(fit_sem)
#to get gamma hat fit index
library(semTools)
moreFitIndices(fit_sem)

#to get x2/df ratio and p-value-thanks to Terry Jorgensen UvA (amsterdam)
fm <- fitMeasures(AssessCheck_fit)
fm[["chisq"]] / fm[["df"]]
fm["pvalue"] 
#to get conventional correlated factors diagram; preferred method
library(semPlot)
semPaths(fit_sem, intercept = FALSE, whatLabel = "std",  residuals = FALSE, exoCov = TRUE) 

#get a picture,standardised scores shown; this not preferred for this model
library(lavaanPlot)
lavaanPlot(model = fit_sem, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"), coefs = TRUE, stand=TRUE)