---
title: "Reliability Estimation"
author: "Tate Paxton"
date: "2025-10-3"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, echo = FALSE, message = FALSE, warning=FALSE}
library(readxl)
library(irr)
library(tidyverse)
library(reshape2)
library(knitr)
```

Set Working Directory, change to the folder where you have stored the three relevant excel documents.
```{r}
# Set working directory
setwd("~/Documents/PhD Files/Externships/CC Assess")

# Optional: print to confirm
getwd()

# Note: library designation chunk and data import chunk suppressed in knitted file

```

```{r, echo = FALSE, message = FALSE, warning=FALSE}
# Read in the data, this chunk suppressed in knitted file
irr_all_xlsx <- "IRR All.xlsx"
gort_xlsx    <- "Assignment 2 Scoring.xlsx"
wj_xlsx      <- "Assignment 3 Scoring WJ-IV.xlsx"

# Function to simplify the data read in from many different sheets
read_mat <- function(path, sheet) {
  as.matrix(readxl::read_excel(path, sheet = sheet))
}

sim2.17  <- read_mat(irr_all_xlsx, "Sim 2.17")
voc2.17  <- read_mat(irr_all_xlsx, "Vocab 2.17")
sim7.25  <- read_mat(irr_all_xlsx, "Sim 7.25")
voc7.25  <- read_mat(irr_all_xlsx, "Vocab 7.25")

# Transposed (subjects × raters)
sim2.17_t <- t(sim2.17)
voc2.17_t <- t(voc2.17)
sim7.25_t <- t(sim7.25)
voc7.25_t <- t(voc7.25)

# Woodcock Johnson read in
wjpc     <- read_mat(wj_xlsx, "PCitem")
wjpctot  <- read_mat(wj_xlsx, "PCtot")
wjwa     <- read_mat(wj_xlsx, "WAitem")
wjwatot  <- read_mat(wj_xlsx, "WAtot")
wjor     <- read_mat(wj_xlsx, "ORitem")
wjortot  <- read_mat(wj_xlsx, "ORtot")
wjcombtot<- read_mat(wj_xlsx, "Rawcomb")

# Transposed
wjpc_t     <- t(wjpc)
wjpctot_t  <- t(wjpctot)
wjwa_t     <- t(wjwa)
wjwatot_t  <- t(wjwatot)
wjor_t     <- t(wjor)
wjortot_t  <- t(wjortot)
wjcombtot_t<- t(wjcombtot)

# GORT read in (more involved due to domain breakout)

gortss  <- read_mat(gort_xlsx, "GORT SS")
gorttot <- read_mat(gort_xlsx, "GORT tot")

# Transposed
gortss_t  <- t(gortss)
gorttot_t <- t(gorttot)


# Read with minimal name repair so duplicate headers ("Rate", "Accuracy", etc.) are preserved
gort_raw <- read_excel(gort_xlsx, sheet = "GORT raw", .name_repair = "minimal")

# If the first row in the sheet is actually headers, keep as-is.
# If the first row is data (no header row), set col_names = FALSE instead:
# gort_raw <- read_excel(gort_xlsx, sheet = "GORT raw", col_names = FALSE)

# Ensure we only keep complete 4-column domain blocks
nc <- ncol(gort_raw)
keep_nc <- nc - (nc %% 4)
if (keep_nc == 0) stop("No complete 4-column (Rate/Accuracy/Fluency/Comp.) blocks found.")
if (keep_nc < nc) gort_raw <- gort_raw[, 1:keep_nc, drop = FALSE]

# Split by fixed positions: 1=Rate, 2=Accuracy, 3=Fluency, 4=Comp.
gort_rate <- as.matrix(gort_raw[, seq(1, ncol(gort_raw), by = 4), drop = FALSE])
gort_acc  <- as.matrix(gort_raw[, seq(2, ncol(gort_raw), by = 4), drop = FALSE])
gort_flu  <- as.matrix(gort_raw[, seq(3, ncol(gort_raw), by = 4), drop = FALSE])
gort_comp <- as.matrix(gort_raw[, seq(4, ncol(gort_raw), by = 4), drop = FALSE])

# Optional: transpose to subjects × raters if your ICC expects that orientation
gort_rate_t <- t(gort_rate)
gort_acc_t  <- t(gort_acc)
gort_flu_t  <- t(gort_flu)
gort_comp_t <- t(gort_comp)

