library(lavaan)
library(survey)
library(lavaan.survey)
library(semPlot)
library(semTools)
library(psych)
library(GPArotation)
library(haven)

library(haven)
asertiv <- read_sav("D:/Disco Local D/2019_publicaciones/articulos para publicar/asertividad y estilos de crianza/asertiv.sav")
View(asertiv)
names(asertiv)

####MODELO BASE DE DOS FACTORES DE ASERTIVIDAD #######################################
model.1 <- "
auto =~ As_1 + As_2 + As_3 + As_4 + As_5 + As_6 + As_7 + As_8 + As_9 + As_10 + 
         As_11 + As_12 + As_13 + As_14 + As_15 + As_16 + As_17 + As_18 + As_19 + As_20  
hetero =~ As_22 + As_23 + As_24 + As_25 + As_26 + As_27 + As_28 + As_29 + As_30 + As_31 + 
         As_32 + As_33 + As_34 + As_35
"

###Robustos de minimos cuadrados no ponderados polychoric with data not specific sin errores correlacionados
fit.model.1 <- cfa(model.1, data=asertiv, estimator="WLSMV", missing = "listwise", 
                   ordered = c("As_1", "As_2", "As_3", "As_4", "As_5", "As_6", "As_7", "As_8",
                               "As_9", "As_10", "As_11", "As_12", "As_13", "As_14", "As_15", "As_16",
                               "As_17", "As_18", "As_19", "As_20", "As_21", "As_22", "As_23", "As_24", "As_25",
                               "As_26", "As_27", "As_28", "As_29", "As_30", "As_31",
                               "As_32", "As_33", "As_34", "As_35"))

summary(fit.model.1, fit.measures=TRUE, standardized=TRUE)

###Indices de modificaci?n Robustos de minimos cuadrados no ponderados polychoric sin erores correlacionados
mi.fit.model.1 <- modindices(fit.model.1)
print(mi.fit.model.1)



semPaths(fit.model.1, whatLabels="std", curvePivot = TRUE, layout="tree", edge.label.cex=1.0, 
         sizeInt = 12, rotation=2, sizeLat=8, sizeMan=10, sizeMan2 = 2, style="lisrel",
         thresholds = F, intercepts = F, residuals = T, groups = "latents",  pastel = TRUE)


############### MODELO DE DOS FACTORES DE ASERTIVIDAD, MODELO 3, MODELO FINAL
model.2 <- "
auto =~ As_4 + As_6 + As_7 + As_8 + As_9 +  
         As_11 + As_12 + As_16 + As_20  
hetero =~ As_22 + As_23 + As_24 + As_27 + As_28 + As_29 + As_30 + As_31 + 
         As_32 + As_34 + As_35
"

###Robustos de minimos cuadrados no ponderados polychoric with data not specific sin errores correlacionados
fit.model.2 <- cfa(model.2, data=asertiv, estimator="WLSMV", missing = "listwise", 
                   ordered = c("As_4", "As_6", "As_7", "As_8",
                               "As_9", "As_11", "As_12", "As_16",
                               "As_20", "As_21", "As_22", "As_23", "As_24", 
                               "As_27", "As_28", "As_29", "As_30", "As_31",
                               "As_32", "As_34", "As_35"))

summary(fit.model.2, fit.measures=TRUE, standardized=TRUE)

###Indices de modificaci?n Robustos de minimos cuadrados no ponderados polychoric sin erores correlacionados
mi.fit.model.2 <- modindices(fit.model.2)
print(mi.fit.model.2)



semPaths(fit.model.2, whatLabels="std", curvePivot = TRUE, layout="tree", edge.label.cex=1.0, 
         sizeInt = 12, rotation=2, sizeLat=8, sizeMan=10, sizeMan2 = 2, style="lisrel",
         thresholds = F, intercepts = F, residuals = T, groups = "latents",  pastel = TRUE)

###################### invarianza factorial

#####invarinza factorial
library(lavaan)
library(semPlot)

#crea el modelo DE INVARIANZA POR SEXO
model.3 <- "
auto =~ As_4 + As_6 + As_7 + As_8 + As_9 +  
         As_11 + As_12 + As_16 + As_20  
hetero =~ As_22 + As_23 + As_24 + As_27 + As_28 + As_29 + As_30 + As_31 + 
         As_32 + As_34 + As_35
"
#Invariance

inv.sex.conf <- measEq.syntax(configural.model = model.3,estimator="WLSMV", ID.fac = "std.lv", parameterization = "theta",  group = "sex", orthogonal=FALSE, data=asertiv,
                              ID.cat = "Wu.Estabrook.2016",return.fit=TRUE,group.equal = c("thresholds"))

summary(inv.sex.conf, fit.measures=TRUE)


