Intraclass Correlation — Rater Reliability

suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))

compute_shared <- function(df, params, col_map = list()) {
  # === SHARED EXPORTS ===
  #   initial_rows/final_rows/rows_removed  $ row accounting (ratings)
  #   subj_h / rater_h / score_h  $ humanized names of the three mapped columns
  #   n_subjects / k_raters       $ complete-crossing design dimensions
  #   rater_levels                $ character — the actual rater level values
  #   n_duplicates                $ duplicate ratings averaged
  #   n_incomplete                $ subjects excluded (not rated by all raters)
  #   incomplete_subjects         $ character — the excluded subject ids
  #   icc / icc_p                 $ named list of ICC values + F-test p for ICC(3,1)
  #   icc_results_df   $ data.frame(icc_type, value, interpretation, use_when)
  #   rater_means_df   $ data.frame(rater_name, mean_score)
  #   agreement_df     $ data.frame(rater_1_score, rater_2_score) — complete subjects
  #   agreement_r      $ numeric — Pearson r between the first two raters
  #   variance_df      $ data.frame(source, variance_pct)
  #   headline_band    $ Koo & Li band for ICC(2,1)
  #   metrics / json_output
  # === /SHARED EXPORTS ===

  subj_h  <- humanize_semantic("subject", col_map)
  rater_h <- humanize_semantic("rater", col_map)
  score_h <- humanize_semantic("score", col_map)

initial_rows <- nrow(df)
  for (req in c("subject", "rater", "score")) {
    if (is.null(df[[req]])) {
      stop(sprintf("column_mapping must map the &#x27;%s' column (%s / %s / %s).",
                   req, subj_h, rater_h, score_h))
    }
  }

v <- df$score
  if (!is.numeric(v)) {
    conv <- suppressWarnings(as.numeric(as.character(v)))
    n_orig <- sum(!is.na(v) & as.character(v) != "")
    if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) {
      v <- conv
    } else {
      stop(sprintf("The rating column(%s) must be numeric — fewer than 95%% of its values convert cleanly.",
                   score_h))
    }
  }
  df$score   <- v
  df$subject <- trimws(as.character(df$subject))
  df$rater   <- trimws(as.character(df$rater))
  keep <- !is.na(df$score) & !is.na(df$subject) & df$subject != "" &
          !is.na(df$rater) & df$rater != ""
  n_invalid <- sum(!keep)
  df <- df[keep, , drop = FALSE]
  if (nrow(df) < 10) {
    stop(sprintf("Only %d usable ratings after removing rows with missing %s, %s, or %s — need at least 10.",
                 nrow(df), subj_h, rater_h, score_h))
  }

  rater_levels <- sort(unique(df$rater))
  k <- length(rater_levels)
  if (k < 2) {
    stop(sprintf("ICC needs at least 2 distinct values in the %s column; found %d. Each row should be one rating by one rater.",
                 rater_h, k))
  }
  if (k > 30) {
    stop(sprintf("The %s column has %d distinct values — that looks like an identifier, not a set of raters. Map the column naming who/what gave each rating.",
                 rater_h, k))
  }

agg <- aggregate(score ~ subject + rater, data = df, FUN = mean)
  n_duplicates <- nrow(df) - nrow(agg)

cnt <- table(agg$subject)
  complete_subjects   <- names(cnt)[cnt == k]
  incomplete_subjects <- names(cnt)[cnt < k]
  n_incomplete <- length(incomplete_subjects)
  n <- length(complete_subjects)
  if (n < 5) {
    stop(sprintf("ICC needs at least 5 %s values rated by every %s; only %d of %d have a complete set of %d ratings.",
                 subj_h, rater_h, n, length(cnt), k))
  }
  dat <- agg[agg$subject %in% complete_subjects, , drop = FALSE]
  if (is.na(var(dat$score)) || isTRUE(var(dat$score) == 0)) {
    stop(sprintf("The %s values are constant — reliability cannot be estimated.", score_h))
  }
  final_rows <- nrow(dat)
  rows_removed <- initial_rows - final_rows

  dat$subject_f <- factor(dat$subject)
  dat$rater_f   <- factor(dat$rater, levels = rater_levels)

a2  <- aov(score ~ subject_f + rater_f, data = dat)
  ms2 <- summary(a2)[[1]][["Mean Sq"]]
  MSR <- ms2[1]; MSC <- ms2[2]; MSE <- ms2[3]

  icc21 <- (MSR - MSE) / (MSR + (k - 1) * MSE + k * (MSC - MSE) / n)
  icc31 <- (MSR - MSE) / (MSR + (k - 1) * MSE)
  icc2k <- (MSR - MSE) / (MSR + (MSC - MSE) / n)
  icc3k <- (MSR - MSE) / MSR

a1  <- aov(score ~ subject_f, data = dat)
  ms1 <- summary(a1)[[1]][["Mean Sq"]]
  MSB <- ms1[1]; MSW <- ms1[2]
  icc11 <- (MSB - MSW) / (MSB + (k - 1) * MSW)

f_stat <- MSR / MSE
  f_df1  <- n - 1
  f_df2  <- (n - 1) * (k - 1)
  f_p    <- pf(f_stat, f_df1, f_df2, lower.tail = FALSE)

