###--------------------------------------------------------------------------###
### (c) 2022 | Hynek Cígler, Ph.D.                                           ###
### cigler@fss.muni.cz                                                       ###
### Psychology Department, Faculty of Social Studies                         ###
### Masaryk University                                                       ###
### Psychometrika: Měření v psychologii                                      ###
### P13: Praktická ukázka analýz a psychometrické kauzistiky                 ###
###                                                                          ###
###------ Sun Dec  4 20:57:14 2022 ------------------------------------------###




## Tento skript je ve formatu UTF-8
## ěščřžýáíéúů
## Pokud na predchozim radku nevidite spravne diakritiku:
## V RStudiu: File > Reopen with Encoding > UTF-8





# Knihovny a načtení dat -----------------------------------------------------------------------------

## Nové balíčky nainstalujete např. pomocí příkazu (uvozovky jsou důležit!):  install.packages("psych")

library(psych)    ## běžné psychometrické analýzy
library(lavaan)   ## Faktorová analýza
library(semTools) ## doplněk lavaanu
library(semPlot)  ## strukturní diagramy
library(mirt)     ## IRT analýza
library(foreign)  ## práce s externími datasety (např. SPSS/SAV aj.)
source("https://raw.githubusercontent.com/hynekcigler/Functions/master/extract.R") ## extrakt indexů dobré shody z lavaanu


## Test intelektového potenciálu
## Multiple-choice test, součást Krátkého inteligenčního testu (k dispozici v KDM)
## Cílem je identifikace správného objektu z řady
tip <- read.csv2("https://is.muni.cz/el/fss/podzim2022/PSYn4790/um/irt/tip.csv") 

## Rosenbergova škála self-esteemu
rses <- read.spss("https://is.muni.cz/el/fss/podzim2021/PSYn4790/um/106208423/rses_data.sav", 
                  use.value.labels = FALSE, to.data.frame = TRUE)
rses$id <- gsub(" ", "", rses$id) ## úprava ID


# IRT: TIP ------------------------------------------------------------------------------------

View(tip)
describe(tip)


## CTT a FA ------------------------------------------------------------------------------------

alpha(tip[4:32]) ## vnitřní konzistence
fa.parallel(tip[4:32], sim = F, n.iter = 50, error.bars = T) ## paralelní analýza
(tip_efa2 <- fa(tip[4:32], nfactors = 2, rotate = "geominQ")) ## EFA se dvěma faktory a šikmou rotací geomin
plot(tip_efa2)
corPlot(tip_efa2)

## konfirmační faktorová analýza
(cfa_mod1 <- paste0("TIP =~ ", paste(names(tip[4:32]), collapse = "+")))  ## definice modelu
tip_cfa1_int <- cfa(cfa_mod1, tip, estimator = "mlr") ## odhad modelu
summary(tip_cfa1_int, fit = T, std = T, ci = T) ## prohlédnutí modelu

## alternativní EFA s využitím lepšího estimátoru MLR (psych package výše používá minimální rezidua)
tip_efa2_lavaan <- efaUnrotate(tip[4:32], estimator = "mlr", nf = 2) 
summary(tip_efa2_lavaan)
tip_efa2_lavaan_geomin <- oblqRotate(tip_efa2_lavaan, method = "geomin")

## srovnání modelů
semtext(list(tip_cfa1_int, tip_efa2_lavaan), output = "table", cfi = T, type = "robust") ## fit indexy
lavTestLRT(tip_cfa1_int, tip_efa2_lavaan) ## signifikance rozdílu modelů
summary(compareFit(tip_cfa1_int, tip_efa2_lavaan, nested = T)) ## kombinace obojího pomocí package semTools

semPaths(tip_efa2_lavaan, what = "std") ## strukturní diagram CFA

## strukturní diagram lepší EFA s estimátorem MLR
windows(10,15)
fa.diagram(fa.results = tip_efa2_lavaan_geomin@loading, 
           Phi = tip_efa2_lavaan_geomin@phi, 
           cut = 0, simple = F, digits = 2, sort = F)


