library(haven)
SCoA_Brazil_NZ_WLSMV <- read_sav("C:/Users/gbro008/Google Drive/2019_08 Moscow summer school/Lectures/resources/SCoA Brazil_NZ WLSMV.sav")
View(SCoA_Brazil_NZ_WLSMV)
#install required libraries
library(lavaan)
library(semPlot)
#to create 8 factor SCoA model, nb no sf4 in Brazil data
SCoA6.modelwls <- '
Bad =~ bd1 + bd2 + bd3 + bd4 + bd5
CE =~ ce1 + ce2 + ce3 + ce4 + ce5 + ce6
IG =~ ig1 + ig2 + ig3
PE =~ pe1 + pe2
SF =~ sf1 + sf2 + sf3
SI =~ si1 + si2 + si3 + si4 + si5
SQ =~ sq1 + sq2
TI =~ ti1 + ti2 + ti3 + ti4 + ti5 + ti6
'

#cfa analysis with WLSmv estimator, ordered command to inform CFA that this is ordered scores
SCoAwls_fit <-cfa (SCoA6.modelwls, data=SCoA_Brazil_NZ_WLSMV, ordered=c("bd1", "bd2", "bd3", "bd4", "bd5", "ce1", "ce2", "ce3", "ce4", "ce5", "ce6", "ig1", "ig2", "ig3", "pe1", "pe2", "sf1", "sf2", "sf3", "si1", "si2", "si3", "si4", "si5", "sq1", "sq2", "ti1", "ti2", "ti3", "ti4", "ti5", "ti6"), estimator="WLS")

#get results of cfa
summary (SCoAwls_fit, standardized=TRUE)
#detailed fit indices
fitmeasures(SCoAwls_fit)
#to get gamma hat fit index
library(semTools)
moreFitIndices(SCoAwls_fit)
#to get x2/df ratio and p-value
fm <- fitMeasures(SCoAwls_fit)
fm[["chisq"]] / fm[["df"]]
fm["pvalue"] 

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

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

#displays MI in decreasing order, and print only those that would substantially improve model fit i.e., changes in x2 square  greater than 20.00  with df = 1)
modind_scoawls <- modindices(SCoAwls_fit, sort. = TRUE)
modind_scoawls[modind_scoawls$mi > 20.00 , ]