###--------------------------------------------------###
### Rosenber Self-Esteem Scale - teaching examlpe    ###
### PSYn4790 Psychometrika: Měření v psychologii     ###
### Faculty of Social Studies, Masaryk University    ###
### Author 1: Hynek Cígler, Ph.D.                    ###
###------ Wed Oct 14 09:17:19 2020 ------------------###



# Load packages -----------------------------------------------------------

library(psych)
library(foreign)
library(lavaan)
library(semTools)
library(semPlot)
source("https://raw.githubusercontent.com/hynekcigler/Functions/master/extract.R")

library(parallel)


# Read and prepare data ---------------------------------------------------------------

dat <- read.spss("rses_data.sav", use.value.labels = FALSE, to.data.frame = TRUE)
dat <- dat[dat$dataset == 9, 2:11]
dat <- na.omit(dat) ## list-wise delete pro zjednodušení


# Descriptives ------------------------------------------------------------

describe(dat)


# CFA model ---------------------------------------------------------------

mod1 <- "
SE =~ rses01 + rses02 + rses03 + rses04 + rses05 + rses06 + rses07 + rses08 + rses09 + rses10
NEG =~ rses02 + rses05 + rses06 + rses08 + rses09
"

cfa1 <- cfa(mod1, dat, estimator = "mlr", orthogonal = T, std.lv = T) ## intervalová CFA
cfa2 <- cfa(mod1, dat, estimator = "WLSMV", ordered = T, orthogonal = T, std.lv = T) ## ordinální CFA (Graded Response IRT model)

semtext(list(cfa1, cfa2), output="table", cfi=T) ## shoda modelu s daty - docela dobrá
semtext(cfa1, cfi=T) ## copy2clipboard

summary(cfa1, std=T, fit=T) ## parametry modelu, který budeme používat pro odhad reliability
semPaths(cfa1, what="std", bifactor = "SE", layout = "tree2", rotation = 2, exoVar = F, edge.label.cex = 0.8, residuals = F, 
         reorder = F)

# FINAL dat a deskriptivy ---------------------------------------------------------------

dat[c(2, 5, 6, 8, 9)] <- 5-dat[c(2, 5, 6, 8, 9)] ## reverse items
x <- alpha(dat)
write.table(cbind(x$item.stats, x$alpha.drop$raw_alpha), "clipboard", sep="\t", dec=",") ## copy2excel

liche <- rowSums(dat[c(1:5)*2-1]) ## liché položky
sude <- rowSums(dat[c(1:5)*2]) ## sudé položky

# Odhady reliability ------------------------------------------------------

alfa <- alpha(dat)$total[1:2] ## Cronbachovo alfa
half <- splitHalf(dat, covar = T) ## různé odhady split-half reliability
print(half, digits=3)
half$maxAB ## nejlepší rozdělení
half$minAB ## nejlepší rozdělení

## odhad split-half nad lichými a sudými položkami
LICH <- c(SB = 2*cor(liche, sude)/(1+cor(liche, sude)), ## Spearman-Brownův vzorec
          lambda4 = 4*cov(liche,sude)/var(rowSums(dat))) ## Gutmanova lambda4


# * GLB koeficienty -------------------------------------------------------

glb(dat) ## různé "domnělé" glb koeficienty
glb.fa(dat) ## Bentlerův koeficient GLB

# * Expolorační omega -----------------------------------------------------


omega1_tot <- omega(dat, nfactors = 1, covar = T)$omega.tot ## McDonaldova omega pro 1faktorové řešení
omega2_tot <- omega(dat, nfactors = 2, covar = T)$omega.tot ## McDonaldova celková omega pro explorační 2faktorové řešení
omega2_h <- omega(dat, nfactors = 2, covar = T)$omega_h ## McDonaldova hierarchická omega pro explorační 2faktorové řešení


# * Model-based reliability -------------------------------------------------

cfa3 <- cfa(mod1, dat, estimator = "mlr", orthogonal = T, std.lv = T) ## CFA na rekodovaných položkách
cfa4 <- cfa(mod1, dat, estimator = "WLSMV", ordered = T, orthogonal = T, std.lv = T) ## ordinální CFA (Graded Response IRT model)

reliability(cfa3, return.total = T)
reliability(cfa4, return.total = T)


# * FSD, odhady faktorových skórů apod. -----------------------------------

maximalRelia(cfa3) ## celkový skór


phi <- lavInspect(cfa3, what="cor.lv")
lambda <- lavInspect(cfa3, what="std")$lambda
sigma <- lavInspect(cfa3, what="cor.ov")
FD <- diag(phi %*% t(lambda) %*% solve(sigma) %*% lambda %*% phi)^0.5
FD ## factor score determinacy

lavInspect(cfa3)