## Ordinální faktorová analýza -----------------------------------------------------------------

tip_cfa1_ord <- cfa(cfa_mod1, tip, estimator = "wlsmv", ordered = T) ## odhad ordinálního CFA modelu
summary(tip_cfa1_ord, std = T, fit = T, ci = T)
semtext(list(tip_cfa1_int, tip_efa2_lavaan, tip_cfa1_ord), output = "table", cfi = T) ## srovnání všech tří modelů



## IRT -----------------------------------------------------------------------------------------

irt_1PL <- mirt(tip[4:32], 1, itemtype = "Rasch", SE = T) ## Raschův model
summary(irt_1PL) ## stadnardizované faktorové náboje
M2(irt_1PL)  ## shoda modelu s daty
coef(irt_1PL, simplify = T, IRTpars = T) ## nestandardizované diskriminační parametry

irt_2PL <- mirt(tip[4:32], 1, itemtype = "2PL", SE = T) ## 2PL model
summary(irt_2PL)
M2(irt_2PL)
coef(irt_2PL, simplify = T, IRTpars = T)

irt_3PL <- mirt(tip[4:32], 1, itemtype = "3PL", SE = T) ## 3PL model
summary(irt_3PL)
M2(irt_3PL)
coef(irt_3PL, simplify = T, IRTpars = T)


## Srovnání parametrů ordinální CFA a IRT modelu
plot(lavInspect(tip_cfa1_ord, what= "est")$lambda, 
     coef(irt_2PL, simplify = T, IRTpars = T)$items[,"a"], 
     xlab = "ordinální CFA", ylab = "2PL IRT model", 
     main = "Srovnání faktorových nábojů CFA vs. 2PL IRT")
plot(lavInspect(tip_cfa1_ord, what= "est")$tau, 
     coef(irt_2PL, simplify = T, IRTpars = T)$items[, "b"], 
     xlab = "ordinální CFA", ylab = "2PL IRT model", 
     main = "Srovnání obtížností položek CFA vs. 2PL IRT")

anova(irt_1PL, irt_2PL, irt_3PL) ## srovnání modelů


## Alternativní IRT model: 3PL model, ve kterém ale není parametr pseudouhádnutelnosti odhadován, ale... 
## ... je nastaven na teroetickou hodnotu 1 správná ze 6 možných
irt_3PLfixed <- mirt(tip[4:32], 1, itemtype = "2PL", SE = T, guess = 1/6) 
M2(irt_3PLfixed) ## shoda s modelem se prakticky neliší, SABIC je dokonce nižší (a tedy se model zdá být lepší)...
anova(irt_3PLfixed, irt_3PL) ## ... ačkoli je statisticky významně horší. Efekt je ale zanedbatelný, osobně bych volil tento model.

plot(coef(irt_3PL, simplify = T, IRTpars = T)$items[, "b"], 
     coef(irt_3PLfixed, simplify = T, IRTpars = T)$items[, "b"], 
     xlab = "3PL model", ylab = "3PL s fixovanou spodní asymptotou", 
     main = "Srovnání obtížností položek obou 3PL modelů")
plot(coef(irt_3PL, simplify = T, IRTpars = T)$items[, "a"], 
     coef(irt_3PLfixed, simplify = T, IRTpars = T)$items[, "a"], 
     xlab = "3PL model", ylab = "3PL s fixovanou spodní asymptotou", 
     main = "Srovnání diskriminačních parametrů položek obou 3PL modelů")



## odhad reliability všech modelů
fscores(irt_1PL, returnER = T)
fscores(irt_2PL, returnER = T)
fscores(irt_3PL, returnER = T)
fscores(irt_3PLfixed, returnER = T)


## odhad míry latentního rysu
scores_dat <- rowSums(irt_1PL@Data$data)
scores_1PL <- fscores(irt_1PL)
scores_2PL <- fscores(irt_2PL)
scores_3PL <- fscores(irt_3PL)

