###----------------------------------------------###
### (c) 2018 | Hynek Cígler                      ###
### Masaryk University, Department of Psychology ###
### PSY028: Data Analysis course                 ###
### Block 3: CFA, special topics                 ###
###------ Thu Dec 06 11:08:12 2018 --------------###




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

setwd("c:/Users/Hynek/Dropbox/MU/Výuka/PSY028 PhD metodologie/PSY028_E - SEM")
library(psych)
library(lavaan)
library(semTools)
library(semPlot)
source("https://raw.githubusercontent.com/hynekcigler/Functions/master/extract.R") ## Funkce pro snadnější extrahování fitu
# R version 3.5.1 (2018-07-02)
# other attached packages:
# [1] semPlot_1.1    semTools_0.5-1 lavaan_0.6-3   psych_1.8.4 

# Téma 1: Vícedimenzionální CFA modely ------------------------------------


# * Data ------------------------------------------------------------------

cormatrix <- read.csv(file = "Blok2 - CFA (Adam)/datamatrix.csv", header = T)
cormatrix <- as.matrix(cormatrix) ## Převedeme data do matice
N <- 696 ## sample size
# Proměnné:
#
# Word Meaning
# Sentence Comprehension
# Odd Words
# Mixed Arithmetic
# Remainders
# Missing Numbers
# Gloves
# Boots
# Hatchets


# * 1dimensional model ----------------------------------------------------

model1F <- "
G =~ WrdMean + SntComp + OddWrds + MxdArit + Remndrs + MissNum + Gloves + Boots + Hatchts  # Vsechny manifestní promenné na jednom faktoru
"

cfa1F <- cfa(model1F, sample.cov = cormatrix, sample.nobs = N, estimator = "ML")
summary(cfa1F, fit=T, std=T, ci=T, rsquare=T)
semtext(cfa1F, type = "standard")
semPaths(cfa1F, what="path", whatLabels = "std", layout = "tree", exoVar=F)



# * 3dimensional orthogonal model -----------------------------------------

model3F <- "
Verbal =~ WrdMean + SntComp + OddWrds
Numerical =~ MxdArit + Remndrs + MissNum
Spatial =~ Gloves + Boots + Hatchts
"

cfa3Forth <- cfa(model3Forth, sample.cov = cormatrix, sample.nobs = N, estimator = "ML", orthogonal = T)
summary(cfa3Forth, fit=T, std=T, ci=T, rsquare=T)
semtext(cfa3Forth, type = "standard")
semPaths(cfa3Forth, what="path", whatLabels = "std", layout = "tree", exoVar=F)


# * 3dimensional oblique model --------------------------------------------

cfa3F <- cfa(model3F, sample.cov = cormatrix, sample.nobs = N, estimator = "ML")
summary(cfa3F, fit=T, std=T, ci=T, rsquare=T)
semtext(cfa3F, type = "standard")
semPaths(cfa3F, what="path", whatLabels = "std", layout = "tree", exoVar=F)


# * model comparison ------------------------------------------------------

semtext(list(cfa1F, cfa3Forth, cfa3F), type = "standard", cfi=T, bic=T)
semtext(list(cfa1F, cfa3Forth, cfa3F), type = "standard", cfi=T, bic=T, output = "table")
lavTestLRT(cfa1F, cfa3F)
## totéž manuálně - srov. s lavTestLRT(cfa1F, cfa3F)
(chi <- fitmeasures(cfa1F, fit.measures = c("chisq"))-fitmeasures(cfa3F, fit.measures = c("chisq"))) ## chi2
(df <- fitmeasures(cfa1F, fit.measures = c("df"))-fitmeasures(cfa3F, fit.measures = c("df")))        ## df
pchisq(chi, df, lower.tail = F) ## p-value

lavTestLRT(cfa3Forth, cfa3F)



# * Hierarchický model ----------------------------------------------------

model3Fhier <- "
Verbal =~ WrdMean + SntComp + OddWrds      ## faktor 1. řádu
Numerical =~ MxdArit + Remndrs + MissNum   ## faktor 1. řádu
Spatial =~ Gloves + Boots + Hatchts        ## faktor 1. řádu
G =~ Verbal + Numerical + Spatial          ## faktor 2. řádu
"
cfa3Fhier <- cfa(model3Fhier, sample.cov = cormatrix, sample.nobs = N, estimator = "ML")
summary(cfa3Fhier, fit=T, std=T, ci=T, rsquare=T)
semtext(cfa3Fhier, type = "standard")
semPaths(cfa3Fhier, what="path", whatLabels = "std", layout = "tree", exoVar=F)