```

# Reliability Estimates

The primary analysis estimated interrater reliability following, where possible, the methods outlined by the test publishers. The WISC-V Technical and Interpretive Manual specifies the use of intraclass correlations (ICC), while the GORT-5 reports product-moment correlations between two raters and multiple administrations. The Woodcock Johnson-IV Tests of Achievement authors provide no mention or estimate of interrater reliability. We adopted ICCs as our primary index because they are the most widely recommended for continuous or ordinal scoring of standardized test items and can be applied consistently across measures. 

We sought additional information as needed within their cited sources (McGraw & Wong, 1996; Shrout & Fleiss, 1979) and more recent papers (Koo and Li, 2016). We used a “Case 3” design for total score reliability [ICC(3k)] and item-level reliability [ICC(3,1)], as “each target is rated by each of the same k judges, who are the only judges of interest” (Shrout & Fleiss, 1979, p. 421). This employs a two-way, random effects ANOVA to estimate intraclass correlations.

We also used Krippendorff’s alpha to confirm item-level ICCs (Hayes and Krippendorff, 2007; Marzi et al., 2024) because it is robust to non-interval scaling and rater imbalance. These values were found to be nearly identical to ICC in each case, though we report both to maintain thoroughness. 

\newpage

## WISC V Reliability Estimates

Total Score ICC across two administrations with 11-13 raters
```{r}
# Read in the data
wisc_data <- read_excel(
  "/Users/tatepaxton/Documents/PhD Files/Externships/CC Assess/IRR All.xlsx",
  sheet = "WISC Raw"
)

# Convert to matrix
wisc_data <- as.matrix(wisc_data)

# Ensure it's numeric
wisc_data <- apply(wisc_data, c(1, 2), as.numeric)

# Transpose: now rows = subjects (subtests), columns = raters
wisc_data_t <- t(wisc_data)

# Run ICC for Similarities (only Sim2_17 and Sim7_25)
sim_icc <- psych::ICC(wisc_data_t[c("Sim2_17", "Sim7_25"), ])

# Run ICC for Vocabulary (Voc2_17 and Voc7_25)
voc_icc <- psych::ICC(wisc_data_t[c("Voc2_17", "Voc7_25"), ])

# Print results
sim_icc
voc_icc
```

Item level ICC
```{r}
## --- WISC item-level ICC(3,1) ---
# Two-way mixed, absolute agreement, SINGLE rater (ICC(3,1))

icc_sim_217 <- irr::icc(sim2.17_t, model = "twoway", type = "agreement", unit = "single")  # ICC(3,1)
icc_sim_725 <- irr::icc(sim7.25_t, model = "twoway", type = "agreement", unit = "single")  # ICC(3,1)
icc_voc_217 <- irr::icc(voc2.17_t, model = "twoway", type = "agreement", unit = "single")  # ICC(3,1)
icc_voc_725 <- irr::icc(voc7.25_t, model = "twoway", type = "agreement", unit = "single")  # ICC(3,1)

cat(
  "\nItem-level ICC(3,1):",
  "\n  Similarities (2/17):", round(icc_sim_217$value, 2),
  "\n  Similarities (7/25):", round(icc_sim_725$value, 2),
  "\n  Vocabulary   (2/17):", round(icc_voc_217$value, 2),
  "\n  Vocabulary   (7/25):", round(icc_voc_725$value, 2), "\n"
)

```

Item-level Krippendorff's alpha
```{r}
# WISC item-level Krippendorff's alpha 

# Krippendorff's alpha (interval) with correct orientation: raters x items ---
ka_sim_217 <- irr::kripp.alpha(t(sim2.17_t), method = "interval")
ka_sim_725 <- irr::kripp.alpha(t(sim7.25_t), method = "interval")
ka_voc_217 <- irr::kripp.alpha(t(voc2.17_t), method = "interval")
ka_voc_725 <- irr::kripp.alpha(t(voc7.25_t), method = "interval")

cat(
  "\nKrippendorff's alpha (interval):",
  "\n  Similarities (2/17):", round(ka_sim_217$value, 2),
  "\n  Similarities (7/25):", round(ka_sim_725$value, 2),
  "\n  Vocabulary   (2/17):", round(ka_voc_217$value, 2),
  "\n  Vocabulary   (7/25):", round(ka_voc_725$value, 2), "\n"
)

