#################################################################
###### Data

read.csv(file = "C:/Adam/datamatrix.csv", header = T)

# Promenné:
#
# Word Meaning
# Sentence Comprehension
# Odd Words
# Mixed Arithmetic
# Remainders
# Missing Numbers
# Gloves
# Boots
# Hatchets

# N = 696

# Ulozíme data do objektu
cormatrix <- read.csv(file = "C:/Adam/datamatrix.csv", header = T)

# Prevedeme objekt na matici
cormatrix <- as.matrix(cormatrix)

#################################################################
###### Packages

library(lavaan)
library(semTools)
library(corrplot)
library(semPlot)

#################################################################
###### Specifikace modelu

# Jednofaktorový model

model1F <- '

### Faktorové náboje
# NA* u první promenné u kazdého faktoru sdeluje lavaanu, ze si neprejeme urcit rozptyl faktoru fixací prvního náboje

F1 =~ NA*WrdMean + SntComp + OddWrds + MxdArit + Remndrs + MissNum + Gloves + Boots + Hatchts  # Vsechny manifestní promenné na jednom faktoru

### Rozptyly latentních promenných a reziduální rozptyly

F1 ~~ 1*F1  # Nastavíme rozptyl prvního faktoru na 1
'

# Defaulty parametru:
# 
# Matice náboju (lambda) je nulová (krome prvního náboje u kazdého faktoru, který je 1), v syntaxi prepisujeme
# Matice faktorových korelací (Phí, lavaan jí ríká Psí) obsahuje volne odhadnutelné parametry, v syntaxi muzeme prepsat
# Matice reziduálních korelací (D-psí, lavaan jí ríká Théta) obsahuje volne odhadnutelné parametry pouze na diagonále, v syntaxi muzeme prepsat


## Faktorové náboje:
# V syntaxi urcujeme, které parametry se volne odhadnou nebo zafixují/omezí jinak nez na nulu
# Pocet faktoru urcujeme tím, ze je jmenujeme a rovnou urcíme, u kterých promenných chceme odhadnout faktorové náboje
# Polozky, které faktor "sytí", uvádíme za operátor '=~', a oddelujeme pomocí "+"

# Kvuli implicitní fixaci prvního náboje na 1 je toto:
# F1 =~ WrdMean + SntComp + OddWrds
# 
# ...ve skutecnosti tímto:
# F1 =~ 1*WrdMean + SntComp + OddWrds 

# ...cemuz zamezíme následovne:
# F1 =~ NA*WrdMean + SntComp + OddWrds
# ...tím "uvolníme" první náboj.

## Faktorové rozptyly, korelace a kovariance:
# V syntaxi urcujeme, které parametry se zafixují na nulu/omezí jinak nez na nulu
# Uvádíme názvy faktoru, kterým jsou v syntaxi prirazeny nejaké faktorové náboje 
# Rozptyly, korelace i kovariance urcujeme operátorem '~~':
#
# Rozptyl:
# F1 ~~ parametr*F1
#
# Korelace / kovariance:
# F1 ~~ F2
#
# F1 ~~ F2
# F2 ~~ F3
# F1 ~~ F3
#
# F1 ~~ F2 + F3

## Korelace / kovariance reziduí:
# V syntaxi urcujeme, které parametry mimo diagonálu se volne odhadnou 
# Muzeme menit i povahu parametru na diagonále, ale to je pomerne vzácná situace
# Uvádíme názvy manifestních promenných, korelace / kovariance urcujeme operátorem '~~', stejne jako faktorové korelace / kovariance:
#
# Korelace / kovariance:
# WrdMean ~~ SntComp
#
# WrdMean ~~ SntComp
# SntComp ~~ OddWrds
# Wrdmean ~~ OddWrds
#
# WrdMean ~~ SntComp + OddWrds

model1F <- '

### Faktorové náboje
# NA* u první promenné u kazdého faktoru sdeluje lavaanu, ze si neprejeme urcit rozptyl faktoru fixací prvního náboje

F1 =~ NA*WrdMean + SntComp + OddWrds + MxdArit + Remndrs + MissNum + Gloves + Boots + Hatchts  # Vsechny manifestní promenné na jednom faktoru

### Rozptyly latentních promenných a reziduální rozptyly

F1 ~~ 1*F1  # Nastavíme rozptyl prvního faktoru na 1
'

#################################################################
###### Odhad parametru

fit1F <- sem(model = model1F, sample.cov = cormatrix, sample.nobs = 696, estimator = "ML")

# Výpis matic s ocíslovanými parametry

lavInspect(fit1F)

# Jednotlivé matice

