############################################
############################################

#Program written by Christopher Moore to apply generalizability theory with R.
#Version: 11/6/12
#Please report errors and suggestions to moor0554 AT umn DOT edu.
#Please visit http://blog.lib.umn.edu/moor0554/canoemoore/ for the most recent version.


############################################
############################################

#Set up R.
options(stringsAsFactors = F, scipen = 10)
windows(record = T)

#Libraries
library(ggplot2); theme_set(theme_bw())
library(plyr)
library(reshape)
library(stringr)
library(lme4)
library(arm)

#A function to examine the generalizability study design.
funk.design.g <- function(g.out, data.g = NULL) {
  if(is.null(data.g) == T) data.g <- attr(g.out, "mer")@frame
  facets <- names(data.g)[names(data.g) != "Score"]
  facets <- as.matrix(outer(facets, facets, paste))
  facets <- facets[lower.tri(facets)]
  funk.design.g.facet <- function(x) {
    data.temp <- x
    data.temp <- str_split_fixed(data.temp, " ", 2)
    data.temp <- data.temp[, order(c(length(levels(data.g[, data.temp[, 1]])), length(levels(data.g[, data.temp[, 2]]))), decreasing = T)]
    formula.design.g <- as.formula(paste("~", paste(data.temp, collapse = " + ")))
    data.temp <- xtabs(data = data.g, formula = formula.design.g)
    name.A <- names(attr(data.temp, "dimnames"))[1]
    design.nested <- sum(colSums(data.temp > 0) > 1) == ncol(data.temp) & sum(rowSums(data.temp == 0) == (ncol(data.temp) - 1)) == nrow(data.temp)
    design.crossed.fully <- sum(colSums(data.temp > 0) == nrow(data.temp)) == ncol(data.temp)
    design.crossed.partially <- sum(colSums(data.temp > 0) == nrow(data.temp)) < ncol(data.temp) & (sum(rowSums(data.temp > 0) == 1) == nrow(data.temp)) == F
    if(design.nested) {
      design.g <- "strictly nested in"
    } else {
      if(design.crossed.fully) {
        design.g <- "fully crossed with"
      } else {
        if(design.crossed.partially) {
          design.g <- "partially crossed with"
        } else {
          design.g <- "no interaction with"
        }
      }
    }
    names(attr(data.temp, "dimnames"))[1] <- paste(names(attr(data.temp, "dimnames"))[1], " (", design.g, " ", names(attr(data.temp, "dimnames"))[2], ")", sep = "")
    design.nested <- sum(rowSums(data.temp > 0) > 1) == nrow(data.temp) & sum(colSums(data.temp == 0) == (nrow(data.temp) - 1)) == ncol(data.temp)
    design.crossed.fully <- sum(rowSums(data.temp > 0) == ncol(data.temp)) == nrow(data.temp)
    design.crossed.partially <- sum(rowSums(data.temp > 0) == ncol(data.temp)) < nrow(data.temp) & (sum(colSums(data.temp > 0) == 1) == ncol(data.temp)) == F
    if(design.nested) {
      design.g <- "strictly nested in"
    } else {
      if(design.crossed.fully) {
        design.g <- "fully crossed with"
      } else {
        if(design.crossed.partially) {
          design.g <- "partially crossed with"
        } else {
          design.g <- "no interaction with"
        }
      }
    }
    names(attr(data.temp, "dimnames"))[2] <- paste(names(attr(data.temp, "dimnames"))[2], " (", design.g, " ", name.A, ")", sep = "")
    data.temp
  }
  alply(facets, 1, funk.design.g.facet)
}

#A function to calculate means of each level of each facet.
funk.means.g <- function(g.out, data.g = NULL) {
  if(is.null(data.g) == T) data.g <- attr(g.out, "mer")@frame
  facets <- names(data.g)[names(data.g) != "Score"]
  funk.means.g.facet <- function(x) {
    data.temp <- x
    data.temp <- ddply(data.g, data.temp, summarize, Score = mean(Score, na.rm = T))
    names(data.temp)[names(data.temp) == x] <- "Level"
    data.frame("Factor" = x, data.temp)
  }
  adply(facets, 1, funk.means.g.facet)[, -1]
}

#A function to conduct a generalizability study (estimate variance components).
funk.g <- function(data.g, formula.g) {
  lmer.out <- lmer(data = data.g, formula = formula.g)
  var.comp <- ldply(VarCorr(lmer.out))
  names(var.comp) <- c("Factor", "Variance")
  var.comp <- rbind(var.comp, data.frame("Factor" = "Residual", "Variance" = attr(VarCorr(lmer.out), "sc") ^ 2))
  var.comp$Percent <- round(var.comp$Variance / sum(var.comp$Variance) * 100, 1)
  attr(var.comp, "mer") <- lmer.out
  class(var.comp) <- c("data.frame", "G")
  var.comp
}

