---
title: "R Notebook"
output:
  html_notebook: default
  html_document:
    df_print: paged
editor_options:
  chunk_output_type: console
---

## Código R dos Exemplos

### Capitulo 8 - AMOS

```{r Cap 8, echo=F}
library(lavaan)

MBI <-  read.table('D:/OneDrive - NOVASBE/AMOS/exercicios/MBI.dat', header = TRUE)
model <-
  '
	  Exaustao =~ It1 + It2 + It3 + It4 + It5
    Descrenca =~ It6 + It7 + It8 + It9
    Realizacao =~ It10 + It11 + It12 + It13 + It14 + It15
	'
fit <- sem(model, data = MBI, ordered = TRUE)
summary(fit, fit.measures = TRUE)
standardizedSolution(fit)

```

### Capitulo 9- Regressão Linear e Path Analysis

#### Exemplo 9.2. Evolução do PIB nos USA

```{r ex 9.2, eval=FALSE, include=F}
# biblioteca para análise de equações estruturais
library(lavaan)
#Define a matriz de covariâncias
cov <-  '
  4256.530
  449.016   52.984
  1535.097  139.449   1114.447
  537.482   53.291    170.024   73.747'
#Converte matriz para formato de leitura do Lavaan
bdCov <- getCov(cov, names = c("gnp", "labor", "capital", "time"))
#Define o modelo. A regressão é indicada por ‘~’
modelo <- ' gnp ~ labor + capital + time '
# Ajusta o modelo usando a função genérica ‘sem’
fitModelo <- sem(modelo, sample.cov = bdCov, sample.nobs = 23)
# Produz as estimativas estandardizadas com ‘std=True’
summary(fitModelo, std = TRUE, rsquare = TRUE)
# Representação gráfica do modelo
library(lavaanPlot)
lavaanPlot(model = fitModelo,
           coefs = TRUE,
           stand = TRUE)
```

#### Exemplo 9.4 Desempenho escolar vs investimento

```{r ex 9.4, eval=FALSE, include=FALSE}
# Biblioteca para ler ficheiros do SPSS
library(haven)
bd <-
  read_sav("D:/OneDrive - NOVASBE/ReportNumber/Site/aee/bds/SATvsCOST.sav")
# Bibiloteca para An. Eq. Estruturais
library(lavaan)
# Modelo
modelo <- '  VERBAL + MATH ~ COST + RATIO + SALARY'
# Ajusta o modelo, com cálculo da ordenada na origem (b0) usando meanstructure=TRUE
fitModelo <- sem(model = modelo,
                 data = bd,
                 meanstructure = TRUE)
# Produz output com estimativas estandardizadas e R2
summary(object = fitModelo, std = T, rsq = T)
# Representação gráfica do modelo
library(lavaanPlot)
lavaanPlot(
  model = fitModelo,
  coefs = TRUE,
  stand = TRUE,
  graph_options = list(rankdir = "LR"),
  stars = "regress"
)

modelo <- lm(MATH ~ COST + RATIO + SALARY, data = bd)
car::vif(modelo)

```

#### Exemplo 9.6 Decisão para realizar estudos

```{r ex 9.6, echo=F}
# Lê base de dados do AMOS
library(haven)
bd <- read_sav("D:/OneDrive - NOVASBE/ReportNumber/Site/aee/bds/Path-Ingram.sav")
# Bibiloteca para An. Eq. Estruturais
library(lavaan)
# Define o modelo saturado
modelo <- '
Intent ~ Attitude + SubNorm + PBC
Behavior ~ Attitude + SubNorm + PBC + Intent
'
modelo <- '
Intent ~ Attitude + SubNorm + PBC
Behavior ~ PBC + Intent
'

# Ajusta o modelo
fitModelo <- sem(modelo, data = bd)
# Produz output com estimativas estandardizadas e R^2
summary(fitModelo, std = T, rsq = T)
# Representação gráfica do modelo
library(lavaanPlot)
lavaanPlot(
  model = fitModelo,
  coefs = TRUE,
  stand = TRUE,
  graph_options = list(rankdir = "LR"),
  stars = "regress",
  digits = 3
)
#Seleção de indíces de qualidade de ajustamento
fitMeasures(fitModelo,
            fit.measures = c("cfi", "tli", "srmr", "rmsea", "rmsea.pvalue"))

modelo <- '
Intent ~ a*Attitude + SubNorm + PBC
Behavior ~ PBC + b*Intent
med := a*b
'
fitModelo <- sem(modelo, data = bd)
# Produz output com estimativas estandardizadas e R^2
summary(fitModelo, std = T, rsq = T)
```

### Exemplo 9.8 Preço das casas

```{r ex 9.8, echo=F}
# Importar bd do SPSS
library(haven)
bd <-
  read_sav("D:/OneDrive - NOVASBE/ReportNumber/Site/aee/bds/PrecoCasas.sav")
# Centrar as variaveis com a função scale()
bd[, 1:4] <- scale(bd[, 1:4], scale = F)
#Adiciona o termo de interação
bd$LocalXArea <- bd$Local * bd$Area
# Bibiloteca para An. Eq. Estruturais
library(lavaan)
# Define o modelo saturado
modelo <- 'Preço ~ Idade + Area + IMI + Local + LocalXArea'
# Ajusta o modelo
fitModelo <- sem(modelo, data = bd)
# Produz output com estimativas estandardizadas e R^2
summary(fitModelo, std = T, rsquare = T)
# Representação gráfica do modelo
library(lavaanPlot)
lavaanPlot(
  model = fitModelo,
  coefs = TRUE,
  stand = TRUE,
  graph_options = list(rankdir = "LR"),
  stars = "regress",
  digits = 3
)

```

