---
title: "R Codes"
author: "Bianca G Martins"
date: "22 de abril de 2020"
output: html_document
---

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

## read database
bd <- read_excel("enter the address for the file here", sheet = "database")


## Define items analysis

it = c("PEVI1","PEVI2","PEVI3","PEVI4","PEVI5","PEVI6","PEVI7","PEVI8","PEVI9","PEVI10","PEVI11","PEVI12","PEVI13","PEVI14","PEVI15")


## Specify Factorial Model to CFA
#NUT = Nutrition; PAC = Physical Activity; PRB = Preventive Behaviors; SOR = Social Relationships; STM = Stress Management; PHY = Physical and PSY = Psychological (second-order factors) 
#Modelo1 = First Order; #Modelo2 = Second Order
modelo1 = '
NUT =~ PEVI1 + PEVI2 + PEVI3
PAC =~ PEVI4 + PEVI5 + PEVI6
PRB =~ PEVI7 + PEVI8 + PEVI9
SOR =~ PEVI10+ PEVI11 + PEVI12
STM =~ PEVI13 + PEVI14 + PEVI15
'

modelo2 = '
NUT =~ PEVI1 + PEVI2 + PEVI3
PAC =~ PEVI4 + PEVI5 + PEVI6
PRB =~ PEVI7 + PEVI8 + PEVI9
SOR =~ PEVI10+ PEVI11 + PEVI12
STM =~ PEVI13 + PEVI14 + PEVI15
PHY =~ NUT + PAC + PRB
PSY =~ SOR + STM
'

##Specify Structural Model
#Modelo 3 = Complete Structural Model
modelo3 = '
NUT =~ PEVI1 + PEVI2 + PEVI3
PAC =~ PEVI4 + PEVI5 + PEVI6
PRB =~ PEVI7 + PEVI8 + PEVI9
SOR =~ PEVI10+ PEVI11 + PEVI12
STM =~ PEVI13 + PEVI14 + PEVI15
PHY =~ NUT + PAC + PRB
PSY =~ SOR + STM
PHY ~ Sex + WSdic + CAreadic + Age
PSY ~ Sex + WSdic + CAreadic + Age
'
```

##Confirmatory Factor Analysis
```{r}
#fit the model CFA
fit1 = cfa(modelo1, data=bd,std.lv=T, ordered = it, estimator="WLSMV", group = "Sample")
fit2 = cfa(modelo2, data=bd,std.lv=T, ordered = it, estimator="WLSMV", group = "Sample")

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

summary(fit2, fit.measures=TRUE, standardized=TRUE)
fitmeasures(fit2, fit.measures = "all", baseline.model = NULL)

```

##Modification indexes (MI)
```{r}
#MI
im1 = modindices(fit1, sort.=TRUE, minimum.value=4)
im2 = modindices(fit2, sort.=TRUE, minimum.value=4)

#MI for fator loadings
selIM=im1[im1$op == "~~",]
knitr::kable(selIM)

selIM2=im2[im2$op == "~~",]
knitr::kable(selIM)
```

##Reliability
```{r}
#Reliability
fiab <- round(reliability(fit1),3)#First order only
knitr::kable(fiab)

```

##Invariance analysis
###Group = Sample (0 - Test, 1 - Validation)
```{r}
grupo ="Sample"
test.seq <- c("thresholds","loadings","regression", "means","residuals")
meq.list <- list()
for (i in 0L:length(test.seq)) {
  if (i == 0L) {
    meq.label <- "configural"
    group.equal <- ""
  } else {
    meq.label <- test.seq[i]
    group.equal <- test.seq[1:i]
  }
  meq.list[[meq.label]] <- measEq.syntax(configural.model = modelo1,
                                         ridge=1e-05,
                                         data = bd,
                                         ordered = it,
                                         parameterization = "theta",
                                         #ID.fac = "ul",
                                         ID.cat = "Wu.Estabrook.2016",
                                         group = grupo,
                                         group.equal = group.equal,
                                         return.fit = TRUE)
}

compareFit(meq.list)
```


##Structural Model Analysis
```{r}
# Fit the Structural Model
fit3 = sem(modelo3, data=bd,std.lv=T, ordered = it, estimator="WLSMV")

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