# * Bifaktorový model -----------------------------------------------------

model3Fbifac <- "
Verbal =~ NA*WrdMean + SntComp + OddWrds    ## specifický faktor
Numerical =~ NA*MxdArit + Remndrs + MissNum ## specifický faktor
Spatial =~ NA*Gloves + Boots + Hatchts      ## specifický faktor
G =~ NA*WrdMean + SntComp + OddWrds + MxdArit + Remndrs + MissNum + Gloves + Boots + Hatchts ## general factor

Verbal ~~ 1*Verbal
Numerical ~~ 1*Numerical
Spatial ~~ 1*Spatial
G ~~ 1*G
"
cfa3Fbifac <- cfa(model3Fbifac, sample.cov = cormatrix, sample.nobs = N, estimator = "ML", orthogonal = T)
summary(cfa3Fbifac, fit=T, std=T, ci=T, rsquare=T)
inspect(cfa3Fbifac, what = "est")$theta ## MxdArit má odhadnutý negativní rozptyl (Heywood case)
semtext(cfa3Fbifac, type = "standard")
semPaths(cfa3Fbifac, what="path", whatLabels = "std", layout = "tree2", exoVar=F, bifactor = "G", exoCov=F)


# * Upravený bifaktorový model (jeden specifický faktor odstraněný) -------

model3Fbifac_mod <- "
Verbal =~ NA*WrdMean + SntComp + OddWrds    ## specifický faktor
Spatial =~ NA*Gloves + Boots + Hatchts      ## specifický faktor
G =~ NA*WrdMean + SntComp + OddWrds + MxdArit + Remndrs + MissNum + Gloves + Boots + Hatchts ## general factor

Verbal ~~ 1*Verbal
Spatial ~~ 1*Spatial
G ~~ 1*G
"
cfa3Fbifac_mod <- cfa(model3Fbifac_mod, sample.cov = cormatrix, sample.nobs = N, estimator = "ML", orthogonal = T)
summary(cfa3Fbifac_mod, fit=T, std=T, ci=T, rsquare=T)
semtext(cfa3Fbifac_mod, type = "standard")
semPaths(cfa3Fbifac_mod, what="path", whatLabels = "std", layout = "tree2", exoVar=F, bifactor = "G", exoCov=F, 
         reorder=F, manifests=c("WrdMean", "SntComp", "OddWrds", "MxdArit", "Remndrs", "MissNum", "Gloves", "Boots", "Hatchts"))



# * Model comparison ------------------------------------------------------

semtext(list(cfa3F, cfa3Fhier, cfa3Fbifac_mod, cfa3Fbifac), type="standard", output = "table", 
        bic=T, cfi=T)
lavTestLRT(cfa3Fhier, cfa3Fbifac_mod, cfa3Fbifac)



# Téma 2: Modifikační indexy ----------------------------------------------

# * 1faktorový model ------------------------------------------------------

modindices(cfa1F, sort. = T, maximum.number = 20)



# Categorical model -------------------------------------------------------


help(ability) ## 4dimenzionální IQ test
describe(ability)


# * Intervalové řešení ----------------------------------------------------

## korelované řešení
catmodel4F <- "
reason =~ reason.4 + reason.16 + reason.17 + reason.19
letter =~ letter.7 + letter.33 + letter.34 + letter.58
matrix =~ matrix.45 + matrix.46 + matrix.47 + matrix.55
rotate =~ rotate.3 + rotate.4 + rotate.6 + rotate.8
"

cat4Fint <- cfa(catmodel4F, ability, estimator = "MLR") ## použili jsme MLR estimátor - viz později
summary(cat4Fint, std=T, fit=T, ci=T)
semPaths(cat4Fint, what="path", whatLabels = "std", exoVar=F)

## hierarchické řešení
catmodelG <- "
reason =~ reason.4 + reason.16 + reason.17 + reason.19
letter =~ letter.7 + letter.33 + letter.34 + letter.58
matrix =~ matrix.45 + matrix.46 + matrix.47 + matrix.55
rotate =~ rotate.3 + rotate.4 + rotate.6 + rotate.8
g =~ reason + letter + matrix + rotate
"

