#### (c) Hynek Cígler                                                                       ####
#### Katedra psychologie, Fakulta sociálních studií Masarykovy univerzity                   ####
#### Studijní skript ke kurzu PSY028_E: Statistická analýza dat v psychologii (jaro 2018)   ####
#### Seminář 4: odhady různých IRT modelů                                                   ####


# Načte balíčky, data -----------------------------------------------------

library(psych)
library(mirt)
library(readxl)
library(lavaan)

aut <- read_excel("c:/Users/Hynek/Dropbox/MU/Výuka/PSY028 PhD metodologie/Blok3/autonomie.xlsx")
aut <- na.omit(aut[, 3:9])
aut$Aaut2 <- 8-aut$Aaut2
aut$Aaut4 <- 8-aut$Aaut4
aut$Aaut7 <- 8-aut$Aaut7


# Odhad IRT modelů ------------------------------------------------------------

A_RSM <- mirt(aut, 1, "rsm", SE=T) ## Rating scale model
A_PCM <- mirt(aut, 1, "Rasch", SE=T) ## Partial credit model
A_GRM <- mirt(aut, 1, "graded", SE=T) ## Graded response model
A_GPCM <- mirt(aut, 1, "gpcm", SE=T) ## Generalized partial credit model
A_NRM <- mirt(aut, 1, "nominal", SE=T) ## Bockův Nominal response model

A_mixed<- mirt(aut, 1, c("nominal",
                         "rsm",
                         "Rasch",
                         "graded",
                         "gpcm",
                         "gpcm",
                         "gpcm")) ## Bockův Nominal response model
summary(A_mixed)


A_RSM@Fit$SABIC
A_PCM@Fit$SABIC
A_GRM@Fit$SABIC
A_GPCM@Fit$SABIC
A_NRM@Fit$SABIC


M2(A_RSM)
M2(A_PCM)
M2(A_GRM)
M2(A_GPCM)
M2(A_NRM)

anova(A_RSM, A_PCM)
anova(A_GPCM, A_PCM)
anova(A_GPCM, A_NRM)
anova(A_GRM, A_NRM)

windows()
itemplot(A_RSM, 1)
windows()
itemplot(A_PCM, 1)
windows()
itemplot(A_GPCM, 1)
windows()
itemplot(A_GRM, 1)
windows()
itemplot(A_NRM, 1)



# Neparametrické modely ---------------------------------------------------

A_spline <- mirt(aut, 1, "spline", SE=T) ## spline
A_monopoly <- mirt(aut, 1, "monopoly", SE=T) ## spline
summary(A_monopoly)
itemplot(A_monopoly, 2, type="score")
coef(A_monopoly)
anova(A_monopoly, A_GRM)

# Odhad CFA modelů --------------------------------------------------------

A_CFA <- cfa(model = "aut =~ Aaut1+Aaut2+Aaut3+Aaut4+Aaut5+Aaut6+Aaut7", 
             data = aut, estimator="mlr")
A_cCFA <- cfa(model = "aut =~ Aaut1+Aaut2+Aaut3+Aaut4+Aaut5+Aaut6+Aaut7", 
             data = aut, ordered = names(aut))


source("https://raw.githubusercontent.com/tomatitito/grm-plots/master/grm_plots.R")
cumulative_probs(A_cCFA, "Aaut1")
cumulative_response(A_cCFA, "Aaut1")
category_response (A_cCFA, "Aaut1", 3)
information_plot(A_cCFA, "Aaut1")
icc(A_cCFA, "Aaut1")

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

fscores(A_RSM, returnER=T)
fscores(A_PCM, returnER=T)
fscores(A_GRM, returnER=T)
fscores(A_GPCM, returnER=T)
fscores(A_NRM, returnER=T)
omega(aut, nfactors = 1)

semTools::maximalRelia(A_CFA)
semTools::maximalRelia(A_cCFA)

itemfit(A_RSM, empirical.table = 1)
windows()
itemfit(A_RSM, empirical.plot = 1, empirical.CI = .95, group.bins = 10)



# Srovnání ----------------------------------------------------------------

odhady <- data.frame(
  raw = rowSums(aut),
  RSM = fscores(A_RSM),
  PCM = fscores(A_PCM),
  GRM = fscores(A_GRM),
  GPCM = fscores(A_GPCM),
  NRM = fscores(A_NRM),
  CFA = predict(A_CFA),
  cCFA = predict(A_cCFA)
)

names(odhady) <- c("raw", "RSM", "PCM", "GRM", "GPCM", "NRM", "CFA", "cCFA")


windows()
pairs.panels(odhady)

head(odhady)