lavInspect(fit1F)$lambda
lavInspect(fit1F)$psi
lavInspect(fit1F)$theta

# Výpis matic s odhadnutými parametry

lavInspect(fit1F, what = "est")

lavInspect(fit1F, what = "est")$lambda
lavInspect(fit1F, what = "est")$psi
lavInspect(fit1F, what = "est")$theta

# Výpis parametru modelu a dalsích informací o modelu:

summary(fit1F, standardized = T)

# Reziduální korelacní matice

residuals(fit1F)

# ...muzeme si ji vizualizovat:

corrplot(residuals(fit1F)$cov)

# Výpis indexu shody modelu s daty a dalsích informací o modelu:

fitmeasures(fit1F, fit.measures = c("rmsea", 
                                     "tli", "df","npar",
                                     "chisq",
                                     "rmsea.ci.lower","rmsea.ci.upper", "cfi", "srmr"))

semPaths(fit1F, what="path", whatLabels = "std", layout = "tree")

# Trojfaktorový model s nekorelovanými faktory

### Hypotéza o existenci trí faktoru
### Verbální schopnost (WrdMean, SntComp a OddWrds)
### Numerická schopnost (MxdArit, Remndrs a MissNum)
### Prostorová schopnost (Gloves, Boots a Hatchts)

### V tuto chvíli predpokladáme, ze tyto tri schopnosti spolu nekorelují


model3F <- '

### Faktorové náboje

Verbal =~ NA*WrdMean + SntComp + OddWrds
Numerical =~ NA*MxdArit + Remndrs + MissNum
Spatial =~ NA*Gloves + Boots + Hatchts

### Rozptyly latentních promenných 

Verbal ~~ 1*Verbal
Numerical ~~ 1*Numerical
Spatial ~~ 1*Spatial

### Korelace / kovariance latentních promenných 

Verbal ~~ 0*Numerical
Verbal ~~ 0*Spatial
Numerical ~~ 0*Spatial
'

fit3F <- sem(model3F, sample.cov = cormatrix, sample.nobs = 696, estimator = "ML")

summary(fit3F)

lavInspect(fit3F)$lambda
lavInspect(fit3F)$psi
lavInspect(fit3F)$theta

lavInspect(fit3F, what = "est")$lambda
lavInspect(fit3F, what = "est")$psi
lavInspect(fit3F, what = "est")$theta

residuals(fit3F)
corrplot(residuals(fit3F)$cov)

fitmeasures(fit3F, fit.measures = c("rmsea", 
                                    "tli", "df","npar",
                                    "chisq",
                                    "rmsea.ci.lower","rmsea.ci.upper", "cfi", "srmr"))

# Trojfaktorový model s korelovanými faktory

### Hypotéza o existenci trí faktoru
### Verbální schopnost (WrdMean, SntComp a OddWrds)
### Numerická schopnost (MxdArit, Remndrs a MissNum)
### Prostorová schopnost (Gloves, Boots a Hatchts)

### V tuto chvíli predpokladáme, ze tyto tri schopnosti spolu korelují

model3Fc <- '

### Faktorové náboje

Verbal =~ NA*WrdMean + SntComp + OddWrds
Numerical =~ NA*MxdArit + Remndrs + MissNum
Spatial =~ NA*Gloves + Boots + Hatchts

### Rozptyly latentních promenných 

Verbal ~~ 1*Verbal
Numerical ~~ 1*Numerical
Spatial ~~ 1*Spatial

### Korelace / kovariance latentních promenných jsou defaultne volné, ale prece:

Verbal ~~ Numerical
Numerical ~~ Spatial
Verbal ~~ Spatial

'
fit3Fc <- cfa(model3Fc, sample.cov = cormatrix, sample.nobs = 696, estimator = "ML")

summary(fit3Fc)

lavInspect(fit3Fc)$lambda
lavInspect(fit3Fc)$psi
lavInspect(fit3Fc)$theta

lavInspect(fit3Fc, what = "est")$lambda
lavInspect(fit3Fc, what = "est")$psi
lavInspect(fit3Fc, what = "est")$theta

residuals(fit3Fc)
corrplot(residuals(fit3Fc)$cov)

fitmeasures(fit3Fc, fit.measures = c("rmsea", 
                                    "tli", "df","npar",
                                    "chisq",
                                    "rmsea.ci.lower","rmsea.ci.upper", "cfi", "srmr"))

#### Srovnávání vnorených modelu - Chi2 test

## Model s nekorelovanými faktory je vnoren v modelu s korelovanými faktory (je jeho speciálním prípadem)

## Modely muzeme srovnat pomocí likelihood-ratio testu

anova(fit3F, fit3Fc)