## srovnání odhadů latentního rysu
windows(15,15); pairs.panels(data.frame(HS = scores_dat, PL1 = c(scores_1PL), 
                                        PL2 = c(scores_2PL), PL3 = c(scores_3PL)), 
                             method = "pearson", main = "Pearsonovy korelace")
windows(15,15); pairs.panels(data.frame(HS = scores_dat, PL1 = c(scores_1PL), 
                                        PL2 = c(scores_2PL), PL3 = c(scores_3PL)), 
                             method = "spearman", main = "Spearmanovy korelace")


## shoda položek s daty
itemfit(irt_1PL)
itemfit(irt_2PL)
itemfit(irt_3PL)

## Srovnání modelovaných a pozorovaných ICC
itemfit(irt_1PL, empirical.plot = 20)
itemfit(irt_2PL, empirical.plot = 20)
itemfit(irt_3PL, empirical.plot = 20)

itemfit(irt_1PL, empirical.plot = 26)
itemfit(irt_2PL, empirical.plot = 26)
itemfit(irt_3PL, empirical.plot = 26)

plot(irt_3PL, type = "infotrace")

## graf ICC s konfidenčním intervalem
itemplot(irt_1PL, item = 25, CE = T)
itemplot(irt_2PL, item = 25, CE = T)
itemplot(irt_3PL, item = 25, CE = T)

## další grafy
plot(irt_2PL, type = "score")
plot(irt_2PL, type = "rxx")




## IRT invariance podle pohlaví --------------------------------------------------------------------

table(tip$zena)
irt_2PL_config <- multipleGroup(tip[4:32], 1, itemtype = "2PL", SE = T, 
                                group = as.factor(tip$zena))
fscores(irt_2PL_config, returnER = T)
coef(irt_2PL_config, simplify=T, IRTpars = T)

irt_2PL_metric <- multipleGroup(tip[4:32], 1, itemtype = "2PL", SE = T, 
                                group = as.factor(tip$zena), 
                                invariance = c("slopes", "free_var"))
coef(irt_2PL_metric, simplify=T, IRTpars = T)
irt_2PL_scalar <- multipleGroup(tip[4:32], 1, itemtype = "2PL", SE = T, 
                                group = as.factor(tip$zena), 
                                invariance = c("slopes", "free_var", 
                                               "intercepts", "free_mean"))
coef(irt_2PL_metric)

irt_2PL_means <- multipleGroup(tip[4:32], 1, itemtype = "2PL", SE = T, 
                               group = as.factor(tip$zena), 
                               invariance = c("slopes", "free_var", 
                                              "intercepts"))

irt_2PL_equal <- multipleGroup(tip[4:32], 1, itemtype = "2PL", SE = T, 
                               group = as.factor(tip$zena), 
                               invariance = c("slopes", "intercepts"))
M2(irt_2PL_config)
M2(irt_2PL_metric)
M2(irt_2PL_scalar)
M2(irt_2PL_means)

anova(irt_2PL_means, irt_2PL_scalar, irt_2PL_metric, irt_2PL_config)


DIF(irt_2PL_scalar, which.par = c("a1", "d"), scheme = "drop", items2test = 15:20)
itemplot(irt_2PL_config, 19)
plot(irt_2PL_config, type = "trace")


# Rosenbergova škála self-esteemu --------------------------------------------------------------------------

table(rses$dataset) ## součástí je 9 datasetů. Budeme pracovat jen s devátým
rses <- rses[rses$dataset == 9, ] ## budeme pracovat 

describe(rses)
names(rses)

## CFA --------------------------------------------------------------------------------------


### Intervalová CFA -----------------------------------------------------------------------------


rses_mod1 <- "
rses =~ rses01 + rses02 + rses03 + rses04 + rses05 + rses06 + rses07 + rses08 + rses09 + rses10
"
rses_mod2 <- "
POS =~ rses01 + rses03 + rses04 + rses07 + rses10
NEG =~ rses02 + rses05 + rses06 + rses08 + rses09
"

rses_mod_bi <- "
rses =~ rses01 + rses02 + rses03 + rses04 + rses05 + rses06 + rses07 + rses08 + rses09 + rses10
spec =~ rses02 + rses05 + rses06 + rses08 + rses09
"