koo_li <- function(x) {
    if (is.na(x)) return("not estimable")
    if (x < 0.5) "poor" else if (x < 0.75) "moderate" else if (x < 0.9) "good" else "excellent"
  }
  headline_band <- koo_li(icc21)

  icc_results_df <- data.frame(
    icc_type = c("ICC(1,1) — one-way random, single rater",
                 "ICC(2,1) — absolute agreement, single rater",
                 "ICC(3,1) — consistency, single rater",
                 "ICC(2,k) — absolute agreement, average of all raters",
                 "ICC(3,k) — consistency, average of all raters"),
    value = round(c(icc11, icc21, icc31, icc2k, icc3k), 3),
    interpretation = sapply(c(icc11, icc21, icc31, icc2k, icc3k), koo_li),
    use_when = c(
      sprintf("Each %s is rated by a DIFFERENT random set of raters(no shared panel).", subj_h),
      sprintf("Raters are a random sample and absolute %s values must match — the default to report.", score_h),
      sprintf("These exact %d raters are the only ones of interest; systematic differences between them are forgiven.", k),
      sprintf("Decisions rest on the AVERAGE %s of all %d raters, and absolute values matter.", score_h, k),
      sprintf("Decisions rest on the average of these exact %d raters; only relative ordering matters.", k)
    ),
    stringsAsFactors = FALSE
  )

var_subject <- max(0, (MSR - MSE) / k)
  var_rater   <- max(0, (MSC - MSE) / n)
  var_error   <- max(0, MSE)
  var_total   <- var_subject + var_rater + var_error
  if (var_total <= 0) var_total <- 1
  variance_df <- data.frame(
    source = c("Between subjects", "Between raters", "Residual"),
    variance_pct = round(100 * c(var_subject, var_rater, var_error) / var_total, 1),
    stringsAsFactors = FALSE
  )

rm_agg <- aggregate(score ~ rater_f, data = dat, FUN = mean)
  rater_means_df <- data.frame(
    rater_name = as.character(rm_agg$rater_f),
    mean_score = round(rm_agg$score, 2),
    stringsAsFactors = FALSE
  )

ok_idx <- which(!is.na(rater_means_df$mean_score))
  hi_idx <- ok_idx[which.max(rater_means_df$mean_score[ok_idx])]
  lo_idx <- ok_idx[which.min(rater_means_df$mean_score[ok_idx])]
  rater_spread <- rater_means_df$mean_score[hi_idx] - rater_means_df$mean_score[lo_idx]

d1 <- dat[dat$rater == rater_levels[1], c("subject", "score")]
  d2 <- dat[dat$rater == rater_levels[2], c("subject", "score")]
  names(d1)[2] <- "rater_1_score"
  names(d2)[2] <- "rater_2_score"
  merged <- merge(d1, d2, by = "subject")
  merged <- merged[order(merged$rater_1_score), , drop = FALSE]
  set.seed(42)
  if (nrow(merged) > 1000) merged <- merged[sample(nrow(merged), 1000), , drop = FALSE]
  agreement_df <- data.frame(
    rater_1_score = round(merged$rater_1_score, 3),
    rater_2_score = round(merged$rater_2_score, 3),
    stringsAsFactors = FALSE
  )
  agreement_r <- if (nrow(agreement_df) >= 3)
    suppressWarnings(cor(agreement_df$rater_1_score, agreement_df$rater_2_score))
  else NA_real_

  icc <- list(icc11 = icc11, icc21 = icc21, icc31 = icc31,
              icc2k = icc2k, icc3k = icc3k)

  metrics <- list(
    `Subjects`     = n,
    `Raters`       = k,
    `Ratings Used` = final_rows,
    `ICC(2,1)`     = round(icc21, 3),
    `Reliability`  = headline_band,
    `F-test p`     = signif(f_p, 3)
  )

  json_output <- list(
    answer = paste0(
      "Intraclass correlation from ", format(final_rows, big.mark = ","),
      " ratings of ", n, " ", subj_h, " values by ", k, " raters(", rater_h,
      "): ICC(2,1) = ", round(icc21, 3), " — ", headline_band,
      " absolute agreement for a single rater(Koo & Li). Consistency ICC(3,1) = ",
      round(icc31, 3), "; average-measure ICC(2,k) = ", round(icc2k, 3),
      ". The F-test for subject discrimination is ",
      ifelse(is.na(f_p), "not estimable",
             ifelse(f_p < 0.05, paste0("significant(p = ", signif(f_p, 3), ")"),
                    paste0("not significant(p = ", signif(f_p, 3), ")"))),
      ". Between-", subj_h, " differences explain ", variance_df$variance_pct[1],
      "% of the variance in ", score_h, "."
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "icc_results",
        "rater_means", "subject_agreement", "variance_breakdown"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows = initial_rows, final_rows = final_rows,
    rows_removed = rows_removed, n_invalid = n_invalid,
    subj_h = subj_h, rater_h = rater_h, score_h = score_h,
    n_subjects = n, k_raters = k, rater_levels = rater_levels,
    n_duplicates = n_duplicates, n_incomplete = n_incomplete,
    incomplete_subjects = incomplete_subjects,
    icc = icc, f_stat = f_stat, f_df1 = f_df1, f_df2 = f_df2, f_p = f_p,
    icc_results_df = icc_results_df,
    rater_means_df = rater_means_df,
    rater_hi = rater_means_df$rater_name[hi_idx],
    rater_lo = rater_means_df$rater_name[lo_idx],
    rater_spread = rater_spread,
    agreement_df = agreement_df, agreement_r = agreement_r,
    variance_df = variance_df,
    headline_band = headline_band,
    metrics = metrics, json_output = json_output
  )
}