###--------------------------------------------------###
### (c) 2021 | Hynek Cígler, Ph.D.                   ###
### hynek.cigler@mail.muni.cz                        ###
### Psychology Department, Faculty of Social Studies ###
### Masaryk University                               ###
### PSYn4790 | Psychometrika: Měření v psychologii   ###
### Přednáška o IRT                                  ###
### Prague, April 9, 2019                            ###
###------ Tue Nov 09 11:31:17 2021 ------------------###


library(mirt)
library(psych)

setwd("c:/Users/Hynek/Dropbox/MU/Výuka/PSY259_479 Psychometrika/PSY479/2021/PPTX/IRT/")



# Demonstrace IRT modelu --------------------------------------------------


# * Příprava dat ----------------------------------------------------------

tip <- read.csv2("tip.csv")
head(tip)
tip <- tip[4:32]


# * Odhad modelů a prohlédnutí parametrů ----------------------------------

## Raschův model
tip1PL <- mirt(tip, 1, itemtype = "Rasch", SE=T)
summary(tip1PL)
coef(tip1PL, simplify=T, IRTpars = T)

## 2PL model
tip2PL <- mirt(tip, 1, itemtype = "2PL", SE=T)
summary(tip2PL)
coef(tip2PL, simplify=T, IRTpars = T)
M2(tip2PL)

## 3PL model
tip3PL <- mirt(tip, 1, itemtype = "3PL", SE=T)
summary(tip3PL)
(tip3PLcoef <- coef(tip3PL, simplify=T, IRTpars = T))
M2(tip3PL)
describe(tip3PLcoef$items[,"g"]) ## deskriptivy spodní asymptoty
1/6 ## očekávaná hodnota spodní asymptoty (je o trochu více než odhad )

## 3PL s přednastavenou spodní asymptotou na očekávanou hodnotu 1/6
tip3PLfixed <- mirt(tip, 1, itemtype = "2PL", SE=T, guess = 1/6) ## zafixování spodní asymptoty na 1/6
summary(tip3PLfixed)
coef(tip3PLfixed, simplify=T, IRTpars = T)
M2(tip3PLfixed)


# * Reliabilita -----------------------------------------------------------

fscores(tip1PL, returnER = T)
fscores(tip2PL, returnER = T)
fscores(tip3PLfixed, returnER = T)
fscores(tip3PL, returnER = T)


# * Srovnání shody modelu s daty ------------------------------------------------

M2(tip1PL)
M2(tip2PL)
M2(tip3PLfixed)
M2(tip3PL)

anova(tip1PL, tip2PL)
anova(tip2PL, tip3PLfixed)
anova(tip3PLfixed, tip3PL)


# * Srovnání shody položek s daty -----------------------------------------

itemfit(tip1PL)
itemfit(tip2PL)

itemfit(tip1PL, fit_stats = "infit", method = "ML")

itemfit(tip1PL, empirical.plot = 16)
itemfit(tip2PL, empirical.plot = 16)
itemfit(tip3PLfixed, empirical.plot = 16)
itemfit(tip3PL, empirical.plot = 16)

itemfit(tip1PL, empirical.plot = 24)
itemfit(tip2PL, empirical.plot = 24)
itemfit(tip3PLfixed, empirical.plot = 24)
itemfit(tip3PL, empirical.plot = 24)


itemfit(tip1PL, empirical.plot = 20)
itemfit(tip2PL, empirical.plot = 20)
itemfit(tip3PLfixed, empirical.plot = 20)
itemfit(tip3PL, empirical.plot = 20)

itemfit(tip1PL, empirical.plot = 28)
itemfit(tip2PL, empirical.plot = 28)
itemfit(tip3PLfixed, empirical.plot = 28)
itemfit(tip3PL, empirical.plot = 28)


itemplot(tip1PL, 20, CE=T)
itemplot(tip2PL, 20, CE=T)
itemplot(tip3PLfixed, 20, CE=T)
itemplot(tip3PL, 20, CE=T)


# * IRT explorační model --------------------------------------------------

fa.parallel(tip, cor = "poly", sim=F) ## zdá se, že máme dva faktory
efa1 <- mirt(tip, 2, "2PL") ## bacha, bude trvat déle!
summary(efa1, rotate = "geominQ") ## standardizované faktorové náboje s šikmou rotací geomin
coef(efa1, simplify = T) ## IRT parametry (bacha, v intercept-slope parametrizaci)
M2(efa1)
M2(tip2PL)
anova(tip2PL, efa1) ## srovnání modelů



# Srovnání 2PL modelu s CFA analýzou ----------------------------------------------------------

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

