###--------------------------------------------------###
### (c) 2019 | Hynek Cígler, Ph.D.                   ###
### hynek.cigler@mail.muni.cz                        ###
### Psychology Department, Faculty of Social Studies ###
### Masaryk University                               ###
### Generalizability Theory workshop                 ###
### Prague, April 9, 2019                            ###
###------ Mon Apr 08 11:44:26 2019 ------------------###


# Load/instal packages ----------------------------------------------------

if(!require(devtools)) install.packages("devtools")
library(devtools)
if(!require(hemp)) install_github("cddesja/hemp")
library(hemp)

library(psych)
library(lme4)
library(reshape2)



# Example 1: Ability ------------------------------------------------------

# * prepare data ----------------------------------------------------------

?ability           ## binary ability items 
head(ability)
describe(ability)

dat1 <- as.data.frame(ability[,1:4]) ## only reasoning items used
alpha(dat1[,1:4])   ## Cronbach alpha

dat1$ID <- c(1:nrow(dat1)) ## add an ID variable
head(dat1)

dat1_long <- melt(dat1, id.vars = "ID") ## convert to long format
head(dat1_long)


# * G-study ---------------------------------------------------------------

## Design: p x i

G1_mod <- lmer(value ~ (1 | ID) + (1 | variable), dat1_long) ## estimate model
summary(G1_mod)  ## summary
G1 <- gstudy(G1_mod) ## extract variance percentage using hemp package
G1
G1var <- as.data.frame(VarCorr(G1_mod)) ## extract variance manually (in case you need it or you don't want to use hemp package)
G1var

# * D-study automatically -------------------------------------------------

## Design: p x I
dstudy(G1, n = c("variable" = 4), unit = "ID") ## Dstudy for 4 variables
ICC(dat1[,1:4])$results ## as default, lmer function (ML estimator) is used. It possible to set lmer=F to use OLS estimator using aov function
dstudy(G1, n = c("variable" = 1), unit = "ID") ## Dstudy for 1 variable

## comparison of alpha and ICC(3,4) - list-wise delete to achieve the same results
alpha(na.omit(dat1[,1:4]))$total
ICC(na.omit(dat1[,1:4]))$results

G1_mod2 <- lmer(value ~ (1 | ID), dat1_long) ## estimate model for ICC1
G1b <- gstudy(G1_mod2) ## extract variance percentage using hemp package as for ICC1
dstudy(G1b, n = c("Residual" = 1), unit = "ID") ## ICC(1,1)
dstudy(G1b, n = c("Residual" = 1), unit = "ID") ## ICC(1,4)

# * Manual D-study --------------------------------------------------------

D1rel <- G1var$vcov[G1var$grp == "Residual"]/4  ## relative D-study error variance
D1abs <- G1var$vcov[G1var$grp == "variable"]/4 + G1var$vcov[G1var$grp == "Residual"]/4 ## absolute D-study error variance

G1var$vcov[G1var$grp == "ID"]/(G1var$vcov[G1var$grp == "ID"] + D1rel) ## generalizability coefficient
G1var$vcov[G1var$grp == "ID"]/(G1var$vcov[G1var$grp == "ID"] + D1abs) ## dependability coefficient

## differences are caused because hemp package is rounding variance components to 4 decimal places
## @&##11!!!@@@ WTF!


# * Visual inspection -----------------------------------------------------

dstudy_plot(G1, unit = "ID", facets = list("variable" = c(1:20)), g_coef = T)
dstudy_plot(G1, unit = "ID", facets = list("variable" = c(1:20)), g_coef = F)


# Example 2: Writing (2facet design) --------------------------------------

# * inspect data ----------------------------------------------------------

?writing
head(writing)
table(writing$students) ## 10 students
table(writing$prompt) ## 5 prompt
table(writing$rater) ## 10 rater


# * G-study ---------------------------------------------------------------

## Design: p x i x r
G2_mod <- lmer(scores ~ (1 | students) + (1 | prompts) + (1 | raters) + 
                 (1 | students:prompts) + (1 | students:raters) + (1 | prompts:raters), 
               writing, REML=T, 
               control = lmerControl(optimizer ="Nelder_Mead")) ## this row helps with some convergence issues
summary(G2_mod)
G2 <- gstudy(G2_mod)
G2

# * D-study automatic -----------------------------------------------------

dstudy(G2, n=c("raters" = 2, "prompts" = 5), unit = "students")
dstudy_plot(G2, unit = "students", 
            facets = list("prompts" = c(1:10), "raters" = c(1:5)), 
            g_coef = F)


# * different manual D-studies --------------------------------------------

G2$gstudy.out

## 2 random raters, 5 random prompts, absolute D-study
D2_e1 <- G2$gstudy.out$Est.Variance[5]/5 + # 5 prompts
  G2$gstudy.out$Est.Variance[6]/2 + ## 2 raters
  G2$gstudy.out$Est.Variance[1]/5 + # 1 student * 5 prompts
  G2$gstudy.out$Est.Variance[3]/(2*5) + # 2 raters * 5 prompts
  G2$gstudy.out$Est.Variance[2]/2 + # 2 raters * 1 student
  G2$gstudy.out$Est.Variance[7]/(1*2*5) # residual