inv.sex.metric<- measEq.syntax(configural.model = model.3,estimator="WLSMV", ID.fac = "std.lv", parameterization = "theta",  group = "sex", orthogonal=FALSE, data=asertiv, parameterization = "theta",
                               ID.cat = "Wu.Estabrook.2016",return.fit=TRUE, group.equal = c("thresholds","loadings"),long.equal  = c("thresholds","loadings"))

summary(inv.sex.metric, fit.measures=TRUE)

inv.sex.scalar<- measEq.syntax(configural.model = model.3,estimator="WLSMV", ID.fac = "std.lv", parameterization = "theta",  group = "sex", orthogonal=FALSE, data=asertiv,
                               ID.cat = "Wu.Estabrook.2016",return.fit=TRUE, group.equal = c("thresholds","loadings","intercepts"),long.equal  = c("thresholds","loadings","intercepts"))

summary(inv.sex.scalar, fit.measures=TRUE)

inv.sex.stric<- measEq.syntax(configural.model = model.3,estimator="WLSMV", ID.fac = "std.lv", parameterization = "theta",  group = "sex", orthogonal=FALSE, data=asertiv,
                              ID.cat = "Wu.Estabrook.2016",return.fit=TRUE, group.equal = c("thresholds","loadings","intercepts","residuals"),long.equal  = c("thresholds","loadings","intercepts","residuals"))


summary(inv.sex.stric, fit.measures=TRUE)


#datos mejor organizados
lavaan::anova(inv.sex.stric,inv.sex.scalar,inv.sex.metric, inv.sex.conf)

fit.stats <- rbind(fitmeasures(inv.sex.conf, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")),
                   fitmeasures(inv.sex.metric, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")),
                   fitmeasures(inv.sex.scalar, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")),
                   fitmeasures(inv.sex.stric, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")))
rownames(fit.stats) <- c("configural", "metric","scalar", "strict")
fit.stats

##### invarianza por edad

model.4 <- "
auto =~ As_4 + As_6 + As_7 + As_8 + As_9 +  
         As_11 + As_12 + As_16 + As_20  
hetero =~ As_22 + As_23 + As_24 + As_27 + As_28 + As_29 + As_30 + As_31 + 
         As_32 + As_34 + As_35
"
#Invariance

inv.sex.conf <- measEq.syntax(configural.model = model.4,estimator="WLSMV", ID.fac = "std.lv", parameterization = "theta",  group = "edd", orthogonal=FALSE, data=asertiv,
                              ID.cat = "Wu.Estabrook.2016",return.fit=TRUE,group.equal = c("thresholds"))

summary(inv.sex.conf, fit.measures=TRUE)


inv.sex.metric<- measEq.syntax(configural.model = model.4,estimator="WLSMV", ID.fac = "std.lv", parameterization = "theta",  group = "edd", orthogonal=FALSE, data=asertiv, parameterization = "theta",
                               ID.cat = "Wu.Estabrook.2016",return.fit=TRUE, group.equal = c("thresholds","loadings"),long.equal  = c("thresholds","loadings"))

summary(inv.sex.metric, fit.measures=TRUE)

inv.sex.scalar<- measEq.syntax(configural.model = model.4,estimator="WLSMV", ID.fac = "std.lv", parameterization = "theta",  group = "edd", orthogonal=FALSE, data=asertiv,
                               ID.cat = "Wu.Estabrook.2016",return.fit=TRUE, group.equal = c("thresholds","loadings","intercepts"),long.equal  = c("thresholds","loadings","intercepts"))

summary(inv.sex.scalar, fit.measures=TRUE)

inv.sex.stric<- measEq.syntax(configural.model = model.4,estimator="WLSMV", ID.fac = "std.lv", parameterization = "theta",  group = "edd", orthogonal=FALSE, data=asertiv,
                              ID.cat = "Wu.Estabrook.2016",return.fit=TRUE, group.equal = c("thresholds","loadings","intercepts","residuals"),long.equal  = c("thresholds","loadings","intercepts","residuals"))

summary(inv.sex.stric, fit.measures=TRUE)

#datos mejor organizados
lavaan::anova(inv.sex.stric,inv.sex.scalar,inv.sex.metric, inv.sex.conf)

fit.stats <- rbind(fitmeasures(inv.sex.conf, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")),
                   fitmeasures(inv.sex.metric, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")),
                   fitmeasures(inv.sex.scalar, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")),
                   fitmeasures(inv.sex.stric, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")))
rownames(fit.stats) <- c("configural", "metric","scalar", "strict")
fit.stats

#Partial Invariance

model6<-cfa(model.4, orthogonal=FALSE, data=satisfaccion_afc, estimator="WLSMV",group="grado", group.equal=c("loadings","intercepts"))
summary(model6,fit.measures=TRUE)
lavTestScore(model6)
d<-parTable(model6)
modindices(model6)
model6partial<-cfa(model.4, orthogonal=FALSE, data=satisfaccion_afc, estimator="WLSMV",group="grado", group.equal=c("loadings","intercepts"), group.partial=c("sv6~1"))
summary(model6partial,fit.measures=TRUE)