## Capitulo 10

### Exemplo 10.2

```{r ex 10.2, echo=F}

#Ler a base de dados em SPSS
library(haven)
bd <-
  read_sav("D:/OneDrive - NOVASBE/ReportNumber/Site/aee/bds/PTSD.sav")
# Bibiloteca para An. Eq. Estruturais
library(lavaan)
# Define o modelo saturado. Note o operado '=~' para definir os fatores
modelo <- '
    Rexperiência =~ ptsd1 + ptsd2 + ptsd3 + ptsd4 + ptsd5
    Evitamento =~ ptsd6 + ptsd7 + ptsd8 + ptsd9 + ptsd10 +ptsd16 + ptsd17
    Agitação =~ ptsd11 +ptsd12 + ptsd13 + ptsd14 + ptsd15
    '
# Ajusta o modelo. Por defeito o estimador é ML
fitModelo <- sem(modelo, data = bd)
# Produz output com estimativas estandardizadas e os indices de ajustamento usuais
summary(fitModelo, std = T, fit.measures = T)
# Representação gráfica do modelo
library(lavaanPlot)
lavaanPlot(
  model = fitModelo,
  coefs = TRUE,
  stand = TRUE,
  covs = TRUE,
  graph_options = list(rankdir = "LR")
)

# Estudo da distribuição dos itens
library(MVN)
result <-
  mvn(data = bd [, 5:21], desc = T) # usando as colunas 5 a 21 da bd[,5:21]
# Estatisticas descritivas dos itens
result$Descriptives[, c(1, 4:10)]
# Testes à normalidade univariada com o teste de Shapiro-Wilk
result$univariateNormality
# Teste à normalidade multivariada com o teste de Mardia
result$multivariateNormality
# Para obter histogramas dos primeiros 6 itens
mvn(bd[, 5:10], univariatePlot = "histogram")


#Estimação do Modelo com Estimador Robusto
fitModeloR <- sem(modelo, data = bd, estimator = "MLR")
# Produz output com estimativas estandardizadas
summary(fitModeloR, std = T)
# Seleção de indíces de qualidade de ajustamento
fitMeasures(fitModeloR,
            c("cfi.robust", "tli.robust", "srmr", "rmsea.robust"))
# Calcula indices de modificação
im = modindices(fitModeloR, sort. = TRUE, minimum.value = 11)
#IM para os pesos fatoriais
im[im$op == "~~", ]

# Modelo com erros correlacionados
modelo <- '
    Rexperiência =~ ptsd1 + ptsd2 + ptsd3 + ptsd4 + ptsd5
    Evitamento =~ ptsd6 + ptsd7 + ptsd8 + ptsd9 + ptsd10 +ptsd16 + ptsd17
    Agitação =~ ptsd11 +ptsd12 + ptsd13 + ptsd14 + ptsd15
    # Erros correlacionados ou itens correlacionados
    ptsd9	~~	ptsd10
    ptsd10	~~	ptsd11
    '
fitModeloR <- sem(modelo, data = bd, estimator = "MLR")
# Produz output com estimativas estandardizadas
summary(fitModeloR, std = T)
# Seleção de indíces de qualidade de ajustamento
fitMeasures(fitModeloR,
            c("cfi.robust", "tli.robust", "srmr", "rmsea.robust"))
```

### Exemplo 10.4 AFC itens Nominais

```{r Análise Fatorial de itens nominais, echo=TRUE}
# Ler ficheiro de MS-Excel, note o argumento na="999" para identificar os missing values
library(readxl)
bd <-
  read_excel("D:/OneDrive - NOVASBE/AMOS/exercicios/MCSDS.xlsx", na = "999")
# Biblioteca lavaan e modelo
library(lavaan)
modelo <-
  'MCSDS =~  MCSDS_1 + MCSDS_2 + MCSDS_3 + MCSDS_4 + MCSDS_5 + MCSDS_6 + MCSDS_7 + MCSDS_8 + MCSDS_9 + MCSDS_10 + MCSDS_11 + MCSDS_12 + MCSDS_13'

#Obtenção Correlações tetracoricas
#lavCor(bd, ordered=names(bd))

# Ajustamento do modelo
fitModelo <- sem(modelo, data = bd, ordered = T)
# Produz output com estimativas estandardizadas
summary(fitModelo, std = T, fit.m = T)
# Seleção de indíces de qualidade de ajustamento
fitMeasures(fitModelo, c("cfi", "tli", "srmr", "rmsea"))
```

### Exemplo 10.6. Pesquisa de especificação modelo para responder sem preconceito

```{r spec search, eval=FALSE, include=FALSE}

library(haven)
bd <- read_sav("D:/OneDrive - NOVASBE/ReportNumber/Site/aee/bds/BDMotiv.sav")

# Modelo
modelo <- 'MotExt =~ Item1 + Item3 + Item6 + Item7 + Item9
           MotInt =~ Item2 + Item4 + Item5 + Item8 + Item10'
#Fit modelo
library(lavaan)
fit.inic <- sem(modelo, data=bd, std.lv=TRUE)
# Prepara tabela para pesquisa especificação
# Código original de Marcoulides & Falk (2018) https://doi.org/10.1080/10705511.2017.1409074
source("SearchFunctions.R") # Funções para pesquisa de especificação
tabela.pesq<-search.prep(fit.inic,loadings=TRUE,fcov=TRUE,errors=FALSE)
# Faz a pesquisa de especificação minimizando o BIC com 50 interações (nint) e registando 10 para comparação (tabu.size=10)

result.pesq<-tabu.sem(fit.inic,tabela.pesq,obj=BIC,niter=50,tabu.size=10)

#Definir outras funções objetivo
obj.func<-function(x){fitMeasures(x)["rmsea"]}
result.pesq<-tabu.sem(fit.inic,tabela.pesq,obj=obj.func,niter=50,tabu.size=10)
# Visualizar o melhor modelo
summary(result.pesq$best.mod, std=T, fit.measures=T)
#Grafico
library(lavaanPlot)
lavaanPlot(model = result.pesq$best.mod, coefs = TRUE, stand = TRUE, covs = TRUE,graph_options = list(rankdir = "LR"))

```

