---
title: "OK"
output:
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
##Libraries
```{r}
# Libraries/Packages
library(foreign)
library(lavaan)
library(semTools)
library(haven)
library(readxl)

```

##Database and analysis
```{r}
#BD
bd0 <- read_excel("enter the address for the file here", sheet = "General Database")


## Define items analysis
#it1 = c("EMAQ1","EMAQ2","EMAQ3","EMAQ4","EMAQ5","EMAQ6","EMAQ7","EMAQ8","EMAQ9","EMAQ10","EMAQ11","EMAQ12","EMAQ13","EMAQ14","EMAQ15","EMAQ16","EMAQ17","EMAQ18","EMAQ19","EMAQ20","EMAQ21","EMAQ22")

## Specify Factorial Model to CFA
#Model 0 = Bifactorial Model (without items 10 and 22)
model0 = '
NEG =~ EMAQ1 + EMAQ2 + EMAQ4 + EMAQ5 + EMAQ7 + EMAQ8 + EMAQ9 + EMAQ13 + EMAQ15 + EMAQ16 + EMAQ17 + EMAQ19 + EMAQ21
POS =~ EMAQ3 + EMAQ6 + EMAQ11 + EMAQ12 + EMAQ14 + EMAQ18 + EMAQ20
'
#Model 1 = Refined Bifactorial Model (without items 8, 10, 13, 18 and 22)
model1 = '
NEG =~ EMAQ1 + EMAQ2 + EMAQ5 + EMAQ7 + EMAQ9 + EMAQ15 + EMAQ16 + EMAQ17 + EMAQ19 + EMAQ21
POS =~ EMAQ3 + EMAQ6 + EMAQ11 + EMAQ12 + EMAQ14 + EMAQ20
'

#fit the model CFA
fit1 = cfa(model1, data=bd0, std.lv=T, estimator="DWLS", se="bootstrap", group="Sample")

#factor loadings
summary(fit1, fit.measures=TRUE, standardized=TRUE)
fitmeasures(fit1, fit.measures = "all", baseline.model = NULL)

#Indices de Modificação
im1 = modindices(fit1, sort.=TRUE, minimum.value=4)
selIM1=im1[im1$op == "~~",]
knitr::kable(selIM1)

## Reliability
reliability(fit1)

##Factorial Invariance Analysis
it = c("EMAQ1","EMAQ2","EMAQ3","EMAQ4","EMAQ5","EMAQ6","EMAQ7","EMAQ8","EMAQ9","EMAQ10","EMAQ11","EMAQ12","EMAQ13","EMAQ14","EMAQ15","EMAQ16","EMAQ17","EMAQ18","EMAQ19","EMAQ20","EMAQ21","EMAQ22")

fit1.config = cfa(model1, data=bd0,std.lv=F, auto.cov.lv.x=T, ordered=it, estimator="DWLS", group = "Sample")
fit1.metric = cfa(model1, data=bd0,std.lv=F, auto.cov.lv.x=T, ordered=it, estimator="DWLS", group = "Sample", group.equal = "loadings")
fit1.scalar = cfa(model1, data=bd0,std.lv=F, auto.cov.lv.x=T, ordered=it, estimator="DWLS", group = "Sample", group.equal = c("loadings","thresholds"))
fit1.strict = cfa(model1, data=bd0,std.lv=F, auto.cov.lv.x=T, ordered=it, estimator="DWLS", group = "Sample", group.equal =c("loadings","thresholds", "residuals"))

measEqOut.1 <- compareFit(fit1.config, fit1.metric, fit1.scalar, fit1.strict)
summary(measEqOut.1)
summary(measEqOut.1, fit.measures = "all")
summary(measEqOut.1, fit.measures = c("cfi", "cfi.scaled"))

##Polychoric correlation between EMAQ and TFEQ factors

bd1 <- read_excel("enter the address for the file here", sheet = "EMAQ TFEQ Sample 1")
bd2 <- read_excel("enter the address for the file here", sheet = "EMAQ TFEQ Sample 2")
bd3 <- read_excel("enter the address for the file here", sheet = "EMAQ TFEQ Sample 3")

it2 = c("EMAQ1","EMAQ2","EMAQ3","EMAQ4","EMAQ5","EMAQ6","EMAQ7","EMAQ8","EMAQ9","EMAQ10","EMAQ11","EMAQ12","EMAQ13","EMAQ14","EMAQ15","EMAQ16","EMAQ17","EMAQ18","EMAQ19","EMAQ20","EMAQ21","EMAQ22","TFEQ1","TFEQ2","TFEQ3","TFEQ4","TFEQ5","TFEQ6","TFEQ7","TFEQ8","TFEQ9","TFEQ10","TFEQ11","TFEQ12","TFEQ13","TFEQ14","TFEQ15","TFEQ16","TFEQ17","TFEQ18")

model2 = '
RC =~ TFEQ2 + TFEQ11 + TFEQ12 + TFEQ15 + TFEQ16 + TFEQ18
DA =~ TFEQ1 + TFEQ4 + TFEQ5 + TFEQ7 + TFEQ8 + TFEQ9 + TFEQ13 + TFEQ14 + TFEQ17
AE =~ TFEQ3 + TFEQ6 + TFEQ10
NEG =~ EMAQ1 + EMAQ2 + EMAQ4 + EMAQ5 + EMAQ7 + EMAQ8 + EMAQ9 + EMAQ10 + EMAQ13 + EMAQ15 + EMAQ16 + EMAQ17 + EMAQ19 + EMAQ21
POS =~ EMAQ3 + EMAQ6 + EMAQ11 + EMAQ12 + EMAQ14 + EMAQ18 + EMAQ20 + EMAQ22
NEG ~~ RC
POS ~~ RC
'
model3 = '
RC =~ TFEQ2 + TFEQ11 + TFEQ12 + TFEQ15 + TFEQ16 + TFEQ18
DA =~ TFEQ1 + TFEQ4 + TFEQ5 + TFEQ7 + TFEQ8 + TFEQ9 + TFEQ13 + TFEQ14 + TFEQ17
AE =~ TFEQ3 + TFEQ6 + TFEQ10
NEG =~ EMAQ1 + EMAQ2 + EMAQ4 + EMAQ5 + EMAQ7 + EMAQ8 + EMAQ9 + EMAQ10 + EMAQ13 + EMAQ15 + EMAQ16 + EMAQ17 + EMAQ19 + EMAQ21
POS =~ EMAQ3 + EMAQ6 + EMAQ11 + EMAQ12 + EMAQ14 + EMAQ18 + EMAQ20 + EMAQ22
NEG ~~ DA
POS ~~ DA
'
model4 = '
RC =~ TFEQ2 + TFEQ11 + TFEQ12 + TFEQ15 + TFEQ16 + TFEQ18
DA =~ TFEQ1 + TFEQ4 + TFEQ5 + TFEQ7 + TFEQ8 + TFEQ9 + TFEQ13 + TFEQ14 + TFEQ17
AE =~ TFEQ3 + TFEQ6 + TFEQ10
NEG =~ EMAQ1 + EMAQ2 + EMAQ4 + EMAQ5 + EMAQ7 + EMAQ8 + EMAQ9 + EMAQ10 + EMAQ13 + EMAQ15 + EMAQ16 + EMAQ17 + EMAQ19 + EMAQ21
POS =~ EMAQ3 + EMAQ6 + EMAQ11 + EMAQ12 + EMAQ14 + EMAQ18 + EMAQ20 + EMAQ22
NEG ~~ AE
POS ~~ AE
'

# Fit the Structural Model
fit2 = sem(model2, data=bd1,std.lv=T, ordered = it2, estimator="WLSMV")

#Factor loadings
summary(fit2, fit.measures=TRUE, standardized=TRUE)
inspect(fit2, 'r2')

```