rses_cfa1_int <- cfa(rses_mod1, rses, estimator = "mlr")
rses_cfa2_int <- cfa(rses_mod2, rses, estimator = "mlr")
rses_cfaBI_int <- cfa(rses_mod_bi, rses, estimator = "mlr", std.lv = T, orthogonal = T)

semtext(list(rses_cfa1_int, rses_cfa2_int, rses_cfa2_bi), 
        output = "table", cfi = T, type = "robust")
summary(rses_cfa2_int, std=T, fit=T)
summary(rses_cfaBI_int, std=T, fit=T)

semPlot::semPaths(rses_cfa1_int, what = "std")
semPlot::semPaths(rses_cfa2_int, what = "std")
semPlot::semPaths(rses_cfaBI_int, what = "std")



### Ordinální CFA -------------------------------------------------------------------------------

rses_cfa1_ord <- cfa(rses_mod1, rses, ordered = T)
rses_cfa2_ord <- cfa(rses_mod2, rses, ordered = T)
rses_cfaBI_ord <- cfa(rses_mod_bi, rses, ordered = T, std.lv = T, orthogonal = T)
summary(rses_cfaBI_ord, std=T, fit=T)

semtext(list(rses_cfa1_int, rses_cfa2_int, rses_cfaBI_int, 
             rses_cfa1_ord, rses_cfa2_ord, rses_cfaBI_ord), 
        output = "table", cfi = T, type = "scaled")
semPaths(rses_cfaBI_ord, what = "std", intercepts = F, 
         thresholdColor = "blue", exoVar = F)




## IRT: 1dimenzionální model ------------------------------------------------------------------------

RSES_irt <- mirt(rses[rses$dataset == 9, 2:11], 1, itemtype = "graded", SE = T)
summary(RSES_irt)
coef(RSES_irt, simplify = T, IRTpars = T)
fscores(RSES_irt, returnER = T)

M2(RSES_irt, na.rm=TRUE)
M2(RSES_irt, na.rm=TRUE, type="C2")


# IRT: 2dimenzionální model -------------------------------------------------------------------



## IRT: S-1 bifaktorový model ----------------------------------------------------------------------

mod_bi <- "
RSES = 1-10
NEG = 2, 5, 6, 8, 9
"

RSES_irt2 <- mirt(rses[rses$dataset == 9, 2:11], mod_bi, itemtype = "graded", SE = T)
summary(RSES_irt2)
M2(RSES_irt2, na.rm=TRUE)
M2(RSES_irt2, na.rm=TRUE, type = "C2")
anova(RSES_irt, RSES_irt2)
fscores(RSES_irt2, returnER = T)


# GT: RSES ----------------------------------------------------------------------------------------

## v prvním kroce je potřeba obrátit reverzní položky
pca(rses[2:11])
rses.inv <- rses
reverse.items <- c("rses02", "rses05", "rses06", "rses08", "rses09", 
                   "rses02_2", "rses05_2", "rses06_2", "rses08_2", "rses09_2")
rses.inv[names(rses.inv) %in% reverse.items] <- 5-rses.inv[names(rses.inv) %in% reverse.items]
describe(rses.inv)
pca(rses.inv[2:11])

library(reshape2)
library(lme4)

names(rses.inv)
dat <- rses.inv[rses.inv$dataset == 5, c(12, 2:11, 15:24)]
dat1 <- data.frame(dat[1:11], occasion = 1)
dat2 <- data.frame(dat[c(1, 12:21)], occasion = 2)
names(dat2) <- names(dat1)
head(dat1); head(dat2)

dat.long <- melt(rbind(dat1, dat2), id.vars = c("id", "occasion"))

rses_G <- lmer(value ~ (1 | id) + (1 | variable) + (1 | occasion) + 
                 (1 | id:occasion) + (1 | id:variable) + (1 | occasion:variable), 
               data = dat.long)
summary(rses_G)
vars <- as.data.frame(VarCorr(rses_G))