## Capítulo 11

### Exemplo 11.2 intenção de desistir do curso

```{r exemplo 11.2}

# Importar a base de dados
library(haven)
bd <- read_sav("D:/OneDrive - NOVASBE/AMOS/exercicios/BdUSEI.sav")
# Modelo do USEI
modelo <- '
#Modelo de Medida de 1.ª ordem
  ECP =~ ECP1 + ECP3 + ECP4 + ECP5 + ECP6
  EEM =~ EE14r + EE15 + EE16 + EE17 + EE19
  ECN =~ ECC22 + ECC25 + ECC26 + ECC28 + ECC32
#Modelo de 2.ª ordem
  USEI =~ ECP + EEM + ECN
#Modelo preditivo da itenção de desistir
  Desistir ~ USEI
'
# Ajustamento do Modelo
library(lavaan)
fitModelo <- sem(modelo, data=bd, ordered=T)
# Resultados
summary(fitModelo, std=T, rsq=T)
fitMeasures(fitModelo, c("cfi", "tli", "srmr", "rmsea"))
#Gráfico
library(lavaanPlot)
lavaanPlot(model = fitModelo, coefs = TRUE, stand = TRUE, covs = TRUE,graph_options = list(rankdir = "LR"))

png("graf1.png", width=6 ,height = 3.5, units = "in", res = 300)
x <- 1:5
pm1 <- 1-pnorm(0.280+0.837*x,lower.tail = T)
plot(x, pm1, xlab="USEI", ylab="P(Desistir > 0)", frame=F)
lines(x,pm1)
dev.off()
```

### Exemplo 11.4: Modelo Satisfação Hoteleira

```{r exemplo 13.4, eval=FALSE, include=FALSE}
# Importar a base de dados
# Os missings foram imputadores previamente
library(haven)
bd <-
  read_sav("D:/Dropbox/ReportNumber/Site/form/AEE-VTM/Downloads/SatLealdHoteleiro.sav")
# Estudo da distribuição dos itens
#install.packages("MVN")
library(MVN)
result <- mvn(data = bd [, 2:35], desc = T) # usando as colunas 2 a 35
# Estatisticas descritivas dos itens
result$Descriptives[, c(1, 4:10)]
# Para obter histogramas dos primeiros 6 itens
mvn(bd[, 2:6], univariatePlot = "histogram")

library(lavaan)
# Define o modelo de Medida (Análise factorial Confirmatória) de primeira ordem
SATHOTmm = '
Pessoal =~ EP1 + EP2 + EP4 + EP6
Receção =~ SR1 + SR3 + SR5
BarRest =~ BR1 + BR3 + BR6 + BR7
Segurança =~ S1 + S3 + S4
Alojamento =~ Q1 + Q4 + Q5 + Q7
Satisfação =~ SAT1 + SAT2 + SAT4
Lealdade =~ SAT6+SAT7
'
# Ajustar o modelo com estimdador MLR porque os itens tem 10 pontos
fit.SAT = sem(
  model = SATHOTmm,
  data = bd,
  estimator = "MLR",
  orthogonal = F
)
# Qualidade do Ajustamento e pesos factoriais
summary(fit.SAT, fit.measures = T, standardized = T)
fitmeasures(fit.SAT, c("cfi.robust", "tli.robust", "srmr", "rmsea.robust"))

# Fiabilidade dos fatores
install.packages("semTools")
library (semTools)
reliability(fit.SAT)

# Gráfico do modelo
lavaanPlot(
  model = fit.SAT,
  graph_options = list(layout = "dot", rankdir = "LR"),
  node_options = list(shape = "box", fontname = "Helvetica"),
  edge_options = list(color = "grey"),
  stand = T,
  coefs = TRUE,
  covs = TRUE,
  stars = TRUE
)

## Modelo Completo
SATHOT <-  '
# Modelo de medida 1a ordem
  Pessoal =~ EP1 + EP2 + EP4 + EP6
  Receção =~ SR1 + SR3 + SR5
  BarRest =~ BR1 + BR3 + BR6 + BR7
  Segurança =~ S1 + S3 + S4
  Alojamento =~ Q1 + Q4 + Q5 + Q7
  Satisfação =~ SAT1 + SAT2 + SAT4
  Lealdade =~ SAT6 + SAT7
# Fator de 2a ordem
  QualServ =~ Pessoal + Receção + BarRest + Segurança + Alojamento
# Modelo de Mediação (o argumento := permite calcular estatísticas novas)
    Satisfação ~ a*QualServ     # o "a" é o nome atribuído ao peso estrutural
    Lealdade ~ b*Satisfação
    Lealdade ~ c*QualServ
  # Efeito Mediação (a*b)
    Indireto := a*b
  # Efeito direto (c)
    Direto := c
  # Efeito total (c+ab)
    Total := Indireto + Direto
'
# Ajustamento do modelo
fit.SATHOT = sem(SATHOT, data = bd, estimator = "MLR")
# Sumário do Modelo
summary(
  fit.SATHOT,
  fit.measures = TRUE,
  standardized = TRUE,
  rsquare = TRUE
)

#Fiabilidade dos fatores
reliability(fit.SATHOT)
reliabilityL2(fit.SATHOT, "QualServ")
#OmegaL1 = Proporção da variancia total explicada pelo fator de 2a ordem
#OmegaL2 = Proporção da variÃ¢ncia dos fatores de 1a ordem explicada pelo fator de 2a ordem
#PartialOmegaL1 = ProporaÃ§Ã£o da variancia total explicada pelo fator de 2a ordem depois de descontar a variancia dos erros dos fatores de 1a ordem

# Grafico do modelo
lavaanPlot(
  model = fit.SATHOT,
  graph_options = list(layout = "dot", rankdir = "LR"),
  node_options = list(shape = "box", fontname = "Helvetica"),
  edge_options = list(color = "grey"),
  stand = T,
  coefs = TRUE,
  covs = TRUE,
  stars = TRUE
)

```

