# Generalizability theory with real data: variance components and D studies for
# cleverness ratings of alternate-uses responses (persons x tasks x raters, fully
# crossed) and for teachers' ratings of student essays (essays x criteria x raters,
# sparse). Needs lme4 (Bates, Maechler, Bolker & Walker, 2015); no login or token.
# Both CSV links are pinned to one version of the IRW data. The essay fit takes
# about half a minute.

## ---- helpers
library(lme4)
options(digits = 7, scipen = 5)  # R's default digits, in case a .Rprofile changes it; no e-notation
# Variance components from a crossed random-effects fit, with their shares of the
# total. `labels` renames lme4's grouping factors (id, item, rater and their
# interactions) to the facets of the design.
components <- function(fit, labels) {
  vc <- as.data.frame(VarCorr(fit))
  v <- setNames(vc$vcov, labels[vc$grp])[labels]
  data.frame(variance = round(v, 3), percent = round(100 * v / sum(v), 1))
}
# G coefficients for a D study that averages over n_i levels of the first facet
# (tasks, or criteria) and n_r raters. v holds the seven components in the order
# p, i, r, pi, pr, ir, residual.
g_coef <- function(v, n_i, n_r) {
  rel <- v[4] / n_i + v[5] / n_r + v[7] / (n_i * n_r)
  abs <- rel + v[2] / n_i + v[3] / n_r + v[6] / (n_i * n_r)
  c(relative = v[1] / (v[1] + rel), absolute = v[1] / (v[1] + abs))
}
labels <- c(id = "person", item = "task", rater = "rater", "id:item" = "person x task",
            "id:rater" = "person x rater", "item:rater" = "task x rater",
            Residual = "residual (p x t x r, e)")
form <- resp ~ 1 + (1 | id) + (1 | item) + (1 | rater) +
  (1 | id:item) + (1 | id:rater) + (1 | item:rater)

## ---- fetch-clever
url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse_2:v22_0.Forthmann-2024-cleverness_ratings/rows?format=csv"
cl <- read.csv(url)[, c("id", "item", "rater", "resp")]
c(respondents = length(unique(cl$id)), tasks = length(unique(cl$item)),
  raters = length(unique(cl$rater)), ratings = nrow(cl),
  missing = length(unique(cl$id)) * 3 * 5 - nrow(cl))
table(cl$resp)   # cleverness, 1 to 5; higher = cleverer

## ---- fit-clever
fit_cl <- lmer(form, data = cl)
vc_cl <- components(fit_cl, labels)
vc_cl

## ---- reml-ml
# REML against plain maximum likelihood, for the two facets with few levels.
fit_ml <- lmer(form, data = cl, REML = FALSE)   # warns of a singular fit: task is 0
facet_var <- function(fit) {
  v <- as.data.frame(VarCorr(fit))
  setNames(v$vcov, v$grp)[c("item", "rater")]
}
round(rbind(REML = facet_var(fit_cl), ML = facet_var(fit_ml)), 4)   # item = task

## ---- dstudy-clever
v_cl <- vc_cl$variance
round(g_coef(v_cl, 3, 5), 2)     # the design as run: 3 tasks, 5 raters
round(g_coef(v_cl, 3, 1), 2)     # 3 tasks, 1 rater
round(g_coef(v_cl, 6, 1), 2)     # 6 tasks, 1 rater
round(g_coef(v_cl, 3, 1e6), 2)   # 3 tasks, as many raters as you like
d_grid <- expand.grid(tasks = 1:10, raters = c(1, 2, 5))
d_grid$relative <- mapply(function(t, r) g_coef(v_cl, t, r)[1], d_grid$tasks, d_grid$raters)

## ---- interrater
# Interrater reliability with the tasks fixed: raters score the same tasks, so the
# person x task effect is part of the universe score. One rater, one task:
round(c(tasks_fixed = (v_cl[1] + v_cl[4]) / (v_cl[1] + v_cl[4] + v_cl[5] + v_cl[7]),
        tasks_random = unname(g_coef(v_cl, 1, 1)[1])), 2)

## ---- plot-clever
cols <- c("1" = "#93c5fd", "2" = "#2780e3", "5" = "#c2410c")
plot(NULL, xlim = c(1, 10), ylim = c(0, 1), xlab = "Tasks averaged over",
     ylab = "Relative G coefficient", las = 1, main = "Cleverness ratings: a D study")
for (r in c(1, 2, 5)) with(d_grid[d_grid$raters == r, ],
                           lines(tasks, relative, type = "b", pch = 19, lwd = 2, col = cols[as.character(r)]))
abline(h = 0.8, col = "#999", lty = 2)
legend("bottomright", paste(c(1, 2, 5), c("rater", "raters", "raters")), col = cols,
       lwd = 2, pch = 19, bty = "n")