catGint <- cfa(catmodelG, ability, estimator = "MLR")
summary(catGint, std=T, fit=T, ci=T)
semPaths(catGint, what="path", whatLabels = "std", exoVar=F)

## srovnání korelovaného a hierarchického řešení
semtext(list(cat4Fint, catGint), output="table", bic=T)
lavTestLRT(cat4Fint, catGint)


# * Kategorické řešení ----------------------------------------------------

cat4F <- cfa(catmodel4F, ability, ordered = colnames(ability), parameterization="delta") 
catG <- cfa(catmodelG, ability, ordered = colnames(ability))

summary(cat4F, std=T, fit=T, ci=T)
summary(catG, std=T, fit=T, ci=T)

lavTestLRT(catG, cat4F)


semtext(list(cat4Fint, cat4F, catGint, catG), output="table")



# * Item Characteristic Function ------------------------------------------

## parametry prvního subtestu z hierarchického modelu
parG <- parameterEstimates(catG)
View(parG)
lam1 <- parG$est[1:4]
th1 <- parG$est[21:24]
res1 <- parG$est[37:40]
sc1 <- parG$est[58:61]
int1 <- parG$est[74:77]
latvar1 <- parG$est[53,57] ## rozptyl 1. specifického a general faktoru
x <- seq(-2, 2, .01)

print(sqrt(latvar1[1] + parG$est[1:4] * latvar1[2])) ## SD faktoru

## pravděpodobnost
p1 <- pnorm(q = int1[1] + th1[1], mean = x, sd = sqrt(res1[1]),
            lower.tail = F)
p2 <- pnorm(q = int1[2] + th1[2], mean = x, sd = sqrt(res1[2]),
            lower.tail = F)
p3 <- pnorm(q = int1[3] + th1[3], mean = x, sd = sqrt(res1[3]),
            lower.tail = F)
p4 <- pnorm(q = int1[4] + th1[4], mean = x, sd = sqrt(res1[4]),
            lower.tail = F)

## plot
cols <- rainbow(4)
plot(x, p1, lwd=2, col=cols[1], main = "ICC pro subtest reason", type = "l", xlab = "nestandardizovaná theta")
lines(x, p2, lwd=2, col=cols[2])
lines(x, p3, lwd=2, col=cols[3])
lines(x, p4, lwd=2, col=cols[4])
legend("bottomright", parG$rhs[1:4], inset = .01,
       lwd = 2, col = cols, title = "Položka:")


# * Theta vs. delta parametrizace -----------------------------------------

cat4Fdelta <- cfa(catmodel4F, ability, ordered = colnames(ability), parameterization="delta") 
cat4Ftheta <- cfa(catmodel4F, ability, ordered = colnames(ability), parameterization="theta") 

summary(cat4Fdelta, std=T, fit=T, ci=T)
summary(cat4Ftheta, std=T, fit=T, ci=T)

semtext(list(cat4Fdelta, cat4Ftheta), output="table")






# Estimátory --------------------------------------------------------------

cat4Fml <- cfa(catmodel4F, ability, estimator="ML", missing = "listwise")
summary(cat4Fml, fit=T, std=T)
fitmeasures(cat4Fml, fit.measures = c("baseline.chisq", "baseline.df"))

cat4Ffiml <- cfa(catmodel4F, ability, estimator="ML", missing = "fiml") 
summary(cat4Ffiml, fit=T, std=T)
fitmeasures(cat4Ffiml, fit.measures = c("baseline.chisq", "baseline.df"))


cat4Fmlm <- cfa(catmodel4F, ability, estimator="MLM") 
summary(cat4Fmlm, fit=T, std=T)
fitmeasures(cat4Fmlm, fit.measures = c("baseline.chisq", "baseline.df"))
fitmeasures(cat4Fmlm)[c("baseline.chisq.scaled", "baseline.df.scaled")]

cat4Fmlr <- cfa(catmodel4F, ability, estimator="MLR", missing = "fiml") 
summary(cat4Fmlr, fit=T, std=T)
fitmeasures(cat4Fmlr, fit.measures = c("baseline.chisq", "baseline.df"))
fitmeasures(cat4Fmlr)[c("baseline.chisq.scaled", "baseline.df.scaled")]