### Exemplo 11.6 Interação Burnout e Envolvimento

```{r exemplo 11.6}

# Importar base de dados
bd <- read.csv("D:/OneDrive - NOVASBE/AMOS/exercicios/BEI.csv", sep=";")

# Cria bd.prod com termos da interação com a função indProd da SemTools
# Por defeito centra os itens antes de multiplicar os itens e faz os produtos em pares
library(semTools)
bd.mod=indProd(data = bd,
var1 = c("MBI1","MBI2","MBI3","MBI6","MBI7","MBI8","MBI9","MBI10","MBI11","MBI12", "MBI13"),
var2 = c("SE1", "SE4","SE5","SE7","SE8","SE9","SE10","SE12","SE13","SE14","SE15"))
# Funde as duas bases de dados
bd.final <- merge(bd,bd.mod)

# Definir os Modelos de Medida e modelo causal
modelo <- '
# Modelo Burnout
Exaustão =~  MBI1 + MBI2 + MBI3 + MBI4 + MBI5
Descrença =~ MBI6 + MBI7 + MBI8 + MBI9
Ineficácia =~ MBI10 + MBI11 + MBI12 + MBI13 + MBI14 + MBI15
MBI6 ~~	MBI7 + MBI8   # Correlação entre erros para melhorar o fit
MBI9 ~~ MBI7 + MBI8   # Correlação entre erros para melhorar o fit
Burnout =~ Exaustão + Descrença + Ineficácia
# Modelo Envolvimento
EnvComp =~ SE1 + SE4 + SE5
EnvEmoc =~ SE7 + SE8 + SE9 + SE10
EnvCog =~ SE12 + SE13 + SE14 + SE15
Env =~ EnvComp + EnvEmoc + EnvCog
# Modelo interação
BvsE =~ MBI1.SE1 + MBI2.SE4 + MBI3.SE5 + MBI6.SE7 + MBI7.SE8 + MBI8.SE9 + MBI9.SE10 + MBI10.SE12 + MBI11.SE13 + MBI12.SE14 + MBI13.SE15
#Modelo Causal
DI ~ Burnout + Env + BvsE
'
# Ajuste do modelo
fitModelo <-  sem(modelo, data=bd.final, ordered="DI")
summary(fitModelo, std=T, rsquare=T)
fitmeasures(fitModelo, c("CFI", "TLI","SRMR", "RMSEA"))


#Ilustra a interação
int <- probe2WayMC(fitModelo,nameX=c("Burnout", "Env", "BvsE"),nameY="DI", modVar="Env", valProbe = c(1,3,5))
int
plotProbe(int,xlim = c(0,6), xlab = "Burnout", ylab = "Desistir [0-Não; 4-Com frequência]", legend=T, frame=F)
#legend(0,3, legend=c("Env = 1", "Env = 3", "Env = 5"), lty = 1:3, box.lty = 0)

```

## Capítulo 12

### Exemplo 12.2 - Invariância de género