G2$gstudy.out$Est.Variance[4]/(G2$gstudy.out$Est.Variance[4] + D2_e1) ## dependability coefficient

## 2 random raters (absolute), 5 fixed prompts (relative) D-study
D2_e2 <- G2$gstudy.out$Est.Variance[6]/2 + ## 2 raters
  G2$gstudy.out$Est.Variance[1]/5 + # 1 student * 5 prompts
  G2$gstudy.out$Est.Variance[3]/(2*5) + # 2 raters * 5 prompts
  G2$gstudy.out$Est.Variance[2]/2 + # 2 raters * 1 student
  G2$gstudy.out$Est.Variance[7]/(1*2*5) # residual
G2$gstudy.out$Est.Variance[4]/(G2$gstudy.out$Est.Variance[4] + D2_e2) 

## 2 fixed raters (relative), 5 random prompts (absolute) D-study
D2_e3 <- G2$gstudy.out$Est.Variance[5]/5 + # 5 prompts
  G2$gstudy.out$Est.Variance[1]/5 + # 1 student * 5 prompts
  G2$gstudy.out$Est.Variance[3]/(2*5) + # 2 raters * 5 prompts
  G2$gstudy.out$Est.Variance[2]/2 + # 2 raters * 1 student
  G2$gstudy.out$Est.Variance[7]/(1*2*5) # residual
G2$gstudy.out$Est.Variance[4]/(G2$gstudy.out$Est.Variance[4] + D2_e3) ## dependability coefficient

## 2 fixed raters, 5 fixed prompts, relative D-study
D2_e4 <- G2$gstudy.out$Est.Variance[1]/5 + # 1 student * 5 prompts
  G2$gstudy.out$Est.Variance[3]/(2*5) + # 2 raters * 5 prompts
  G2$gstudy.out$Est.Variance[2]/2 + # 2 raters * 1 student
  G2$gstudy.out$Est.Variance[7]/(1*2*5) # residual
G2$gstudy.out$Est.Variance[4]/(G2$gstudy.out$Est.Variance[4] + D2_e4) ## generalizability coefficient

## 1 random rater, 5 fixed prompts, absolute D-study
D2_e5 <- G2$gstudy.out$Est.Variance[6]/1 + ## 1 raters
  G2$gstudy.out$Est.Variance[1]/5 + # 1 student * 5 prompts
  G2$gstudy.out$Est.Variance[3]/(1*5) + # 1 raters * 5 prompts
  G2$gstudy.out$Est.Variance[2]/1 + # 1 raters * 1 student
  G2$gstudy.out$Est.Variance[7]/(1*1*5) # residual
G2$gstudy.out$Est.Variance[4]/(G2$gstudy.out$Est.Variance[4] + D2_e5) ## dependability coefficient

## 1 fixed rater, 5 fixed prompts, relative D-study
D2_e6 <- G2$gstudy.out$Est.Variance[1]/5 + # 1 student * 5 prompts
  G2$gstudy.out$Est.Variance[7]/(1*1*5) # residual
G2$gstudy.out$Est.Variance[4]/(G2$gstudy.out$Est.Variance[4] + D2_e6) ## dependability coefficient


##     The equations above are biased estimators as the sum of variance is not 100%, ... 
## ... but the difference to an unbiased estimator should be usually very small. ...
## ... An unbiased estimator should extend the universe score variance by some ...
## ... interaction terms. 



# Example 3: controling for age -------------------------------------------

iq <- read.csv2("iq.csv")
describe(iq)
head(iq)


# * Cronbach alpha --------------------------------------------------------

# total alpha
alpha(iq[3:22])$total

# separate alpha for cohorts is much lower
x <- NULL
for (i in c(5:10)) {
  cat(paste0("\n", "Age: ", i, "\n"))
  y <- alpha(iq[iq$age == i, 3:22], warnings = F)$total
  print(y)
  x <- c(x, as.numeric(y[1]))
}
x
describe(x) ## mean, median, min/max alphas for cohorts
plot(c(5:10), x, type = "b", xlab="age", ylab="Cronbach's alpha", 
     main = "Cronbach's alpha by age") ## plot alpha by age

# * prepare data ----------------------------------------------------------

iq_long <- melt(iq, id.vars = c("ID", "age")) ## convert to long format
head(iq_long)


# * G-study ---------------------------------------------------------------

## design: ((p x i)|age)

g3_mod <- lmer(value ~ I(scale(age)**1) + I(scale(age)**2) + (1 | variable) + (1 | ID), 
               data = iq_long) ## partial a variance explained by age and squared age of random effect
summary(g3_mod)

g3 <- gstudy(g3_mod)
g3


# * D-study ---------------------------------------------------------------

dstudy(g3, n=c("variable" = 20), unit = "ID") ## compare with the previous plot

# Example 4: Unit of measurement is not respondent ------------------------


# * read and prepare data -------------------------------------------------

dat4 <- read.csv("matematika.csv", sep="\t")
describe(dat4) ## unique is a unique identificator of class*school

