# Agreement between two raters on 40 slides. Python 3.13, NumPy, SciPy, pandas.
import numpy as np, pandas as pd
from scipy import stats

d = pd.read_csv("scores.csv")                        # columns: slide, grade_a, grade_b, area_a, area_b
a, b = d.grade_a.to_numpy(), d.grade_b.to_numpy(); n = len(d); K = 5
tab = np.zeros((K, K)); np.add.at(tab, (a, b), 1); P = tab / n
def kappa(w):                                        # w: agreement weights (1 on the diagonal)
    po = (w * P).sum(); pe = (w * np.outer(P.sum(1), P.sum(0))).sum()
    return po, pe, (po - pe) / (1 - pe)
i, j = np.indices((K, K))
po, pe, k0 = kappa((i == j).astype(float))
se = np.sqrt(po * (1 - po) / (n * (1 - pe) ** 2))  # Cohen's (1960) approximate standard error
print(f"Grades: exact agreement {po:.1%}; Cohen's kappa {k0:.3f} (approximate 95% CI {k0 - 1.96 * se:.3f} to {k0 + 1.96 * se:.3f})")
_, _, kq = kappa(1 - (i - j) ** 2 / (K - 1) ** 2)
print(f"Grades: quadratic-weighted kappa {kq:.3f}; mean difference (B - A) {np.mean(b - a):.2f} grades")
print(f"Grades: Spearman correlation {stats.spearmanr(a, b).statistic:.3f}")

x = np.column_stack([d.area_a, d.area_b]); k = 2
gm = x.mean(); msr = k * ((x.mean(1) - gm) ** 2).sum() / (n - 1); msc = n * ((x.mean(0) - gm) ** 2).sum() / (k - 1)
sse = ((x - x.mean(1, keepdims=True) - x.mean(0) + gm) ** 2).sum(); 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 = stats.f.ppf(0.975, n - 1, (n - 1) * (k - 1)); fu = stats.f.ppf(0.975, (n - 1) * (k - 1), n - 1); f0 = msr / mse
lr = np.log(x[:, 1] / x[:, 0]); m, sd = lr.mean(), lr.std(ddof=1)
print(f"Lesion area: ratio B/A {np.exp(m):.3f} (95% limits of agreement {np.exp(m - 1.96 * sd):.3f} to {np.exp(m + 1.96 * sd):.3f}, on the log scale)")
print(f"Lesion area: Pearson r {stats.pearsonr(x[:, 0], x[:, 1]).statistic:.3f}; mean difference (B - A) {np.mean(x[:, 1] - x[:, 0]):.3f} mm2")
print(f"Lesion area: ICC consistency, single rater, ICC(3,1) = {icc_c:.3f} (95% CI {(f0 / fl - 1) / (f0 / fl + k - 1):.3f} to {(f0 * fu - 1) / (f0 * fu + k - 1):.3f})")
fj = msc / mse; r = icc_a                            # CI for ICC(2,1) (Shrout & Fleiss / McGraw & Wong approximation)
v = (k - 1) * (n - 1) * (k * r * fj + n * (1 + (k - 1) * r) - k * r) ** 2 / ((n - 1) * k ** 2 * r ** 2 * fj ** 2 + (n * (1 + (k - 1) * r) - k * r) ** 2)
f1 = stats.f.ppf(0.975, n - 1, v); f2 = stats.f.ppf(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)
print(f"Lesion area: ICC absolute agreement, single rater, ICC(2,1) = {icc_a:.3f} (95% CI {lo:.3f} to {hi:.3f})")