```{r exemplo 12.2}

library(haven)
bd <- read_sav("D:/OneDrive - NOVASBE/AMOS/exercicios/PEVI.sav")
bd$SexR <- as.numeric(bd$SexR)    # converte a variável SexR de chr para numeric
# Modelo pentafatorial do PEVI
modelo <- '
NUT =~ PEVI1 + PEVI2 + PEVI3        # Nutrição
ATF =~ PEVI4 + PEVI5 + PEVI6        # Actividade Fisica
CPR =~ PEVI7 + PEVI8 + PEVI9        # Comportamentos Preventivos
RSO =~ PEVI10+ PEVI11 + PEVI12      # Relações Sociais
GES =~ PEVI13 + PEVI14 + PEVI15     # Gestão de Stress
'
#Ajusta o modelo
library(lavaan)
fit <-  sem(modelo, data=bd, ordered = T)
summary(fit, fit.measures=T, standardized=T)
fitmeasures(fit, c("CFI", "TLI","SRMR", "RMSEA"))

# Análise de Invariância de Genero (ver texto para detalhes)
library (semTools)

# Análise antiga

#If strict = FALSE, the following four models are tested in order:
#Model 1: configural invariance. The same factor structure is imposed on all groups.
#Model 2: weak invariance. The factor loadings are constrained to be equal across groups.
#Model 3: strong invariance. The factor loadings and intercepts are constrained to be equal across groups.
#Model 4: The factor loadings, intercepts and means are constrained to be equal across groups.

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

measurementInvarianceCat(model = modelo, data = bd, group = "SexR",ordered=it, strict = T, fit.measures = "default", missing="listwise")

test.seq <- c("thresholds", "loadings", "intercepts","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 = modelo,            # modelo inicial sem restrições, para inv. config
    data = bd,
    ordered = TRUE,                       # indica que todos os itens são ordinais. 
    parameterization = "theta",           # identificação das variâncias dos residuos
    ID.cat = "Wu.Estabrook.2016",         # identificação dos thresholds
    ID.fac = "std.lv",                    # identifica as variáveis lantentess com M=0, SD=1
    group = "SexR",                       # variável que define os grupos
    group.equal = group.equal,            # parâmetros a fixar: thresholds, loadings,...
    return.fit = TRUE)
}
tabela <- compareFit(meq.list)                    # compara os modelos ajustados anteriormente
print(x=tabela)

### alternativa

fit_config <- sem(modelo, data=bd, group="SexR", ordered = T, std.lv=T)
fit_thresholds <- sem(modelo, data=bd, group="SexR", ordered=T, std.lv=T, group.equal = c("thresholds"))
fit_metric <- sem(modelo, data=bd, group="SexR", ordered=T, std.lv=T, group.equal = c("thresholds","loadings"))
fit_scalar <- sem(modelo, data=bd, group="SexR", ordered=T, std.lv=T, group.equal = c("thresholds","loadings", "intercepts"))
fit_means <- sem(modelo, data=bd, group="SexR", ordered=T, std.lv=T, group.equal = c("thresholds","loadings", "intercepts", "means"))

summary.FitDiff<-function(object){
          print(object@nested)
          return(object@fit)
}

s <- compareFit(fit_config, fit_thresholds, fit_metric, fit_scalar, fit_means)
summary(s)

```

### Exemplo 12. Relação entre a Satisfação e Burnout

```{r exemplo 12.4}

#bd <- sijs[, c(3,11:15,32:46,47:55)]
#saveRDS(bd,"D:/OneDrive - NOVASBE/AMOS/exercicios/SIJS.rds")
# Lê a base de dados no formato nativo do R (rds)
bd <- readRDS("D:/OneDrive - NOVASBE/AMOS/exercicios/SIJS.rds")

# com variancias negativas

Modelo <- '
#Satisfação no Trabalho (SIJS)
SatTrb =~ SIJS1 + SIJS2 + SIJS3 + SIJS4 + SIJS5
#Envolvimento (UWES3)
Env =~ UWES1 + UWES3 + UWES6
#Burnout (OLBI)
Distanc =~ OLBI1 + OLBI3 + OLBI6 + OLBI7 + OLBI9 + OLBI11 + OLBI15
Exaust =~ OLBI2 + OLBI4 + OLBI5 + OLBI8 + OLBI10 + OLBI12 + OLBI14 + OLBI16
Burnout =~ Distanc + Exaust '
library(lavaan)
fit <- sem(Modelo, data=bd, ordered=T)

# Encontrar variâncias estandardizadas dos erros (matriz psi)
round(diag(lavInspect(fit, "std")$psi),3)


#Modelos
Modelo <- '
#Satisfação no Trabalho
SatTrb =~ SIJS1 + SIJS2 + SIJS3 + SIJS4 + SIJS5
#Envolvimento
Env =~ UWES1 + UWES3 + UWES6
#Burnout
Distanc =~ OLBI1 + OLBI3 + OLBI6 + OLBI7 + OLBI9 + OLBI11 +  OLBI15
Exaust =~ OLBI2 + OLBI4 + OLBI5 + OLBI8 + OLBI10 + OLBI12 + OLBI14 + OLBI16
Burnout =~  Exaust + Distanc
Distanc ~~ 0.01*Distanc
#Regressão
Burnout ~ SatTrb + Env'
# Ajuste modelo
library(lavaan)
fit <- sem(Modelo, data=bd, ordered=T, group="Country")
summary(fit, standardized=T, rsquare=T, fit.measures=T)
fitmeasures(fit, c("cfi.scaled", "tli.scaled","srmr", "rmsea.scaled"))
#standardizedSolution(fit)

parameterEstimates(fit)[fit@ParTable$op == "~",]


#Comparação declives

Z=(-0.718+.456)/sqrt(0.075^2+0.038^2)
2*pnorm(abs(Z), lower.tail = F)

# Encontrar variâncias estandardizadas dos erros (matriz psi)
round(diag(lavInspect(fit, "std")$psi),3)

#Análise invariância do modelo
library(semTools)
test.seq <- c("loadings", "intercepts", "regressions")
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 = Modelo,        
    data = bd,
    estimator = "MLR",
    ID.fac = "ul",
    group = "Country",                    
    group.equal = group.equal,           
    return.fit = TRUE)
}
compareFit(meq.list) 

#Para quando dá o erro S3 class...
 summary.FitDiff<-function(object){
          print(object@nested)
          return(object@fit)}
a<-semTools::compareFit(meq.list) 
summary(a)
```

### Exemplo 12.5 comparação de Médias de variáveis manifestas