#A function to plot random intercepts as either a dot plot with confidence intervals or a histogram (i.e., vertical ruler plots).
funk.plot.ruler <- function(g.out, type.ruler = "dotplot"){
  lmer.out <- attr(g.out, "mer")
  funk <- function(x) {
    data.frame(
      "ID" = row.names(x), 
      "Estimate" = x$"(Intercept)"
    )
  }
  data.plot <- data.frame(ldply(ranef(lmer.out), funk), "se" = unlist(se.ranef(lmer.out)))
  data.plot$Lower <- with(data.plot, Estimate - 2 * se)
  data.plot$Upper <- with(data.plot, Estimate + 2 * se)
  data.plot$Name <- str_replace_all(paste(paste(str_split_fixed(data.plot$.id, ":", 2)[, 1], str_split_fixed(data.plot$ID, ":", 2)[, 1]), paste(str_split_fixed(data.plot$.id, ":", 2)[, 2], str_split_fixed(data.plot$ID, ":", 2)[, 2]), sep = ", "), ",  ", "")
  data.plot <- subset(data.plot, str_detect(.id, ":") == F)
  data.plot$ID <- factor(data.plot$ID, levels = unique(data.plot$ID), ordered = T)
  if(type.ruler == "histogram") {
    ggplot.out <- ggplot(data = data.plot, aes(x = Estimate)) + 
      geom_histogram(color = "darkgrey", fill = "lightblue") + 
      facet_wrap( ~ .id, nrow = 1) + 
      coord_flip() + 
      labs(title = "Vertical ruler: Random intercepts", y = "")
  } else {
    ggplot.out <- ggplot(data = data.plot, aes(x = ID, y = Estimate)) + 
      geom_errorbar(width = 0, aes(ymax = Upper, ymin = Lower), color = "blue") + 
      geom_point(size = 8, color = "white", shape = 18) + 
      geom_text(size = 4, aes(label = ID)) + 
      facet_grid(. ~ .id, scales = "free_x", space = "free") + 
      scale_x_discrete(breaks = NULL) + 
      labs(title = "Vertical ruler: Random intercepts", x = "") + 
      theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
  }
  ggplot.out
}

#A function to perform decision studies.
funk.d <- function(g.out, N = g.out$N, object = "Person", type = "both") {
  if(is.null(N) | length(N) != nrow(g.out)) {
    stop("N was not specified correctly.")
  } else {
    if(type %in% c("both", "absolute")) {
      cases.keep <- g.out$Factor != object
      AbsoluteErrorVariance <- g.out$Variance[cases.keep] / N[cases.keep]
      AbsoluteErrorVariance <- sum(AbsoluteErrorVariance)
      AbsoluteSEM <- sqrt(AbsoluteErrorVariance)
      Dependability <- g.out$Variance[! cases.keep] / (g.out$Variance[! cases.keep] + AbsoluteErrorVariance)
    } else {
      AbsoluteErrorVariance <- NA
      AbsoluteSEM <- NA
      Dependability <- NA
    }
    if(type %in% c("both", "relative")) {
      cases.keep <- g.out$Factor != object & (str_detect(g.out$Factor, object) | g.out$Factor == "Residual")
      RelativeErrorVariance <- g.out$Variance[cases.keep] / N[cases.keep]
      RelativeErrorVariance <- sum(RelativeErrorVariance)
      RelativeSEM <- sqrt(RelativeErrorVariance)
      Generalizability <- g.out$Variance[g.out$Factor == object] / (g.out$Variance[g.out$Factor == object] + RelativeErrorVariance)
    } else {
      RelativeErrorVariance <- NA
      RelativeSEM <- NA
      Generalizability <- NA
    }
    data.frame(RelativeErrorVariance, RelativeSEM, Generalizability, AbsoluteErrorVariance, AbsoluteSEM, Dependability)
  }
}


############################################
############################################