mod1 <- paste0("TIP =~ ", paste(names(tip), collapse = "+")) ## definice 1dimenzionálního modelu
cfa1int <- cfa(mod1, tip, estimator = "mlr", meanstructure = T) ## 1dimenzionální intervalová CFA
summary(cfa1int, std = T, fit = T)
cfa1cat <- cfa(mod1, tip, ordered = T) ## 1dimenzionální kategorická CFA
summary(cfa1cat, std = T, fit = T)

## srovnání shody modelů
semtext(list(cfa1int, cfa1cat), cfi=T, output = "table") ## obě CFA
M2(tom2PL) ## 2PL IRT

## srovnání parametrů 2PL a ordCFA
cfa_par_cat <- lavInspect(cfa1cat, "est")
cfa_par_int <- lavInspect(cfa1int, "est")

irt_par <- coef(tip2PL, simplify=T, IRTpars = F)


plot(cfa_par_cat$lambda[,1], irt_par$items[, "a1"], xlab = "kategorická CFA", ylab = "2PL IRT", 
     main = "Srovnání faktorových nábojů")
plot(cfa_par_int$lambda[,1], irt_par$items[, "a1"], xlab = "intervalová CFA", ylab = "2PL IRT", 
     main = "Srovnání faktorových nábojů")

plot(cfa_par_cat$tau[,1], irt_par$items[, "d"], xlab = "kategorická CFA", ylab = "2PL IRT", 
     main = "Srovnání thresholdů")
plot(cfa_par_int$nu[,1], irt_par$items[, "d"], 
     xlab = "intervalová CFA (intercepty)", ylab = "2PL IRT (thresholdy)", 
     main = "Srovnání thresholdů a interceptů")



## Ukázka fungování IRT s daty z devátého subtestu baterie TOMAL-SE.
## Test obsahuje 40 binárních položek. Respondent odpovídá ano/ne na otázku, zda už viděl příslušný obličej. 
## Za správnou odpověď je udělen 1 bod.

# Srovnání TCC Raschova a 2PL modelu --------------------------------------
## Tato část skriptu produkuje slidy srovnávající IRT odhady z 1PL a 2PL modelu použité v přednášce

dat <- expand.table(LSAT7) ## načtení dat 

HS <- rowSums(dat) ## pozorované hrubé skóry

irt1 <- mirt(dat, 1, itemtype = "Rasch") ## odhad 1dimenzionálního Raschova modelu
coef(irt1, simplify=T)
summary(irt1)
est1a <- fscores(irt1, "ML") ## maximum-likelihood (ML) odhad
est1b <- fscores(irt1, "WLE") ## weighted ML odhad (aka mírná regrese k průměru)
plot(irt1, type ="score")


irt2 <- mirt(dat, 1, itemtype = "2PL") ## odhad 1dimenzionálního 2PL modelu
coef(irt2, simplify=T)
summary(irt2)
plot(irt2, type ="score")
est2a <- fscores(irt2, "ML") ## maximum-likelihood (ML) odhad 
est2b <- fscores(irt2, "WLE") ## weighted ML odhad (aka mírná regrese k průměru)

theta <- matrix(seq(-6,6, by=.01), ncol=1) ## hodnoty theta pro účely grafu
TCC1 <- expected.test(irt1, theta) ## TCC Raschova modelu
TCC2 <- expected.test(irt2, theta) ## TCC 2PL modelu

plot(theta, TCC1, type = "l", lwd = 2, las = 1,
     xlab = "theta", ylab = "hrubé skóre", main = "TCC (Raschův model)")
points(est1a, HS, cex=1.2, pch=16, col = "darkgreen")
# odpoznámkování následujících dvou řádků nakreslí jak ML, tak i WLE odhady
# points(est1b, HS, cex=1.2, pch=16, col = "darkblue")
# legend("bottomright", c("ML", "WLE"), pch=16, col=c("darkgreen", "darkblue"), pt.cex=1.2,
#        inset = .02, ncol = 2, title = "Způsob odhadu")


plot(theta, TCC2, type = "l", lwd = 2, las = 1,
     xlab = "theta", ylab = "hrubé skóre", main = "TCC (2PL)")
points(est2a, HS, cex=1.2, pch=16, col = "darkgreen")
# odpoznámkování následujících dvou řádků nakreslí jak ML, tak i WLE odhady
# points(est2b, HS, cex=1.2, pch=16, col = "darkblue")
# legend("bottomright", c("ML", "WLE"), pch=16, col=c("darkgreen", "darkblue"), pt.cex=1.2,
#        inset = .02, ncol = 2, title = "Způsob odhadu")