```{r exemplo 12.5, eval=T, include=T}

bd <- data.frame(
          Group = c(1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,
                 3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,2L,
                 2L,2L,2L,2L,2L,2L,2L,2L,2L,2L,2L,2L,2L,2L,2L,
                 2L),
          HIP_D = c(2763,3241.5,2740.5,3318,3501,3634.5,
                 4324.5,2970,3168,3307.5,3145.5,3639,3325.5,3295.5,3066,
                 3711,2488.5,3333,3138,2752.5,3096,2785.5,3204,3781.5,
                 3207,2742,3082.5,2665.5,2737.5,2664,2743.5,2695.5,3244.5,
                 3024,2422.5,1605,2943,2314.5,2971.5,3187,2673,3024,
                 2992.5),
          HIP_E = c(2505,3292.5,2646.7,3457.5,3505.5,
                 3628.5,4152,2686.5,2943,2913,3283.5,3546,3172.5,3090,3211.5,
                 3880.5,1725,3147,3055.5,2776.5,2904,2785.5,3159,
                 3418.5,3240,2334,2634,2370,2736,2503.5,2325,2409,3057,
                 2692.5,2443.5,1645.5,2922,1818,2679,2742,2465.5,2911.5,
                 2890.5)
)


library(haven)
bd <- read_sav("D:/OneDrive - NOVASBE/AMOS/exercicios/MCI_NMR.sav")
bd$Diagnostico <- factor(bd$Diagnostico, levels=c(1,2,3), labels=c("Controlo", "PreMCI", "MCI"))
#É necessário rescalonar as variáveis para produzir estimativas corretas no lavaan
bd$HIP_D <- bd$HIP_D/100
bd$HIP_E <- bd$HIP_E/100

#correlation::correlation(bd)
library(lavaan)
#Modelo para comparar variâncias
# Variâncias livres
modVL <- '
HIP_D ~~ HIP_D
HIP_E ~~ HIP_E
HIP_D ~~ HIP_E
'
#Para usar o estimador não enviesado (divide por n-1) da variância usar likelihood = "wishart"
fitVL <- sem(modVL, data=bd, group="Diagnostico")
summary(fitVL, standardized=T)
#fitMeasures(fitVL)
#Variâncias/covariâncias homogéneas
modVC <- '
HIP_D ~~ v*HIP_D
HIP_E ~~ v*HIP_E
HIP_D ~~ cov*HIP_E
'
fitVC <- sem(modVC, data=bd,group="Diagnostico")
summary(fitVC, standardize=T)

# Teste do Qui^2 para a igual de de variancias e covariancias
lavTestLRT(fitVL, fitVC)

#constrained

#Teste para as médias
modVCMC <- '
#Var e Cov iguais
HIP_D ~~ v*HIP_D
HIP_E ~~ v*HIP_E
HIP_D ~~ cov*HIP_E
#Médias dos 3 grupos iguais
HIP_D ~ m*1               # o  ~1 ajusta um modelo só com intercepto
HIP_E ~ m*1
'
fitVCMC <- sem(modVCMC, data=bd, group="Diagnostico")
summary(fitVCMC, standardize=T)
lavTestLRT(fitVC, fitVCMC)

#Compara grupos na HIP_D
# Variâncias homogéneas e médias livres
modHPD_ML <- '
HIP_D ~~ v*HIP_D
HIP_D ~ 1' 
fitHPD_ML <- sem(modHPD_ML, data=bd, group="Diagnostico")
summary(fitHPD_ML, standardize=T)
#Variâncias homogéneas e médias iguais
modHPD_MC <- '
HIP_D ~~ v*HIP_D
HIP_D ~ m*1' 
fitHPD_MC <- sem(modHPD_MC, data=bd, group="Diagnostico")
summary(fitHPD_MC, standardize=T)
# Teste diferença de Qui^2
lavTestLRT(fitHPD_MC, fitHPD_ML)

#Compara grupo 1 com grupo 2
modHPD_MC12 <- '
HIP_D ~~ v*HIP_D
HIP_D ~ c("m", "m", NA)*1' 
fitHPD_MC12 <- sem(modHPD_MC12, data=bd, group="Diagnostico")
lavTestLRT(fitHPD_MC12, fitHPD_ML)

#Compara grupo 1 com grupo 3
modHPD_MC13 <- '
HIP_D ~~ v*HIP_D
HIP_D ~ c("m", NA, "m")*1' 
fitHPD_MC13 <- sem(modHPD_MC13, data=bd, group="Diagnostico")
lavTestLRT(fitHPD_MC13, fitHPD_ML)

#Compara grupo 2 com grupo 3
modHPD_MC23 <- '
HIP_D ~~ v*HIP_D
HIP_D ~ c(NA, "m", "m")*1' 
fitHPD_MC13 <- sem(modHPD_MC23, data=bd, group="Diagnostico")
lavTestLRT(fitHPD_MC13, fitHPD_ML)

#Teste de Tukey sem variâncias homogéneas
library(semTools)
tukeySEM(m1=31.710, m2=27.641, var1=12.951, var2=14.527, n1=15, n2=17,ng=3)
tukeySEM(m1=32.831, m2=27.641, var1=17.765, var2=14.527, n1=43, n2=43, ng=3)


```

### Exemplo 12.7 Comparação medias OTOCS

