###--------------------------------------------------###
### (c) 2020 | 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 3: Teorie zobecnitelnosti              ###
###------ Fri Oct 16 21:42:30 2020 ------------------###


# Knihovny ----------------------------------------------------------------

library(psych)
library(lme4)

## nainstaluje a načte hemp package
if(!require("hemp")) {
  if(!require("devtools")) {
    install.packages("devtools")
    library(devtools)
  }
  devtools::install_github("cddesja/hemp")
  library("hemp")
}


# Data --------------------------------------------------------------------

dat <- read.csv2("PZ2020.csv")
head(dat, 20)
describe(dat)

# G-studie ----------------------------------------------------------------

## odhad modelu
gstudie1 <- lmer("value ~ 1 + (1 | hodnotitel) + (1 | jmeno) + (1 | variable) + 
                (1 | hodnotitel:jmeno) + (1 | hodnotitel:variable) + (1 | jmeno:variable)", dat)
summary(gstudie1)

## extrakce rozptylových komponent
G1_var <- as.data.frame(VarCorr(gstudie1))
G1_var$perc <- G1_var$vcov/sum(G1_var$vcov)

## automatizovaná extrakce hemp balíčkem
G1 <- gstudy(gstudie1)


# Vizuální inspekce - predikce D-studie -----------------------------------

dstudy_plot(G1, unit = "jmeno", facets = list("variable" = c(6, 11, 16), 
                                              "hodnotitel" = c(1:5)), 
            g_coef = T)
dstudy_plot(G1, unit = "jmeno", facets = list("variable" = c(6, 11, 16), 
                                              "hodnotitel" = c(1:5)), 
            g_coef = F)



# D-studie ----------------------------------------------------------------

G1_var

(true_var <- G1_var$vcov[4]) ## rozptyl pravého skóre


# * Relativní D-studie ----------------------------------------------------

## Položky jsou relativní, hodnotitelé rovněž (např. protože jejich "přísnost" se odečte od celkového skóre)
## Usuzujeme na hodnocení právě na těchto položkách "průměrným" hodnotitelem
## 2 hodnotitelé, 11 položek
err_rel <- G1_var$vcov[1]/11 + G1_var$vcov[2]/(2*11) + G1_var$vcov[3]/2 + G1_var$vcov[7]/(2*11)
true_var/(true_var + err_rel) ## koeficient zobecnitelnosti


# * Absolutní D-studie ----------------------------------------------------
## Všechny fasety jsou absolutní.
## Usuzujeme na hodnocení libovolným hodnotitelem v libovolně zvolených kritériích.
## 2 hodnotitelé, 11 položek
err_abs <- G1_var$vcov[5]/2 + G1_var$vcov[6]/11 + G1_var$vcov[1]/11 + G1_var$vcov[2]/(2*11) + G1_var$vcov[3]/2 + G1_var$vcov[7]/(2*11)
true_var/(true_var + err_abs) ## koeficient závislosti


# * Reálná D-studie ----------------------------------------------------
## Reliabilitu součtu právě těchto 11 kritérií (a ne jiných) posouzených náhodně vybranými dvěma hodnotiteli.
## toto by bylo reálně použito.
## Hodnotitelé jsou absolutní, položky relativní fasetou
## 2 hodnotitelé, 11 položek
err_semi <- G1_var$vcov[5]/2 + G1_var$vcov[1]/11 + G1_var$vcov[2]/(2*11) + G1_var$vcov[3]/2 + G1_var$vcov[7]/(2*11)
true_var/(true_var + err_semi) ## koeficient zobecnitelnosti



# * Reliabilita přísnosti odhadu hodnotitele ------------------------------
## Právě v těchto 11 kritériích
## 10 náhodně vybraných (absolutních) osob, 11 definovaných (relativních) položek
rater_var <- G1_var$vcov[5]
rater_err_rel <- G1_var$vcov[4]/10 + G1_var$vcov[1]/(10*11) + G1_var$vcov[2]/11 + G1_var$vcov[3]/11 + G1_var$vcov[7]/(10*11)
rater_var/(rater_var + rater_err_rel) ## koeficient zobecnitelnosti



# * Odhad optimálního počtu hodnotitelů -----------------------------------

y <- NULL
r1 <- c(1:15) ## počet hodnotitelů
i1 <- c(5, 11, 15)
for (r in r1) {
  x <- NULL
  for (i in i1) {
    x <- c(x, G1_var$vcov[5]/r + G1_var$vcov[1]/i + G1_var$vcov[2]/(r*i) + G1_var$vcov[3]/r + G1_var$vcov[7]/(r*i))
  }
  y <- rbind(y, x)
}

y <- as.data.frame(y)
dimnames(y) <- list(paste0("r", c(1:15)), c("i5", "i11", "i15"))
y

y <- true_var/(true_var+y)

plot(r1, y$i5, type="b", xlim=c(1, 15), ylim=c(.5,1), xaxt = "n", col=rainbow(3)[1], lwd=3,
     xlab = "Počet hodnotitelů", ylab="koeficient zobecnitelnosti", 
     main = "Přijímací zkoušky do NMGR psychologie PS2020")
lines(r1, y$i11, type="b", col=rainbow(3)[2], lwd=3)
lines(r1, y$i15, type="b", col=rainbow(3)[3], lwd=3)
abline(h=c(y$i11[2], 1), lty=2, col="gray"); abline(v=2, lty=2, col="gray")
points(r1[2], y$i11[2], pch=4, lwd=3, cex=2)
axis(1, c(1:15))
legend("bottomright", c("5", "11", "15"), col=rainbow(3), lwd=3, pch=1, inset=.02, title = "Počet položek")
