Intraclass Correlation (ICC)

# ICC(2,1): Inter-rater reliability for single rater from random effects model
# Two-way random effects: subjects are random, raters are random

# Install/load packages
library(irr)         # For ICC computation
library(psych)       # Alternative ICC with detailed output
library(tidyverse)   # Data manipulation
library(reshape2)    # Wide/long format conversion

# Simulate realistic HAM-D ratings (or load: data <- read.csv("hamd_ratings.csv"))
set.seed(2025)
n_subjects <- 45
n_raters <- 4

# True depression severity for each subject (0-52 scale)
true_severity <- rnorm(n_subjects, mean=24, sd=8)
true_severity <- pmax(0, pmin(52, true_severity))  # Bound to 0-52

# Each rater adds measurement error
data_wide <- data.frame(
  subject = 1:n_subjects,
  rater1 = round(true_severity + rnorm(n_subjects, 0, 3)),
  rater2 = round(true_severity + rnorm(n_subjects, 0, 3)),
  rater3 = round(true_severity + rnorm(n_subjects, 0, 3.5)),  # Slightly less reliable
  rater4 = round(true_severity + rnorm(n_subjects, 0, 2.5))   # Slightly more reliable
)

# Bound ratings to valid range
data_wide[,2:5] <- apply(data_wide[,2:5], 2, function(x) pmax(0, pmin(52, x)))

head(data_wide)

# === STEP 1: Check Assumptions ===

# 1. Normality of ratings per rater
par(mfrow=c(2,2))
for (i in 2:5) {
  qqnorm(data_wide[,i], main=paste("Rater", i-1))
  qqline(data_wide[,i])
}

for (i in 2:5) {
  test <- shapiro.test(data_wide[,i])
  cat(sprintf("Rater %d: W=%.3f, p=%.3f\n", i-1, test$statistic, test$p.value))
}

# 2. Homogeneity of variance across raters
library(car)
data_long <- melt(data_wide, id.vars="subject", variable.name="rater", value.name="rating")
leveneTest(rating ~ rater, data=data_long)
cat("\nStandard deviations by rater:\n")
data_wide %>% 
  summarise(across(starts_with("rater"), sd))

# 3. Check for systematic rater bias
cat("\nMean ratings by rater(check for bias):\n")
data_wide %>% 
  summarise(across(starts_with("rater"), mean))

# Boxplots
ggplot(data_long, aes(x=rater, y=rating, fill=rater)) +
  geom_boxplot() +
  labs(title="HAM-D Ratings by Rater(check for bias & variance)",
       y="HAM-D Score", x="Rater") +
  theme_classic()

# === STEP 2: Compute ICC(2,1) - Single Rater, Random Effects ===

# Method 1: Using irr package
icc_result <- icc(data_wide[,2:5], model="twoway", type="agreement", unit="single")
print(icc_result)
# Output: ICC(2,1) = 0.XX, 95% CI [0.XX, 0.XX], F-test p < .001

# Method 2: Using psych package (more detailed)
icc_psych <- ICC(data_wide[,2:5])
print(icc_psych)
# Shows ICC(1,1), ICC(2,1), ICC(3,1) and average measures versions

# === STEP 3: Interpretation Benchmarks ===
cat("\n=== ICC Interpretation Guidelines(Koo & Li, 2016) ===\n")
cat("< 0.50: Poor reliability\n")
cat("0.50-0.75: Moderate reliability\n")
cat("0.75-0.90: Good reliability\n")
cat("> 0.90: Excellent reliability\n\n")

icc_value <- icc_result$value
cat(sprintf("ICC(2,1) = %.3f\n", icc_value))
if (icc_value < 0.50) {
  cat("Interpretation: Poor single-rater reliability. Ratings highly variable.\n")
} else if (icc_value < 0.75) {
  cat("Interpretation: Moderate single-rater reliability. Consider averaging raters.\n")
} else if (icc_value < 0.90) {
  cat("Interpretation: Good single-rater reliability. Acceptable for research.\n")
} else {
  cat("Interpretation: Excellent single-rater reliability.\n")
}

# === STEP 4: Compare with Average Measures ICC(2,k) ===
icc_avg <- icc(data_wide[,2:5], model="twoway", type="agreement", unit="average")
cat(sprintf("\nICC(2,k) = %.3f [average of %d raters]\n", icc_avg$value, n_raters))
cat("ICC(2,k) is higher because averaging reduces measurement error(Spearman-Brown prophecy formula).\n")

# === STEP 5: Visualize Reliability ===

# Scatterplot matrix to visualize consistency
pairs(data_wide[,2:5], main="Scatterplot Matrix: Rater Consistency",
      pch=19, col=rgb(0,0,1,0.3))

# Bland-Altman style plot (Rater 1 vs Rater 2 example)
mean_12 <- (data_wide$rater1 + data_wide$rater2) / 2
diff_12 <- data_wide$rater1 - data_wide$rater2
mean_diff <- mean(diff_12)
sd_diff <- sd(diff_12)

ggplot(data.frame(mean=mean_12, diff=diff_12), aes(x=mean, y=diff)) +
  geom_point(alpha=0.5) +
  geom_hline(yintercept=mean_diff, color="blue", linetype="dashed") +
  geom_hline(yintercept=mean_diff + 1.96*sd_diff, color="red", linetype="dashed") +
  geom_hline(yintercept=mean_diff - 1.96*sd_diff, color="red", linetype="dashed") +
  labs(title="Bland-Altman Plot: Rater 1 vs Rater 2",
       x="Mean HAM-D Score", y="Difference(Rater1 - Rater2)",
       subtitle=paste("Mean diff:", round(mean_diff,2), "± 95% LOA")) +
  theme_classic()

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===\n")
cat(sprintf(
  "Inter-rater reliability for single clinician ratings of depression severity 
(HAM-D) was assessed using ICC(2,1) based on a two-way random effects model 
with absolute agreement. Forty-five patients were each rated by four randomly 
selected clinicians from a pool of 12 trained raters. The ICC was %.2f 
(95%% CI [%.2f, %.2f]), indicating %s single-rater reliability. The F-test 
confirmed significant between-subject variance, F(%d, %d) = %.2f, p < .001. 
When averaging across all four raters, reliability improved to ICC(2,4) = %.2f, 
suggesting that composite scores from multiple raters provide %s reliability 
for clinical and research use.\n",
  icc_result$value,
  icc_result$lbound,
  icc_result$ubound,
  ifelse(icc_result$value >= 0.75, "good", ifelse(icc_result$value >= 0.50, "moderate", "poor")),
  icc_result$df1,
  icc_result$df2,
  icc_result$Fvalue,
  icc_avg$value,
  ifelse(icc_avg$value >= 0.75, "good to excellent", "moderate to good")
))