dat4_long <- melt(dat4, id.vars = c("ID", "class", "school", "unique")) ## convert to long format
head(dat4_long)
alpha(dat4[3:14])$total

# * G-study ---------------------------------------------------------------

# Model: I x (P:C)
g4_model <- lmer(value ~ (1 | ID) + (1 | variable) + (1 | unique) + 
                   (1 | variable:unique), 
                 data = dat4_long, REML=T, 
                 control = lmerControl(optimizer ="Nelder_Mead"))
summary(g4_model)
g4 <- gstudy(g4_model)
g4


# * D-study ---------------------------------------------------------------

## for person within class only
dstudy(g4, n=c("variable" = 12), unit = "ID")

## for person relative total (across classes)
(x <- g4$gstudy.out)

d4_rel <- x$Est.Variance[5]/(1*12*1) + ## 1 children * 12 residuals * 1 class
  x$Est.Variance[2]/(12*1) ## 12 variables * 1 class
d4_abs <- x$Est.Variance[5]/(1*12*1) + ## 1 children * 12 residuals * 1 class
  x$Est.Variance[2]/(12*1) + ## 12 variables * 1 class
  x$Est.Variance[4]/(1*12*1) ## 1 children * 12 variables * 1 class

(x$Est.Variance[1] + x$Est.Variance[3])/
    (x$Est.Variance[1] + x$Est.Variance[3] + d4_rel) ## generalizability coefficient
(x$Est.Variance[1] + x$Est.Variance[3])/
  (x$Est.Variance[1] + x$Est.Variance[3] + d4_abs) ## dependability coefficient


# ## for class ------------------------------------------------------------

(y <- mean(table(dat4_long$unique)/12)) ## average 21.6 children per class

d4_class_rel <- x$Est.Variance[1]/(1*y) + ## 1 class * 21.6 children
  x$Est.Variance[2]/(1*4) + ## 1 class * 4 variable
  x$Est.Variance[5]/(1*y*12) ## 1 class * 21.6 children * 12 variables
d4_class_abs <- x$Est.Variance[1]/(1*y) + ## 1 class * 21.6 children
  x$Est.Variance[4]/(1*4) + ## 1 class * 4 variable
  x$Est.Variance[2]/(1*4) + ## 1 class * 4 variable
  x$Est.Variance[5]/(1*y*12) ## 1 class * 21.6 children * 12 variables

x$Est.Variance[3]/(x$Est.Variance[3] + d4_class_rel) ## generalizability coefficient
x$Est.Variance[3]/(x$Est.Variance[3] + d4_class_abs) ## generalizability coefficient



# Example 5: Partialy nested design --------------------------------------------------

?writing2
g5_model1 <- lmer(scores ~ (1 | students/raters) + (1 | prompts) + (1 | students:prompts), 
                  data = writing2) ## explicitely specified nestedness
g5_model2 <- lmer(scores ~ (1 | students) + (1 | raters:students) + (1 | prompts) + (1 | students:prompts), 
                  data = writing2) ## specify nestedness manually
g5_model3 <- lmer(scores ~ (1 | students) + (1 | raters) + (1 | prompts) + (1 | students:prompts), 
                  data = writing2) ## ignore nestedness and specify the basic model without interaction
summary(g5_model1)
summary(g5_model2)
summary(g5_model3)
g5_1 <- gstudy(g5_model1)
g5_2 <- gstudy(g5_model2)
g5_3 <- gstudy(g5_model3)
g5_1
g5_2
g5_3
## NOTE: untill the same codes are used in the lowest level, there is no need to explicitely specify nestedness
## NOTE: just don't include meaningless variance components like students:raters (no information about it)

dstudy(g5_1, n=c("prompts"=5, "raters"=1), unit = "students")



# Example 6: Confirmatory multidimensional GT -----------------------------

dat6 <- as.data.frame(ability[,1:8]) ## only reasoning and letter items used

dat6$ID <- c(1:nrow(dat6)) ## add an ID variable
head(dat6)

dat6_long <- melt(dat6, id.vars = "ID") ## convert to long format
head(dat6_long, n=20)
dat6_long$trait <- 0

dat6_long$trait[grep("letter", dat6_long$variable)] <- 1



# * g-study ---------------------------------------------------------------

## In this model, two correlated specific traits are defined.
G6_mod_a <- lmer(value ~ 1 + trait + (0 + as.factor(trait) | ID) + (1 | variable), 
                 dat6_long, REML=T, 
                 control = lmerControl(optimizer ="Nelder_Mead")) ## estimate model
summary(G6_mod_a)

dat6_longB <- dat6_long
dat6_longB$trait[dat6_longB$trait == 0] <- -.5
dat6_longB$trait[dat6_longB$trait == 1] <-  .5

## In this model, the general trait is defined alongside with orthogonal specific ability
G6_mod_b <- lmer(value ~ 1 + trait + (1 | ID) + (0 + trait | ID) + (1 | variable), 
                 dat6_longB, REML=T, 
                 control = lmerControl(optimizer ="Nelder_Mead")) ## estimate model
summary(G6_mod_b)
