###--------------------------------------------------###
### (c) 2021 | Hynek Cígler, Ph.D.                   ###
### hynek.cigler@mail.muni.cz                        ###
### https://is.muni.cz/osoba/hynek.cigler            ###
### Psychology Department, Faculty of Social Studies ###
### Masaryk University                               ###
### Reliability models in Classical Test Theory      ###
###        and latent trait theories                 ###
### Prague, March 23, 2021                           ###
###------ Tue Mar 23 11:46:40 2021 ------------------###


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

spb <- function(r,m=2) {
  m*r/(1+(m-1)*r)
}

dat <- psychTools::holzinger.swineford
dat <- dat[c(8:9,11, 13:14, 16:17, 19:20)]
head(dat)
describe(dat)

fa.parallel(dat, sim=F) ## parallel analysis
corPlot(dat, las=2) ## correlation diagram
fa(dat, nfactors = 3) ## factor analysis with suggested number of factors

## partial correlation network with model selection
parcor <- bootnet::estimateNetwork(dat, "ggmModSelect", corMethod = "spearman")
plot(parcor)



# Model free reliability --------------------------------------------------

## Cronbach's alpha
(alpha.est <- alpha(dat))
alpha.est$total

## Split-half (Revelle's beta plus glb lambda 4)
(split_est <- splitHalf(dat, brute = T, covar = T, check.keys = F)) ## Warning: The default is using correlation, not covariance matrix!
split_est$minrb ## Revelle beta
split_est$maxrb ## maximalized lambda 4

## Bentler's glb
(glb_est <- glb(dat))
glb_est$glb.IC; glb_est$glb.Km; glb_est$glb.Fa ## three alhoritms for Bentler's glb estimation
glb_est$glb.max ## The best of these algoritms

## Spearman-Brown
(corAB <- cor(rowSums(dat[as.logical(glb_est$keys[,"K1"])]), 
             rowSums(dat[as.logical(glb_est$keys[,"K2"])])))
spb(corAB, m = 2)

(corCD <- cor(rowSums(dat[c(1, 3, 5, 7, 9)]), 
              rowSums(dat[c(2, 4, 6, 8)])))
spb(corCD, m = 2)


# Sometging between -------------------------------------------------------

## WARNING: omega function uses as default correlation, not covariance matrix
(omega1 <- omega(dat, nfactors = 1, flip = F, digits = 3, covar = T)) ## unidimensional omega - higher than alpha (tau-equivalency relaxed)
(omega3 <- omega(dat, nfactors = 3, flip = F, digits = 3, covar = T, sl = T)) ## three-dimensional omega - higher than alpha (tau-equivalency relaxed + handling with multidimensionality)
omega1$omega.tot
omega3$omega_h ## hierarchical omega (reliability of general factor)
omega3$omega.tot ## total omega (total explained variance)



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

# * unidimensional model --------------------------------------------------

mod0 <- paste0("G =~ ", paste(names(dat), collapse = " + "))
cat(mod0)
cfa0 <- cfa(mod0, dat, estimator = "uls", orthogonal = F, std.lv = T) ## unweighted least square estimator is chosen for the best comparability to psych::omega function
semtext(cfa0) ## did not fit data
cfa0_rel <- reliability(cfa0) ## McDonald's omega (omega3) is therefore biased
cfa0_rel

# * Three-dimensional model -----------------------------------------------

mod1 <- paste0("spatial =~ ", paste(names(dat)[1:3], collapse = " + "), 
               "\nverbal =~ ", paste(names(dat)[4:6], collapse = " + "), 
               "\ncount =~ ", paste(names(dat)[7:9], collapse = " + "))
cat(mod1)

cfa1 <- cfa(mod1, dat, estimator = "mlr", orthogonal = F, std.lv = T)
summary(cfa1, std=T, fit=F)
semtext(cfa1, type="robust")
(cfa1_rel <- reliability(cfa1, return.total = TRUE))
cfa1_rel[4, ] ## McDonald's omegas hierarchical for each subtest plus omega total


# * Hierarchical model -----------------------------------------------

mod2 <- paste(mod1, "G =~ spatial + verbal + count", sep = "\n")
cat(mod2)

cfa2 <- cfa(mod2, dat, estimator = "mlr", orthogonal = T, std.lv = T)
summary(cfa2, std=T, fit=F)
semtext(list(cfa1, cfa2), type="robust") ## models are same
(cfa2_rel <- reliability(cfa2, return.total = TRUE)) ## reliabilities are same

(cfa2_relH <- reliabilityL2(cfa2, "G")) ## second order reliability: 
##          L1 - items' explained variances, 
##          L2 - 1st order factor explained variances, 
##          partial L1 - items' explained variances controlling 1st order factors



# * S-1 Bifactor model --------------------------------------------------------
## The full bifactor model does not make a sense and provides negative residual variances

mod3 <- paste0(mod1, "\n", "G =~ ", paste(names(dat), collapse = " + "))
cat(mod3)

mod3 <- "
verbal =~ t06_paracomp + t07_sentcomp + t09_wordmean
count =~ t10_addition + t12_countdot + t13_sccaps
G =~ t01_visperc + t02_cubes + t04_lozenges + t06_paracomp + t07_sentcomp + t09_wordmean + t10_addition + t12_countdot + t13_sccaps
"

cfa3 <- cfa(mod3, dat, estimator = "mlr", orthogonal = T, std.lv = T)
summary(cfa3, std=T, fit=F)
semtext(list(cfa1, cfa2, cfa3), type="robust") ## model fits much better

(cfa3_rel <- reliability(cfa3, return.total = T))



# Comparison --------------------------------------------------------------

result <- c(
  alpha = alpha.est$total$raw_alpha, ## Cronbach alpha
  std.alpha = alpha.est$total$std.alpha,  ## std. Cronbach alpha
  beta = split_est$minrb,  ## Revelle beta
  lambda4.max = split_est$maxrb, ## maximalized lamgda4
  glb = glb_est$glb.max,  ## Bentler's glb
  omega1.psych = omega1$omega.tot, ## uni-dimensional omega from psych package
  omega3.psych_h = omega3$omega_h, ## hierarchical omega from psych package (3 factors)
  omega3.psych_t = omega3$omega.tot,  ## total omega from psych package (3 factors)
  omega1.sem = cfa0_rel[3,], ## omega of unidimensional model from lavaan
  omega3.sem_tot = cfa1_rel["omega3", "total"], ## omega total from 3-dimensional model
  omega3.sem_h = cfa2_relH[1],  ## omega hierarchical of general factor
  omega_bi.sem_h = cfa3_rel["omega3", "G"], ## hierarchical omega of general factor in S-1 model
  omega_bi.sem_tot = cfa3_rel["omega3", "total"] ## total omega in S-1 model
)


# * print results ---------------------------------------------------------

sort(result, decreasing  = T) ## all the results

## dimension free results
dimension.free <- c(T, T, F, T, T, F,  F, F, T, T, F, F, T); names(dimension.free) <- names(result)
sort(result[dimension.free], decreasing  = T)

## model based total variance explained
model.based.tot <- c(F, F, F, F, F,  T, F, T, T, T, F, F, T); names(model.based.tot) <- names(result)
sort(result[model.based.tot], decreasing  = T)

## model based general factor variance explained
model.based.h <- c(F, F, T, F, F, F,  T, F, F, F, T, T, F); names(model.based.h) <- names(result)
sort(result[model.based.h], decreasing  = T) ## psych package underestimates omega hierarchical as it uses two-step estimator