```{r}

#Lê o ficheirod e dados

bdOriginal <-  readRDS("D:/OneDrive - NOVASBE/AMOS/exercicios/df.OTOCS5.rds")

bd <- bdOriginal[, c("Country", "OTCS.OTCS1.", "OTCS.OTCS2.","OTCS.OTCS3.","OTCS.OTCS4.", "OTCS.OTCS5.")]

colnames(bd) <- c("País", "OTOCS1","OTOCS2","OTOCS3","OTOCS4","OTOCS5")

saveRDS(bd,"D:/OneDrive - NOVASBE/AMOS/exercicios/otocs.rds")


#Lê ficheiro de dados em formato Rds
bd <-  readRDS("D:/OneDrive - NOVASBE/AMOS/exercicios/otocs.rds")

#Define o modelo com dois erros correlacionados
library(lavaan)
modelo <- 'OTOCS =~ OTOCS1 + OTOCS2 + OTOCS3 + OTOCS5
          OTOCS1 ~~ OTOCS2'
#ajusta o modelo
fit <- sem(model=modelo, data=bd, std.lv=T, ordered=T)
summary(fit, fit.measures=T, standardized=T)
fitmeasures(fit, c("cfi.scaled", "tli.scaled","srmr", "rmsea.scaled"))
modificationindices(fit, sort=T)

#Análise de  invariância
library(semTools)
test.seq <- c("thresholds","loadings","intercepts","means")
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 = modelo,        
    data = bd,
    ordered=T,
    group = "País",                    
    group.equal = group.equal,           
    return.fit = TRUE)
}
print(compareFit(meq.list))

# Comparação de Médias

fitModML <-  sem(modelo, data=bd, group="País", ordered=T, group.equal=c("thresholds", "loadings", "intercepts"))  #Estima as médias livres
summary(fitModML, standardized=T, fit.measures=T)
fitModMC <-  sem(modelo, data=bd, group="País", ordered=T, group.equal=c("thresholds", "loadings", "intercepts", "means"))  #Estima as médias livres
summary(fitModMC, standardized=T, fit.measures=T)
#Testa diferenças dos dois modelos (ML vs MC)
lavTestLRT(fitModMC, fitModML)

#Modelo MIMIC
modeloM <- 'OTOCS =~ OTOCS1 + OTOCS2 + OTOCS3 + OTOCS5
          OTOCS1 ~~ OTOCS2
OTOCS ~ País  # Regressão Probit da OTOCS em função do País'
fitMM <- sem(modeloM, data=bd, ordered = T, std.lv=T)
summary(fitMM, standardized=T, rsquare=T)

#Recode dummies
bd <-  readRDS("D:/OneDrive - NOVASBE/AMOS/exercicios/otocs.rds")
bd$País <-  as.numeric(bd$País)
bd$País[bd$País==1] <- 0
bd$País[bd$País==2] <- 1
table(bd$País)


```

## Capítulo 13

### Exemplo 13.3 Trajetórias de Investimento

```{r exemplo 13.3, eval=FALSE, include=FALSE}

# Lê base de Dados
bd <-  read.csv2("D:/OneDrive - NOVASBE/AMOS/exercicios/ivp_br.csv", dec=".")

# Boxplots do IGS

png("Plot13_3_1.png", width=10 ,height = 5.5, units = "in", res = 300)
par(mfrow=c(1,2)) #Define a area dos graficos num array 1l*2C
 boxplot(bd[,7:11], xlab="Ano", ylab="Indice Gasto Social", frame=F)
 boxplot(bd[,17:21], xlab="Ano", ylab="Indice Des. Humano", frame=F)
dev.off()

# Define os modelos de crescimento latente
modCLDH <- '
# Gasto Social
interc.GS =~ 1*igs1992 + 1*igs1996 + 1*igs2000 + 1*igs2004 + 1*igs2008
decliv.GS =~ 0*igs1992 + igs1996 + igs2000 + 1*igs2004 + igs2008
# Desenvolvimento Humano
interc.DH =~ 1*DH.1992 + 1*DH.1996 + 1*DH.2000 + 1*DH.2004 + 1*DH.2008
decliv.DH =~ 0*DH.1992 + DH.1996 + DH.2000 + DH.2004 + 1*DH.2008
# Regressão do crescimento do DH na DS
decliv.DH ~ decliv.GS
'
# Ajustamento do Modelo
library (lavaan)
fitmodCLDH <- growth(model=modCLDH, data=bd)
summary(fitmodCLDH, standardize=T, fit.measures=T)

# Gráfico do modelo
library(semPlot)
png("Plot13_14_2.png", width=3 ,height = 3.5, units = "in", res = 300)
semPaths(modCLDH, rotation=2, "std", style="lisrel",fixedStyle = c("black",1), freeStyle = c("black",1), asize=2, esize=1,edge.label.cex=1.5)
dev.off()



```

### Exemplo 13.4 Evolução da democracia liberal