cat4Fwls <- cfa(catmodel4F, ability, estimator="WLS") 
summary(cat4Fwls, fit=T, std=T)
fitmeasures(cat4Fwls, fit.measures = c("baseline.chisq", "baseline.df"))
fitmeasures(cat4Fwls)

cat4Fdwls <- cfa(catmodel4F, ability, estimator="DWLS") 
summary(cat4Fdwls, fit=T, std=T)
fitmeasures(cat4Fdwls, fit.measures = c("baseline.chisq", "baseline.df"))
fitmeasures(cat4Fdwls)

cat4Fwlsmv <- cfa(catmodel4F, ability, estimator="WLSMV") 
summary(cat4Fwlsmv, fit=T, std=T)
fitmeasures(cat4Fwlsmv, fit.measures = c("baseline.chisq", "baseline.df"))
fitmeasures(cat4Fwlsmv)[c("baseline.chisq.scaled", "baseline.df.scaled")]

semtext(list(cat4Fml, cat4Ffiml, cat4Fmlm, cat4Fmlr, cat4Fwls, cat4Fdwls, cat4Fwlsmv), output = "table", cfi=T, type="standard")
semtext(list(cat4Fmlm, cat4Fmlr, cat4Fwlsmv), output = "table", cfi=T, type="scaled")
semtext(list(cat4Fmlm, cat4Fmlr, cat4Fwlsmv), output = "table", cfi=T, type="robust")




# Reliabilita -------------------------------------------------------------

reliability(cat4Fint)
reliability(cat4F)

reliability(catGint)
reliability(catG) ## nefunguje

reliabilityL2(catGint, "g")
reliabilityL2(catG, "g") ## zřejmě není korigované oproti thresholdům





# Invariance --------------------------------------------------------------

help(datCat)
describe(datCat)
head(datCat)


help(HolzingerSwineford1939)
head(HolzingerSwineford1939)
describe(HolzingerSwineford1939)


dat <- HolzingerSwineford1939
names(dat)[7:15] <- colnames(cormatrix)
head(dat)


# * automatic -------------------------------------------------------------

## Pozor, measurementInvariance je "depracted" a asi bude brzy odstraněna
measurementInvariance(model=model3F, data=dat, group="school", estimator="mlr", strict=T)

(invar1 <- measurementInvariance(model=model3F, data=dat, group="school", estimator="mlr", strict=T))

(invar2 <- measurementInvariance(model=model3F, data=dat, group="school", estimator="mlr", strict=T,
                      fit.measures = c("cfi.robust", "tli.robust", "rmsea.robust", "srmr")))
semtext(invar1, type="robust", output="table", cfi=T, bic=T)
semtext(invar1, type="robust", output="table", cfi=T, bic=T, 
        srmr.type = "srmr_bentler_nomean") ## stejné SRMR, jako používá v doporučení Chen
lavTestLRT(invar1$fit.configural,
           invar1$fit.loadings,
           invar1$fit.intercepts,
           invar1$fit.residuals,
           invar1$fit.means)


# * Manuální --------------------------------------------------------------

fit_config <- measEq.syntax(configural.model=model3F, data = dat, estimator = "mlr", 
                            group = "school", ID.fac = "std.lv", return.fit=T,
                            group.equal = "")
fit_metric <- measEq.syntax(configural.model=model3F, data = dat, estimator = "mlr", 
                            group = "school", ID.fac = "std.lv", return.fit=T,
                            group.equal = "loadings")
fit_scalar <- measEq.syntax(configural.model=model3F, data = dat, estimator = "mlr", 
                            group = "school", ID.fac = "std.lv", return.fit=T,
                            group.equal = c("loadings", "intercepts"))
fit_residual <- measEq.syntax(configural.model=model3F, data = dat, estimator = "mlr", 
                            group = "school", ID.fac = "std.lv", return.fit=T,
                            group.equal = c("loadings", "intercepts", "residuals"))
fit_mean <- measEq.syntax(configural.model=model3F, data = dat, estimator = "mlr", 
                              group = "school", ID.fac = "std.lv", return.fit=T,
                              group.equal = c("loadings", "intercepts", "means"))