#Read and process synthetic data set #4 from Brennan's (2001) _Generalizability Theory_ in order to replicate GENOVA output on page 469.
data.g <- read.csv(text = "
Rater,Person,Item1,Item2,Item3,Item4,Item5,Item6,Item7,Item8,Item9,Item10,Item11,Item12
1,1,5,6,5,5,,,,,,,,
1,2,9,3,7,7,,,,,,,,
1,3,3,4,3,3,,,,,,,,
1,4,7,5,5,3,,,,,,,,
1,5,9,2,9,7,,,,,,,,
1,6,3,4,3,5,,,,,,,,
1,7,7,3,7,7,,,,,,,,
1,8,5,8,5,7,,,,,,,,
1,9,9,9,8,8,,,,,,,,
1,10,4,4,4,3,,,,,,,,
2,1,,,,,5,3,4,5,,,,
2,2,,,,,7,5,5,5,,,,
2,3,,,,,5,3,3,5,,,,
2,4,,,,,3,1,4,3,,,,
2,5,,,,,7,7,3,7,,,,
2,6,,,,,3,3,6,3,,,,
2,7,,,,,7,5,5,7,,,,
2,8,,,,,7,5,5,4,,,,
2,9,,,,,6,6,6,5,,,,
2,10,,,,,3,5,6,5,,,,
3,1,,,,,,,,,6,7,3,3
3,2,,,,,,,,,7,7,5,2
3,3,,,,,,,,,6,5,1,6
3,4,,,,,,,,,5,3,3,5
3,5,,,,,,,,,2,7,5,3
3,6,,,,,,,,,4,5,1,2
3,7,,,,,,,,,5,5,5,4
3,8,,,,,,,,,3,2,1,1
3,9,,,,,,,,,5,8,1,1
3,10,,,,,,,,,5,7,1,1
")
data.g
data.g <- melt(data.g, id.vars = c("Rater", "Person"), variable_name = "Item")
data.g <- na.omit(data.g)
head(data.g)
data.g$Rater <- factor(data.g$Rater)
data.g$Person <- factor(data.g$Person)
data.g$Item <- factor(as.integer(str_replace(as.character(data.g$Item), fixed("Item"), "")))
names(data.g)[names(data.g) == "value"] <- "Score"

#Examine design.
funk.design.g(data.g = data.g)

#Examine means.
funk.means.g(data.g = data.g)

#Conduct G study (estimate variance components).
formula.g <- as.formula("Score ~ (1 | Person) + (1 | Person : Rater) + (1 | Rater) + (1 | Rater : Item)")
g.out <- funk.g(data.g, formula.g)
g.out

#Plot vertical rulers.
funk.plot.ruler(g.out)
funk.plot.ruler(g.out, type.ruler = "histogram")

#Conduct D study by "looping" over combinations of possible counts.
array.denominator <- expand.grid(list("Items" = 1 : 8, "Raters" = 1 : 6))
head(array.denominator)
tail(array.denominator)
funk.N <- function(x) {
  Items <- x$Items
  Raters <- x$Raters
  N <- ifelse(str_detect(g.out$Factor, "Rater") | str_detect(g.out$Factor, "Residual"), Raters, 1) * 
    ifelse(str_detect(g.out$Factor, "Item") | str_detect(g.out$Factor, "Residual"), Items, 1)
  data.frame(g.out, N)
}
list.g.out <- alply(array.denominator, 1, funk.N)
list.g.out[1:2]
d.out <- ldply(list.g.out, function(x) funk.d(g.out = x))
head(d.out)
tail(d.out)

#Plot results of D study.
data.plot <- d.out
data.plot$Raters <- factor(data.plot$Raters)
data.plot.observed <- subset(data.plot, Items == 4 & Raters == 3)
funk.plot.d <- function(y, title, ylab) {
  ggplot(data = data.plot, aes_string(x = "Items", y = y)) + 
    geom_line(lwd = 1, aes(linetype = Raters, color = Raters)) + 
    geom_point(data = data.plot.observed, size = 3, aes_string(x = "Items", y = y)) + 
    geom_text(data = data.plot.observed, size = 3, label = "Observed", vjust = -1, aes_string(x = "Items", y = y)) + 
    scale_color_brewer(palette = "YlOrRd") + 
    labs(title = title, y = ylab) + 
    theme(legend.box = "horizontal", legend.position = "bottom")
}
funk.plot.d(y = "RelativeErrorVariance", title = "Decision study: Relative error variance", ylab = expression(paste(hat(sigma) ^ 2, "(", delta, ")")))
funk.plot.d(y = "RelativeSEM", title = "Decision study: Relative standard error of measurement", ylab = expression(paste(hat(sigma), "(", delta, ")")))
funk.plot.d(y = "Generalizability", title = "Decision study: Generalizability of scores for relative decisions", ylab = expression(hat(rho)^2))
funk.plot.d(y = "AbsoluteErrorVariance", title = "Decision study: Absolute error variance", ylab = expression(paste(hat(sigma) ^ 2, "(", Delta, ")")))
funk.plot.d(y = "AbsoluteSEM", title = "Decision study: Absolute standard error of measurement", ylab = expression(paste(hat(sigma), "(", Delta, ")")))
funk.plot.d(y = "Dependability", title = "Decision study: Dependability of scores for absolute decisions", ylab = expression(hat(phi)))


############################################
############################################

#Read and process table 5.1 from Shavelson and Webb's (1991) _Generalizability Theory: A Primer_.
data.g <- read.csv(text = "
Subject,Person,Rater1,Rater2,Rater3
Mathematics,1,4,4,4
Mathematics,2,6,7,6
Mathematics,3,8,7,7
Mathematics,4,6,8,7
Mathematics,5,2,1,1
Mathematics,6,5,4,4
Mathematics,7,4,5,6
Mathematics,8,7,7,6
Reading,1,5,5,6
Reading,2,7,9,5
Reading,3,4,3,2
Reading,4,9,11,7
Reading,5,5,5,3
Reading,6,7,6,5
Reading,7,6,8,9
Reading,8,5,9,9
")
data.g
data.g <- melt(data.g, id.vars = c("Subject", "Person"), variable_name = "Rater")
data.g <- na.omit(data.g)
head(data.g)
data.g$Person <- factor(data.g$Person)
data.g$Rater <- factor(as.integer(str_replace(as.character(data.g$Rater), "Rater", "")))
data.g$Subject <- factor(data.g$Subject)
names(data.g)[names(data.g) == "value"] <- "Score"

#Examine design.
funk.design.g(data.g = data.g)

#Examine means.
funk.means.g(data.g = data.g)

#Conduct G study (estimate variance components).
formula.g <- as.formula("Score ~ (1 | Person) + (1 | Rater) + (1 | Subject) + (1 | Person : Rater) + (1 | Person : Subject) + (1 | Subject : Rater)")
g.out <- funk.g(data.g, formula.g)
g.out

#Plot vertical rulers.
funk.plot.ruler(g.out)
funk.plot.ruler(g.out, type.ruler = "histogram")

#Average over conditions of the fixed facet.
facet.fixed <- "Subject"
n.fixed <- length(levels(data.g[, facet.fixed]))
facets.random <- c(names(data.g)[names(data.g) %in% c(facet.fixed, "Score") == F], "Residual")
g.out.orig <- g.out
g.out <- data.frame(
  "Factor" = c("Person", "Rater", "Residual"), 
  "Variance" = rbind(
    subset(g.out, Factor == "Person", select = "Variance") + subset(g.out, Factor == "Person:Subject", select = "Variance") / n.fixed, 
    subset(g.out, Factor == "Rater", select = "Variance") + subset(g.out, Factor == "Subject:Rater", select = "Variance") / n.fixed, 
    subset(g.out, Factor == "Person:Rater", select = "Variance") + subset(g.out, Factor == "Residual", select = "Variance") / n.fixed
  )
) 
g.out$Percent <- round(g.out$Variance / sum(g.out$Variance) * 100, 1)
attr(g.out, "mer") <- attr(g.out.orig, "mer")
g.out

#Conduct D study by "looping" over combinations of possible counts.
array.denominator <- expand.grid(list("Raters" = 1 : 6))
head(array.denominator)
tail(array.denominator)
funk.N <- function(x) {
  Raters <- x$Raters
  N <- ifelse(str_detect(g.out$Factor, "Rater") | str_detect(g.out$Factor, "Residual"), Raters, 1)
  data.frame(g.out, N)
}
list.g.out <- alply(array.denominator, 1, funk.N)
list.g.out[1:2]
d.out <- ldply(list.g.out, function(x) funk.d(g.out = x))
head(d.out)
tail(d.out)

#Plot results of D study.
data.plot <- d.out
data.plot.observed <- subset(data.plot, Raters == 3)
funk.plot.d <- function(y, title, ylab) {
  ggplot(data = data.plot, aes_string(x = "Raters", y = y)) + 
    geom_line(lwd = 1, color = "blue") + 
    geom_point(data = data.plot.observed, size = 3, aes_string(x = "Raters", y = y)) + 
    geom_text(data = data.plot.observed, size = 3, label = "Observed", vjust = -1, aes_string(x = "Raters", y = y)) + 
    labs(title = title, y = ylab) + 
    theme(legend.box = "horizontal", legend.position = "bottom")
}
funk.plot.d(y = "RelativeErrorVariance", title = "Decision study: Relative error variance", ylab = expression(paste(hat(sigma) ^ 2, "(", delta, ")")))
funk.plot.d(y = "RelativeSEM", title = "Decision study: Relative standard error of measurement", ylab = expression(paste(hat(sigma), "(", delta, ")")))
funk.plot.d(y = "Generalizability", title = "Decision study: Generalizability of scores for relative decisions", ylab = expression(hat(rho)^2))
funk.plot.d(y = "AbsoluteErrorVariance", title = "Decision study: Absolute error variance", ylab = expression(paste(hat(sigma) ^ 2, "(", Delta, ")")))
funk.plot.d(y = "AbsoluteSEM", title = "Decision study: Absolute standard error of measurement", ylab = expression(paste(hat(sigma), "(", Delta, ")")))
funk.plot.d(y = "Dependability", title = "Decision study: Dependability of scores for absolute decisions", ylab = expression(hat(phi)))