Cohen's Kappa (κ)

# ============================================================================
# Comprehensive Cohen's Kappa Analysis for Diagnostic Agreement
# Complete workflow: simulation, kappa calculation, diagnostics, visualization
# ============================================================================

library(irr)          # For kappa calculations
library(psych)        # Additional statistics
library(vcd)          # For agreement visualization
library(ggplot2)      # For plotting
library(gridExtra)    # For multi-panel plots
library(boot)         # For bootstrap confidence intervals

# Seed for reproducibility
set.seed(2025)

# ============================================================================
# 1. DATA SIMULATION: 120 patients, 2 psychiatrists, 4 diagnostic categories
# ============================================================================

n <- 120
categories <- c("Depression", "Anxiety", "Bipolar", "No Disorder")

cat("\n========== COHEN'S KAPPA: DIAGNOSTIC AGREEMENT ==========\n")
cat("Sample size:", n, "patients\n")
cat("Raters: 2 psychiatrists\n")
cat("Categories:", length(categories), "diagnostic categories\n\n")

# Simulate true diagnoses (latent variable)
true_diagnosis <- sample(categories, n, replace=TRUE, 
                         prob=c(0.35, 0.25, 0.15, 0.25))

# Psychiatrist 1: 85% accuracy
rater1 <- sapply(true_diagnosis, function(x) {
  if(runif(1) < 0.85) return(x)
  sample(setdiff(categories, x), 1)
})

# Psychiatrist 2: 80% accuracy (slightly lower)
rater2 <- sapply(true_diagnosis, function(x) {
  if(runif(1) < 0.80) return(x)
  sample(setdiff(categories, x), 1)
})

# Create data frame
diag_data <- data.frame(
  Rater1 = factor(rater1, levels=categories),
  Rater2 = factor(rater2, levels=categories)
)

# ============================================================================
# 2. CALCULATE COHEN'S KAPPA
# ============================================================================

cat("--- Cohen's Kappa Calculation ---\n")
kappa_result <- kappa2(diag_data, weight="unweighted")
print(kappa_result)

kappa_value <- kappa_result$value
kappa_se <- sqrt(kappa_result$var.kappa)

cat("\nCohen's Kappa: κ =", round(kappa_value, 3))
cat("\nStandard Error:", round(kappa_se, 3))

# ============================================================================
# 3. CONFIDENCE INTERVALS
# ============================================================================

cat("\n\n--- 95% Confidence Interval ---\n")
# Asymptotic CI
ci_lower <- kappa_value - 1.96 * kappa_se
ci_upper <- kappa_value + 1.96 * kappa_se
cat("Asymptotic 95% CI: [", round(ci_lower, 3), ",", round(ci_upper, 3), "]\n")

# Bootstrap CI
cat("\nBootstrap 95% CI(500 iterations):\n")
boot_kappas <- replicate(500, {
  indices <- sample(1:n, replace=TRUE)
  boot_data <- diag_data[indices, ]
  kappa2(boot_data, weight="unweighted")$value
})
boot_ci <- quantile(boot_kappas, c(0.025, 0.975))
cat("Bootstrap CI: [", round(boot_ci[1], 3), ",", round(boot_ci[2], 3), "]\n")

# ============================================================================
# 4. STATISTICAL SIGNIFICANCE TEST
# ============================================================================

cat("\n--- Statistical Significance Test ---\n")
z_score <- kappa_value / kappa_se
p_value <- 2 * (1 - pnorm(abs(z_score)))  # Two-tailed
cat("H₀: κ = 0 (no agreement beyond chance)\n")
cat("Z-score:", round(z_score, 3), "\n")
cat("P-value:", format.pval(p_value, digits=3), "\n")
cat(ifelse(p_value < 0.001, "Result: Highly significant(p < .001)\n",
           ifelse(p_value < 0.05, "Result: Significant(p < .05)\n",
                  "Result: Not significant\n")))

# ============================================================================
# 5. CONFUSION MATRIX AND AGREEMENT STATISTICS
# ============================================================================

