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
'

#Model 1 configural invariance by country
SCoAwls_configfit <- 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", "sf4", "si1", "si2", "si3", "si4", "si5", "sq1", "sq2", "ti1", "ti2", "ti3", "ti4", "ti5", "ti6"), meanstructure=TRUE, std.lv=TRUE, group="Country") 
summary(SCoAwls_configfit, fit.measures=TRUE)

#to get gamma hat fit index
library(semTools)
moreFitIndices(SCoAwls_configfit)

#to get x2/df ratio and p-value-thanks to Terry Jorgensen UvA (amsterdam)
fm <- fitMeasures(SCoAwls_configfit)
fm[["chisq"]] / fm[["df"]]
fm["pvalue"] 

#if weak invariance not proven, then look at MI to see if changes could be made to make the model invariant; then run again 
modind_scoawlsconfig <- modindices(SCoAwls_configfit, sort. = TRUE)
modind_scoawlsconfig [modind_scoawlsconfig$mi > 20.00 , ]


###MODEL 2 weak invariance model
SCoAwls_metricfit <- 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", "sf4", "si1", "si2", "si3", "si4", "si5", "sq1", "sq2", "ti1", "ti2", "ti3", "ti4", "ti5", "ti6"), meanstructure=TRUE, std.lv=TRUE, group="Country", group.equal=c("loadings")) 
summary(SCoAwls_metricfit,fit.measures=TRUE)

#checking the significance of chi-square change in R, model 2 v.s model 1:
anova(SCoAwls_metricfit, SCoAwls_configfit)

###MODEL 3 strong invariance model
SCoAwls_scalarfit <- 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", "sf4", "si1", "si2", "si3", "si4", "si5", "sq1", "sq2", "ti1", "ti2", "ti3", "ti4", "ti5", "ti6"), meanstructure=TRUE, std.lv=TRUE, group="Country", group.equal=c("loadings","intercepts")) 
summary(SCoAwls_scalarfit,fit.measures=TRUE)

#checking the significance of chi-square change in R, model 2 v.s model 1:
anova(SCoAwls_scalarfit, SCoAwls_metricfit)

###MODEL 4 strict invariance model
SCoAwls_strictfit <- 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", "sf4", "si1", "si2", "si3", "si4", "si5", "sq1", "sq2", "ti1", "ti2", "ti3", "ti4", "ti5", "ti6"), meanstructure=TRUE, std.lv=TRUE, group="Country", group.equal=c("loadings","intercepts", "residuals")) 
summary(SCoAwls_strictfit,fit.measures=TRUE)
#checking the significance of chi-square change in R, model 2 v.s model 1:
anova(SCoAwls_strictfit, SCoAwls_scalarfit)

#detailed fit indices
fitmeasures(SCoAwls_configfit)

#Model last full invariance
SCoAwls_fullfit <- 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", "sf4", "si1", "si2", "si3", "si4", "si5", "sq1", "sq2", "ti1", "ti2", "ti3", "ti4", "ti5", "ti6"), meanstructure=TRUE, std.lv=TRUE, group="Country", group.equal=c(group.equal=c("loadings","intercepts","residuals", "lv.variances","lv.covariances","means")
summary(SCoAwls_fullfit,fit.measures=TRUE)