```

\newpage

## WJ-IV Reliability Estimates

The WJ-IV does not have multiple adiminstrations rated. Nor were we able to convert the raw total scores into standard scores to estimate the total score reliability on a shared scale. Therefore, we will focus on item-level reliability using ICC and K's alpha.

ICC
```{r}
## --- WJ-IV item-level ICC(3,1) ---
# Two-way mixed, absolute agreement, SINGLE rater (ICC(3,1))
# Inputs are items x raters matrices you already made: wjpc_t, wjwa_t, wjor_t

icc_wj_pc <- irr::icc(wjpc_t, model = "twoway", type = "agreement", unit = "single")  # Passage Comprehension
icc_wj_wa <- irr::icc(wjwa_t, model = "twoway", type = "agreement", unit = "single")  # Word Attack
icc_wj_or <- irr::icc(wjor_t, model = "twoway", type = "agreement", unit = "single")  # Oral Reading

cat(
  "\nWJ-IV Item-level ICC(3,1):",
  "\n  Passage Comprehension:", round(icc_wj_pc$value, 2),
  "\n  Word Attack:          ", round(icc_wj_wa$value, 2),
  "\n  Oral Reading:         ", round(icc_wj_or$value, 2), "\n"
)

```

K's Alpha
```{r}
# WJ-IV item-level Krippendorff's alpha (interval) 

ka_wj_pc <- irr::kripp.alpha(t(wjpc_t), method = "interval")
ka_wj_wa <- irr::kripp.alpha(t(wjwa_t), method = "interval")
ka_wj_or <- irr::kripp.alpha(t(wjor_t), method = "interval")

cat(
  "\nWJ-IV Krippendorff's alpha (interval):",
  "\n  Passage Comprehension:", round(ka_wj_pc$value, 2),
  "\n  Word Attack:          ", round(ka_wj_wa$value, 2),
  "\n  Oral Reading:         ", round(ka_wj_or$value, 2), "\n"
)

```

\newpage

## GORT-5 Reliability Estimates

For the total score reliability, we used the scaled scores across the four domains, which gave a shared metric on which we can fairly estimate reliability across different domains. Of note, this means that the total score reliability is a general reliability across all four domains. We calculated item-level reliability for each domain using the raw scores of each story. We did this for both ICC(3,1) and K's alpha estimates. 

ICC
Total Score reliability
```{r}
## --- GORT total-score ICC(3,k) using scaled scores (GORT SS tab) ---
# Rows = subjects, cols = raters
gort_icc_tbl <- psych::ICC(gortss_t)

gort_icc_tbl

```

Item-level
```{r}
## --- GORT item-level ICC(3,1) by domain ---
icc_gort_rate <- irr::icc(gort_rate_t, model = "twoway", type = "agreement", unit = "single")
icc_gort_acc  <- irr::icc(gort_acc_t,  model = "twoway", type = "agreement", unit = "single")
icc_gort_flu  <- irr::icc(gort_flu_t,  model = "twoway", type = "agreement", unit = "single")
icc_gort_comp <- irr::icc(gort_comp_t, model = "twoway", type = "agreement", unit = "single")

cat(
  "\nGORT Item-level ICC(3,1):",
  "\n  Rate:        ", round(icc_gort_rate$value, 2),
  "\n  Accuracy:    ", round(icc_gort_acc$value, 2),
  "\n  Fluency:     ", round(icc_gort_flu$value, 2),
  "\n  Comprehension:", round(icc_gort_comp$value, 2), "\n"
)

```

K's alpha
```{r}
## --- GORT item-level Krippendorff's alpha (interval) by domain ---
ka_gort_rate <- irr::kripp.alpha(t(gort_rate_t), method = "interval")
ka_gort_acc  <- irr::kripp.alpha(t(gort_acc_t),  method = "interval")
ka_gort_flu  <- irr::kripp.alpha(t(gort_flu_t),  method = "interval")
ka_gort_comp <- irr::kripp.alpha(t(gort_comp_t), method = "interval")

cat(
  "\nGORT Krippendorff's alpha (interval):",
  "\n  Rate:        ", round(ka_gort_rate$value, 2),
  "\n  Accuracy:    ", round(ka_gort_acc$value, 2),
  "\n  Fluency:     ", round(ka_gort_flu$value, 2),
  "\n  Comprehension:", round(ka_gort_comp$value, 2), "\n"
)

```