semtext(list(fit_config, fit_metric, fit_scalar, fit_residual, fit_mean), 
        output="table", bic=T, cfi=T)
lavTestLRT(fit_config, fit_metric, fit_scalar, fit_mean)
lavTestLRT(fit_config, fit_metric, fit_scalar, fit_residual)


# * Problém s různými SRMR ------------------------------------------------

fitmeasures(fit_config)[c("srmr", "srmr_bentler", "srmr_bentler_nomean", "srmr_mplus", "srmr_mplus_nomean")]
fitmeasures(fit_metric)[c("srmr", "srmr_bentler", "srmr_bentler_nomean", "srmr_mplus", "srmr_mplus_nomean")]
fitmeasures(fit_scalar)[c("srmr", "srmr_bentler", "srmr_bentler_nomean", "srmr_mplus", "srmr_mplus_nomean")]

semtext(list(fit_config, fit_metric, fit_scalar, fit_residual, fit_mean), output="table", srmr.type = "srmr")
semtext(list(fit_config, fit_metric, fit_scalar, fit_residual, fit_mean), output="table", srmr.type = "srmr_bentler")
semtext(list(fit_config, fit_metric, fit_scalar, fit_residual, fit_mean), output="table", srmr.type = "srmr_bentler_nomean")
semtext(list(fit_config, fit_metric, fit_scalar, fit_residual, fit_mean), output="table", srmr.type = "srmr_mplus")
semtext(list(fit_config, fit_metric, fit_scalar, fit_residual, fit_mean), output="table", srmr.type = "srmr_mplus_nomean")




# * Kategorická invariance ------------------------------------------------

mod.cat <- ' FU1 =~ u1 + u2 + u3 + u4
             FU2 =~ u5 + u6 + u7 + u8 '


# * * Automatic -----------------------------------------------------------

measurementInvarianceCat(model=mod.cat, data = datCat, ordered = paste0("u", 1:8), group="g")

measurementInvarianceCat(model=mod.cat, data = datCat, ordered = paste0("u", 1:8), group="g", 
                         parameterization="theta", strict=T)



# * * Manuálně ------------------------------------------------------------

cat_config <- measEq.syntax(configural.model=mod.cat, data = datCat, ordered = paste0("u", 1:8), 
                            parameterization="theta", group = "g", ID.fac = "std.lv", return.fit=T,
                            group.equal = "")
cat_load <- measEq.syntax(configural.model=mod.cat, data = datCat, ordered = paste0("u", 1:8), 
                            parameterization="theta", group = "g", ID.fac = "std.lv", return.fit=T,
                            group.equal = c("loadings"))
cat_thresh <- measEq.syntax(configural.model=mod.cat, data = datCat, ordered = paste0("u", 1:8), 
                          parameterization="theta", group = "g", ID.fac = "std.lv", return.fit=T,
                          group.equal = c("thresholds"))
cat_metric <- measEq.syntax(configural.model=mod.cat, data = datCat, ordered = paste0("u", 1:8), 
                          parameterization="theta", group = "g", ID.fac = "std.lv", return.fit=T,
                          group.equal = c("loadings", "thresholds"))
cat_scalar <- measEq.syntax(configural.model=mod.cat, data = datCat, ordered = paste0("u", 1:8), 
                            parameterization="theta", group = "g", ID.fac = "std.lv", return.fit=T,
                            group.equal = c("loadings", "thresholds", "intercepts"))
cat_strict <- measEq.syntax(configural.model=mod.cat, data = datCat, ordered = paste0("u", 1:8), 
                            parameterization="theta", group = "g", ID.fac = "std.lv", return.fit=T,
                            group.equal = c("loadings", "thresholds", "intercepts", "residuals"))
cat_means <- measEq.syntax(configural.model=mod.cat, data = datCat, ordered = paste0("u", 1:8), 
                            parameterization="theta", group = "g", ID.fac = "std.lv", return.fit=T,
                            group.equal = c("loadings", "thresholds", "intercepts", "means"))

semtext(list(cat_config, cat_load, cat_thresh, cat_metric, cat_scalar, cat_strict, cat_means), 
        output = "table", cfi=T)

lavTestLRT(cat_config, cat_metric, cat_scalar, cat_strict)
lavTestLRT(cat_config, cat_metric, cat_scalar, cat_means)