## ---- fetch-essays
url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse_4:v7_0.teacherjudgements_lohmann_2026_essayratings/rows?format=csv"
es <- read.csv(url)
# The table also holds an expert benchmark score for every essay and criterion;
# the G study is about the teachers, so set it aside.
es <- es[es$rater != "expert_benchmark", c("id", "item", "rater", "resp")]
c(essays = length(unique(es$id)), criteria = length(unique(es$item)),
  teachers = length(unique(es$rater)), ratings = nrow(es))
table(es$resp)   # 1 (very low quality) to 7 (very high quality)
pairs <- unique(es[, c("id", "rater")])
table(table(pairs$rater))   # essays per teacher
table(table(pairs$id))      # teachers per essay

## ---- fit-essays
# The default optimizer stops just short of lme4's convergence tolerance on this
# fit; bobyqa converges and gives the same estimates to three decimals.
fit_es <- lmer(form, data = es, control = lmerControl(optimizer = "bobyqa"))
vc_es <- components(fit_es, sub("t x r", "c x r", sub("task", "criterion", labels)))
vc_es

## ---- dstudy-essays
v_es <- vc_es$variance
sapply(c(1, 2, 3, 5), function(r) round(g_coef(v_es, 4, r), 2))   # 4 criteria, 1-5 raters
round(g_coef(v_es, 1e6, 1), 2)    # one rater, as many criteria as you like
round(sqrt(v_es[c(1, 3)]), 2)     # SD of essays' universe scores; SD of teachers' severity

## ---- balanced-essays
# A balanced part of the essay data: the essays rated by exactly three teachers,
# each on all four criteria. Treating each essay's teachers as its own (raters
# nested within essays) makes the design (r:p) x c, the classic balanced case.
n_teach <- table(pairs$id)
bal <- es[es$id %in% names(n_teach)[n_teach == 3], ]
c(essays = length(unique(bal$id)), teachers = length(unique(bal$rater)), ratings = nrow(bal))
table(table(bal$id))   # ratings per essay: 3 teachers x 4 criteria = 12 for every essay

## ---- anova-essays
# The classic route, by hand: sums of squares, mean squares, then the
# expected-mean-squares equations solved for the five components (see the callout).
# ave() returns each row's group mean, so the sums run over all ratings.
n_p <- length(unique(bal$id)); n_r <- 3; n_c <- 4
m    <- mean(bal$resp)
m_p  <- ave(bal$resp, bal$id)                 # essay means
m_c  <- ave(bal$resp, bal$item)               # criterion means
m_rp <- ave(bal$resp, bal$id, bal$rater)      # teacher-within-essay means
m_pc <- ave(bal$resp, bal$id, bal$item)       # essay x criterion means
ss <- c(p = sum((m_p - m)^2), c = sum((m_c - m)^2), "r:p" = sum((m_rp - m_p)^2),
        pc = sum((m_pc - m_p - m_c + m)^2), "rc:p,e" = sum((bal$resp - m_rp - m_pc + m_p)^2))
df <- c(n_p - 1, n_c - 1, n_p * (n_r - 1), (n_p - 1) * (n_c - 1), n_p * (n_r - 1) * (n_c - 1))
ms <- ss / df
round(cbind(df, ms), 3)
v_anova <- c(essay = (ms[["p"]] - ms[["r:p"]] - ms[["pc"]] + ms[["rc:p,e"]]) / (n_r * n_c),
             criterion = (ms[["c"]] - ms[["pc"]]) / (n_p * n_r),
             "rater within essay" = (ms[["r:p"]] - ms[["rc:p,e"]]) / n_c,
             "essay x criterion" = (ms[["pc"]] - ms[["rc:p,e"]]) / n_r,
             residual = ms[["rc:p,e"]])

## ---- reml-balanced
# The same design by REML: id:rater is a teacher within an essay.
fit_bal <- lmer(resp ~ 1 + (1 | id) + (1 | item) + (1 | id:rater) + (1 | id:item), data = bal)
vb <- as.data.frame(VarCorr(fit_bal))
v_reml <- setNames(vb$vcov, vb$grp)[c("id", "item", "id:rater", "id:item", "Residual")]
# The full sparse fit, its seven components combined to match the nested design:
# rater within essay = rater + essay x rater; residual = criterion x rater + residual.
vf <- as.data.frame(VarCorr(fit_es))
vf <- setNames(vf$vcov, vf$grp)
v_full <- c(vf[["id"]], vf[["item"]], vf[["rater"]] + vf[["id:rater"]], vf[["id:item"]],
            vf[["item:rater"]] + vf[["Residual"]])
round(cbind(ANOVA = v_anova, REML = v_reml, "full sparse fit" = v_full), 3)
signif(max(abs(v_anova - v_reml)), 2)   # largest gap between ANOVA and REML