# Agreement between two raters on 40 slides. R 4.3.3, base R.
d <- read.csv("scores.csv"); a <- d$grade_a; b <- d$grade_b; n <- nrow(d); K <- 5
P <- table(factor(a, levels = 0:4), factor(b, levels = 0:4)) / n; P <- matrix(P, K, K)
kappa <- function(w) { po <- sum(w * P); pe <- sum(w * outer(rowSums(P), colSums(P))); c(po, pe, (po - pe) / (1 - pe)) }
I <- row(P) - 1; J <- col(P) - 1
k0 <- kappa((I == J) * 1); se <- sqrt(k0[1] * (1 - k0[1]) / (n * (1 - k0[2])^2))
cat(sprintf("Grades: exact agreement %.1f%%; Cohen's kappa %.3f (approximate 95%% CI %.3f to %.3f)\n", 100 * k0[1], k0[3], k0[3] - 1.96 * se, k0[3] + 1.96 * se))
kq <- kappa(1 - (I - J)^2 / (K - 1)^2)
cat(sprintf("Grades: quadratic-weighted kappa %.3f; mean difference (B - A) %.2f grades\n", kq[3], mean(b - a)))
cat(sprintf("Grades: Spearman correlation %.3f\n", cor(a, b, method = "spearman")))
x <- cbind(d$area_a, d$area_b); k <- 2; gm <- mean(x)
msr <- k * sum((rowMeans(x) - gm)^2) / (n - 1); msc <- n * sum((colMeans(x) - gm)^2) / (k - 1)
sse <- sum((x - outer(rowMeans(x), rep(1, k)) - outer(rep(1, n), colMeans(x)) + gm)^2); mse <- sse / ((n - 1) * (k - 1))
icc_c <- (msr - mse) / (msr + (k - 1) * mse); icc_a <- (msr - mse) / (msr + (k - 1) * mse + k * (msc - mse) / n)
fl <- qf(0.975, n - 1, (n - 1) * (k - 1)); fu <- qf(0.975, (n - 1) * (k - 1), n - 1); f0 <- msr / mse
lr <- log(x[, 2] / x[, 1]); m <- mean(lr); sdl <- sd(lr)
cat(sprintf("Lesion area: ratio B/A %.3f (95%% limits of agreement %.3f to %.3f, on the log scale)\n", exp(m), exp(m - 1.96 * sdl), exp(m + 1.96 * sdl)))
cat(sprintf("Lesion area: Pearson r %.3f; mean difference (B - A) %.3f mm2\n", cor(x[, 1], x[, 2]), mean(x[, 2] - x[, 1])))
cat(sprintf("Lesion area: ICC consistency, single rater, ICC(3,1) = %.3f (95%% CI %.3f to %.3f)\n", icc_c, (f0 / fl - 1) / (f0 / fl + k - 1), (f0 * fu - 1) / (f0 * fu + k - 1)))
fj <- msc / mse; rr <- icc_a                         # CI for ICC(2,1) (Shrout & Fleiss / McGraw & Wong approximation)
v <- (k - 1) * (n - 1) * (k * rr * fj + n * (1 + (k - 1) * rr) - k * rr)^2 / ((n - 1) * k^2 * rr^2 * fj^2 + (n * (1 + (k - 1) * rr) - k * rr)^2)
f1 <- qf(0.975, n - 1, v); f2 <- qf(0.975, v, n - 1)
lo <- n * (msr - f1 * mse) / (f1 * (k * msc + (k * n - k - n) * mse) + n * msr)
hi <- n * (f2 * msr - mse) / (k * msc + (k * n - k - n) * mse + n * f2 * msr)
cat(sprintf("Lesion area: ICC absolute agreement, single rater, ICC(2,1) = %.3f (95%% CI %.3f to %.3f)\n", icc_a, lo, hi))