```{r ex 13.4, eval=FALSE, include=FALSE}

# Ler base de dados
bd <- read.delim(file="D:/OneDrive - NOVASBE/AMOS/exercicios/ld.txt", na.strings=".")
#reduzir a escala das libdem para eliminar problemas com estimação das variancias
bd$libd80 = bd$libdem80 / 10
bd$libd83 = bd$libdem83 / 10
bd$libd86 = bd$libdem86 / 10

library(MVN)
result <- mvn(data = bd [,-10], desc=T)
# Estatisticas descritivas dos itens
result$Descriptives[,c(1, 4:10)]
# Testes à normalidade univariada com o teste de Shapiro-Wilk
result$univariateNormality
# Teste à normalidade multivariada com o teste de Mardia
result$multivariateNormality
# Para obter histogramas dos primeiros 6 itens
mvn(bd[,-10], univariatePlot = "histogram")

library(lavaan)
# Modelo de Perceção Democracia Liberal
ModDL <- '
  DL80 =~ libd80 + gdemoc80 + vanhan80 
  DL83 =~ libd83 + gdemoc83 + vanhan83 
  DL86 =~ libd86 + gdemoc86 + vanhan86
  #Correlaciona os erros das medidas longitudinais
  libd80 ~~ libd83 + libd86
  libd83 ~~ libd86
  gdemoc80 ~~ gdemoc83 + gdemoc86
  gdemoc83 ~~ gdemoc86
  vanhan80 ~~ vanhan83 + vanhan86
  vanhan83 ~~ vanhan86
'
fitModDL <- sem(model=ModDL, data=bd, estimator="MLR")
summary(fitModDL, standardize=T, fit.measures=T)


library(semPlot)
png("Plot13_14.png", width=3 ,height = 3, units = "in", res = 300)
semPaths(fitModDL, rotation=2, "std", style="lisrel",fixedStyle = c("black",1), freeStyle = c("black",1), asize=2, esize=1,edge.label.cex=1.5)
dev.off()

#Modelo de Crescimento Latente
ModCL <- '
  DL80 =~ libd80 + gdemoc80 + vanhan80 
  DL83 =~ libd83 + gdemoc83 + vanhan83 
  DL86 =~ libd86 + gdemoc86 + vanhan86
  #Correlaciona os erros das medidas longitudinais
  libd80 ~~ libd83 + libd86
  libd83 ~~ libd86
  gdemoc80 ~~ gdemoc83 + gdemoc86
  gdemoc83 ~~ gdemoc86
  vanhan80 ~~ vanhan83 + vanhan86
  vanhan83 ~~ vanhan86
  #Fixa intercept dos itens
  libd80 ~ 0*1
  libd83 ~ 0*1
  libd86 ~ 0*1
  gdemoc80 ~ 0*1
  gdemoc83 ~ 0*1
  gdemoc86 ~ 0*1
  vanhan80 ~ 0*1
  vanhan83 ~ 0*1
  vanhan86 ~ 0*1
  #Fixa intercepts das v. latente
  DL80 ~ 0*1
  DL83 ~ 0*1
  DL86 ~ 0*1
  #Fixa variancias de erros negativas
  libd80 ~~ 0.1*libd80
  libd83 ~~ 0.1*libd83
  libd86 ~~ 0.1*libd86
  DL80 ~~ 0.1*DL80
  #Intercepto
  intercpt =~ 1*DL80 + 1*DL83 + 1*DL86
  #Declive
  declive = ~ 0*DL80 + 1*DL83 + 2*DL86
  '
fitModCL <- growth(model=ModCL, data=bd, missing="fiml")
summary(fitModCL, standardize=T, fit.measures=T)
resid(fitModCL)  # to investigate the negative variances and convariances residuals
lavInspect(fitModCL, "theta")

png("Plot13_14_2.png", width=3 ,height = 3.5, units = "in", res = 300)
semPaths(fitModCL, rotation=2, "std", style="lisrel",fixedStyle = c("black",1), freeStyle = c("black",1), asize=2, esize=1,edge.label.cex=1.5)
dev.off()

fitModCL <- growth(model=ModCL, data=bd, missing="fiml", group="Continente")
summary(fitModCL, standardize=T, fit.measures=T)

```

### Exemplo 13.6 Invariância longitudinal

```{r ex 13.6}
bd <- read.delim(file="D:/OneDrive - NOVASBE/AMOS/exercicios/ld.txt", na.strings=".")
#reduzir a escala das libdem para eliminar problemas com estimação das variancias
bd$libd80 = bd$libdem80 / 10
bd$libd83 = bd$libdem83 / 10
bd$libd86 = bd$libdem86 / 10

library(MVN)
result <- mvn(data = bd [,-10], desc=T)
# Estatisticas descritivas dos itens
result$Descriptives[,c(1, 4:10)]
# Testes à normalidade univariada com o teste de Shapiro-Wilk
result$univariateNormality
# Teste à normalidade multivariada com o teste de Mardia
result$multivariateNormality
# Para obter histogramas dos primeiros 6 itens
mvn(bd[,-10], univariatePlot = "histogram")

library(lavaan)
# Modelo de Perceção Democracia Liberal
ModDL <- '
  DL80 =~ libd80 + gdemoc80 + vanhan80 
  DL83 =~ libd83 + gdemoc83 + vanhan83 
  DL86 =~ libd86 + gdemoc86 + vanhan86
  #Correlaciona os erros das medidas longitudinais
  libd80 ~~ libd83 + libd86
  libd83 ~~ libd86
  gdemoc80 ~~ gdemoc83 + gdemoc86
  gdemoc83 ~~ gdemoc86
  vanhan80 ~~ vanhan83 + vanhan86
  vanhan83 ~~ vanhan86
'
library(lavaan)
fitModDL <- sem(model=ModDL, data=bd, estimator="MLR")
summary(fitModDL, standardize=T, fit.measures=T)

# Análise de Invariância Longitudinal
library(semTools)
NomeFactLong <- list(FU = c("DL80","DL83", "DL86"))
#Compara vários modelos aninhados
test.seq <- c("loadings","intercepts","means","residuals")
meq.list <- list()
for (i in 0:length(test.seq)) {
  if (i == 0L) {
    meq.label <- "configural"
    long.equal <- ""
  } else {
    meq.label <- test.seq[i]
    long.equal <- test.seq[1:i]
  }
  meq.list[[meq.label]] <- measEq.syntax(configural.model = ModDL,
                          data = bd,parameterization = "theta",
                          ID.fac = "std.lv",ID.cat = "Wu",
                          longFacNames = NomeFactLong, long.equal = long.equal,
                          return.fit = TRUE)
}

tabInvLon <- compareFit(meq.list)
print(tabInvLon)
```
