# Packages, data, working directory #######################################

setwd("C:/FSS/SEM")
library(psych)
library(lavaan)
library(semTools)
library(semPlot)
library(corrplot)
source("https://raw.githubusercontent.com/hynekcigler/Functions/master/extract.R")

Pedhazurdata <- read.csv("pedhazur_latent2.csv", header = T)
head(Pedhazurdata)

# Prakticke ukazky ####################################################

# Model 1 - IQ a nAch  ####################################################

modelIQ <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
IQ ~~ 1*IQ
"

cfamodelIQ <- cfa(modelIQ, data = Pedhazurdata, estimator = "ML")
summary(cfamodelIQ, fit=T, std=T)
semtext(cfamodelIQ, type = "standard")
semPaths(cfamodelIQ, what="path", whatLabels = "std", layout = "tree", exoVar=F)

modelnAch <- "
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10
nAch ~~ 1*nAch
"

cfamodelnAch <- cfa(modelnAch, data = Pedhazurdata, estimator = "ML")
summary(cfamodelnAch, fit=T, std=T)
semtext(cfamodelnAch, type = "standard")
semPaths(cfamodelnAch, what="path", whatLabels = "std", layout = "tree", exoVar=F)

model1 <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 

nAch ~ IQ

nAch ~~ 1*nAch
IQ ~~ 1*IQ
"

modelsem1 <- cfa(model1, data = Pedhazurdata, estimator = "ML")
summary(modelsem1, fit=T, std=T, rsquare = T)
semtext(modelsem1, type = "standard")

semPaths(modelsem1, what="path", whatLabels = "std", layout = "tree", exoVar=F)

lavInspect(modelsem1)
lavInspect(modelsem1, what = "est")
lavInspect(modelsem1, what = "std")
corrplot(residuals(modelsem1, type = "cor")$cov)

# Model 2 - IQ, SES a nAch  ################################################


model2 <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

nAch ~ a*IQ
IQ ~ b*SES
nAch ~ c*SES

ab := a*b

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"

modelsem2 <- cfa(model2, data = Pedhazurdata, estimator = "ML")
summary(modelsem2, fit=T, std=T, rsquare = T)
semtext(modelsem2, type = "standard")

semPaths(modelsem2, what="path", whatLabels = "std", layout = "tree", exoVar=F)

lavInspect(modelsem2)
lavInspect(modelsem2, what = "est")
lavInspect(modelsem2, what = "std")
corrplot(residuals(modelsem2, type = "cor")$cov)

# Znovu s par ukazkami syntaxe

model2_loadings <- "
IQ =~ NA*IQ1 + l*IQ1 + l*IQ2 + l*IQ3 + l*IQ4 + l*IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

nAch ~ a*IQ
IQ ~ b*SES
nAch ~ c*SES

ab := a*b

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"
modelsem2_loadings <- cfa(model2_loadings, data = Pedhazurdata, estimator = "ML")
summary(modelsem2_loadings, fit=T, std=T, rsquare = T)
semtext(modelsem2_loadings, type = "standard")

model2_loadings_SES <- "
IQ =~ NA*IQ1 + l*IQ1 + l*IQ2 + l*IQ3 + l*IQ4 + l*IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

nAch ~ a*IQ
IQ ~ b*SES
nAch ~ 0*SES

ab := a*b

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"
modelsem2_loadings_SES <- cfa(model2_loadings_SES, data = Pedhazurdata, estimator = "ML")
summary(modelsem2_loadings_SES, fit=T, std=T, rsquare = T)
semtext(modelsem2_loadings_SES, type = "standard")

# Srovnani dvou modelu pomoci LRT - s predpokladem primeho efektu SES na nAch a bez tohoto predpokladu

lavTestLRT(modelsem2_loadings, modelsem2_loadings_SES)

# ^--- zda se, ze model s primym efektem SES na nAch sedi lepe na data

# Model 3 - IQ, SES, nAch, GPA  ################################################

model3 <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

GPA ~ IQ + nAch + SES

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"

modelsem3 <- cfa(model3, data = Pedhazurdata, estimator = "ML")
summary(modelsem3, fit=T, std=T, rsquare = T)
semtext(modelsem3, type = "standard")

semPaths(modelsem3, what="path", whatLabels = "std", layout = "tree", exoVar=F)

lavInspect(modelsem3)
lavInspect(modelsem3, what = "est")
lavInspect(modelsem3, what = "std")
corrplot(residuals(modelsem3, type = "cor")$cov)

# Model 4 - Pribyva vysvetleny rozptyl?   ################################################

model4a <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

GPA ~ IQ + 0*nAch + 0*SES

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"
modelsem4a <- cfa(model4a, data = Pedhazurdata, estimator = "ML")
model4b <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

GPA ~ IQ + nAch + 0*SES

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"
modelsem4b <- cfa(model4b, data = Pedhazurdata, estimator = "ML")
model4c <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

GPA ~ IQ + nAch + SES

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"
modelsem4c <- cfa(model4c, data = Pedhazurdata, estimator = "ML")

lavTestLRT(modelsem4a, modelsem4b)
lavTestLRT(modelsem4b, modelsem4c)

# Model 5 - Moderace - GPA  ################################################

# z GPA uděláme indikátor skupiny

table(Pedhazurdata$GPA)
Pedhazurdata$GPA_group <- ifelse(Pedhazurdata$GPA > 2, "Vysoke GPA", "Nizke GPA")

model5 <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

IQ ~ c(a,a)*nAch + SES

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"

modelsem5 <- cfa(model5, data = Pedhazurdata, estimator = "ML", group = "GPA_group", group.equal = c("loadings","lv.covariances","regressions"))
summary(modelsem5, fit=T, std=T, rsquare = T)
semtext(modelsem5, type = "standard")

model5b <- "
IQ =~ NA*IQ1 + IQ2 + IQ3 + IQ4 + IQ5
nAch =~ NA*nAch1 + nAch2 + nAch3 + nAch4 + nAch5 + nAch6 + nAch7 + nAch8 + nAch9 + nAch10 
SES =~ NA*SES1 + NA*SES2 + NA*SES3

IQ ~ c(a,b)*nAch + SES

nAch ~~ 1*nAch
IQ ~~ 1*IQ
SES ~~ 1*SES
"
modelsem5b <- cfa(model5b, data = Pedhazurdata, estimator = "ML", group = "GPA_group", group.equal = c("loadings","lv.covariances","regressions"))
summary(modelsem5b, fit=T, std=T, rsquare = T)
semtext(modelsem5b, type = "standard")

lavTestLRT(modelsem5, modelsem5b)