cat("\n--- Confusion Matrix ---\n")
conf_matrix <- table(diag_data$Rater1, diag_data$Rater2)
print(addmargins(conf_matrix))

cat("\n--- Agreement Statistics ---\n")
# Observed agreement
po <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("Observed agreement(Po):", round(po, 3), "(", 
    round(po*100, 1), "%)\n")

# Expected agreement by chance
rater1_marginal <- rowSums(conf_matrix) / sum(conf_matrix)
rater2_marginal <- colSums(conf_matrix) / sum(conf_matrix)
pe <- sum(rater1_marginal * rater2_marginal)
cat("Expected agreement(Pe):", round(pe, 3), "(", 
    round(pe*100, 1), "%)\n")
cat("Agreement beyond chance:", round((po - pe)*100, 1), "%\n")

# Kappa formula verification
kappa_manual <- (po - pe) / (1 - pe)
cat("\nKappa(manual): κ =", round(kappa_manual, 3), "\n")

# ============================================================================
# 6. PREVALENCE AND BIAS INDICES
# ============================================================================

cat("\n--- Prevalence and Bias Indices ---\n")
# Prevalence index
prevalence_idx <- abs(sum(conf_matrix[1,]) - sum(conf_matrix[,1])) / n
cat("Prevalence Index:", round(prevalence_idx, 3))
cat("\n(Measures asymmetry in marginal totals)\n")

# Bias index
bias_idx <- (sum(conf_matrix[1,]) - sum(conf_matrix[,1])) / n
cat("Bias Index:", round(bias_idx, 3))
cat("\n(Positive: Rater1 diagnoses more; Negative: Rater2 diagnoses more)\n")

# PABAK (Prevalence-Adjusted Bias-Adjusted Kappa)
pabak <- 2 * po - 1
cat("\nPABAK:", round(pabak, 3))
cat("\n(Less affected by prevalence than standard kappa)\n")

# ============================================================================
# 7. CATEGORY-SPECIFIC ANALYSIS
# ============================================================================

cat("\n--- Category-Specific Kappas ---\n")
for(cat_name in categories) {
  # Binary kappa for this category vs all others
  binary_r1 <- ifelse(rater1 == cat_name, cat_name, "Other")
  binary_r2 <- ifelse(rater2 == cat_name, cat_name, "Other")
  binary_data <- data.frame(
    R1 = factor(binary_r1, levels=c(cat_name, "Other")),
    R2 = factor(binary_r2, levels=c(cat_name, "Other"))
  )
  cat_kappa <- kappa2(binary_data)$value
  cat("  ", cat_name, ": κ =", round(cat_kappa, 3), "\n")
}

# ============================================================================
# 8. INTERPRETATION
# ============================================================================

cat("\n========== INTERPRETATION ==========\n")
if(kappa_value < 0) {
  interpretation <- "Poor(less than chance agreement)"
} else if(kappa_value < 0.20) {
  interpretation <- "Slight"
} else if(kappa_value < 0.40) {
  interpretation <- "Fair"
} else if(kappa_value < 0.60) {
  interpretation <- "Moderate"
} else if(kappa_value < 0.80) {
  interpretation <- "Substantial"
} else {
  interpretation <- "Almost Perfect"
}

cat("\nLandis & Koch(1977) Classification:", interpretation)
cat("\n\nThe two psychiatrists showed", tolower(interpretation), "agreement")
cat("\n(κ =", round(kappa_value, 3), ", 95% CI [", 
    round(ci_lower, 3), ",", round(ci_upper, 3), "]).\n")
cat("\nObserved agreement was", round(po*100, 1), "%, which is")
cat("\n", round((po-pe)*100, 1), "% better than expected by chance alone.\n")

if(kappa_value >= 0.60) {
  cat("\nThis level of agreement is acceptable for research purposes.\n")
  if(kappa_value >= 0.75) {
    cat("Agreement is sufficient for clinical decision-making.\n")
  }
} else {
  cat("\nThis level of agreement may be insufficient for clinical use.\n")
  cat("Consider additional rater training or clearer diagnostic criteria.\n")
}

cat("\n========== ANALYSIS COMPLETE ==========\n")