Kappa Cohena i ICC w SPSS – zgodność oceniających krok po kroku

* Zgodność oceniających: 80 raportów, oceniający A, B i C.
* Składnia zgodna z opisem CROSSTABS, WEIGHTED KAPPA i RELIABILITY w Command Syntax Reference; nie była uruchamiana w SPSS.

* 1. Kappa Cohena: główny problem raportu (kategorie nominalne).
CROSSTABS
  /TABLES=kat_A BY kat_B
  /STATISTICS=KAPPA
  /CELLS=COUNT TOTAL.

* 2. Kappa ważona: ocena z rubryki 1–4 (SPSS 27 lub nowszy), wagi liniowe i kwadratowe.
WEIGHTED KAPPA poz_A BY poz_B
  /CRITERIA ASYMPTOTICCILEVEL=95 WEIGHTING=LINEAR
  /PRINT TABULATION(ALL).
WEIGHTED KAPPA poz_A BY poz_B
  /CRITERIA ASYMPTOTICCILEVEL=95 WEIGHTING=QUADRATIC.

* 3. ICC: punkty 0–40, model dwuczynnikowy losowy, zgodność bezwzględna.
RELIABILITY
  /VARIABLES=pkt_A pkt_B pkt_C
  /SCALE(PUNKTY)=ALL
  /MODEL=ALPHA
  /ICC=MODEL(RANDOM) TYPE(ABSOLUTE) CIN=95 TESTVAL=0.

* Dla porównania: model mieszany, spójność.
RELIABILITY
  /VARIABLES=pkt_A pkt_B pkt_C
  /SCALE(PUNKTY)=ALL
  /MODEL=ALPHA
  /ICC=MODEL(MIXED) TYPE(CONSISTENCY) CIN=95 TESTVAL=0.

# B311 – kappa Cohena, kappa ważona i ICC w SPSS: liczby do tabel w układzie SPSS
# Dane symulowane, set.seed(3110): 80 raportów z ćwiczeń laboratoryjnych oceniło niezależnie troje
# prowadzących (A, B, C). A i B przypisali raportowi główny problem (kategorie nominalne) i ogólną ocenę
# z rubryki 1–4 (kategorie porządkowe); A, B i C przyznali punkty 0–40. Oceniający C punktuje surowiej.
# Metody jak w SPSS: kappa z tabeli kwadratowej z asymptotycznym błędem standardowym
# (Fleiss i in., 1969) i T = kappa / błąd przy H0; kappa ważona z wagami liniowymi (domyślne) i kwadratowymi;
# ICC ze średnich kwadratów dwuczynnikowej ANOVA (raport × oceniający) z testem F i przedziałami
# McGrawa i Wonga (1996); kontrola: psych::ICC() i DescTools::CohenKappa().
suppressPackageStartupMessages({ library(haven); library(psych); library(DescTools) })
set.seed(3110)
n <- 80
# 1. główny problem raportu: 1 = metoda, 2 = analiza danych, 3 = wnioski
prawda <- sample(1:3, n, replace = TRUE, prob = c(.40, .35, .25))
ocen_kat <- function(traf) ifelse(runif(n) < traf, prawda, sapply(prawda, function(k) sample(setdiff(1:3, k), 1)))
kat_A <- ocen_kat(.85); kat_B <- ocen_kat(.80)
# 2. ogólna ocena z rubryki: 1 = niewystarczający, 2 = dostateczny, 3 = dobry, 4 = bardzo dobry
jakosc <- rnorm(n)
progi <- c(-Inf, -0.9, 0.1, 1.0, Inf)
poz_A <- as.integer(cut(jakosc + rnorm(n, 0, .45), progi))
poz_B <- as.integer(cut(jakosc + rnorm(n, 0, .45), progi))
# 3. punkty 0–40 od trojga oceniających; C przyznaje średnio ok. 3 punkty mniej
wynik <- 26 + 5 * jakosc
ogr <- function(x) pmin(40, pmax(0, round(x)))
pkt_A <- ogr(wynik + rnorm(n, 0, 2)); pkt_B <- ogr(wynik + 0.5 + rnorm(n, 0, 2)); pkt_C <- ogr(wynik - 3 + rnorm(n, 0, 2))
dane <- data.frame(id = 1:n, kat_A, kat_B, poz_A, poz_B, pkt_A, pkt_B, pkt_C)

# Krok 1. Plik .sav z etykietami zmiennych i wartości – i ten sam zbiór w CSV
sav <- dane
et_kat <- c("metoda" = 1, "analiza danych" = 2, "wnioski" = 3)
et_poz <- c("niewystarczający" = 1, "dostateczny" = 2, "dobry" = 3, "bardzo dobry" = 4)
for (v in c("kat_A", "kat_B")) sav[[v]] <- labelled(sav[[v]], labels = et_kat)
for (v in c("poz_A", "poz_B")) sav[[v]] <- labelled(sav[[v]], labels = et_poz)
etykiety <- c(id = "Numer raportu", kat_A = "Główny problem – oceniający A", kat_B = "Główny problem – oceniający B",
              poz_A = "Ocena z rubryki 1–4 – oceniający A", poz_B = "Ocena z rubryki 1–4 – oceniający B",
              pkt_A = "Punkty 0–40 – oceniający A", pkt_B = "Punkty 0–40 – oceniający B", pkt_C = "Punkty 0–40 – oceniający C")
for (v in names(etykiety)) attr(sav[[v]], "label") <- etykiety[[v]]
for (v in names(sav)) attr(sav[[v]], "format.spss") <- "F2.0"
write_sav(sav, "../kappa-icc-spss-dane.sav")
write.csv(dane, "../kappa-icc-spss-dane.csv", row.names = FALSE, fileEncoding = "UTF-8")
sapply(read_sav("../kappa-icc-spss-dane.sav"), function(x) names(attributes(x)))   # tylko etykiety i format

# Funkcja: kappa (ważona albo nie) z błędami standardowymi Fleissa, Cohena i Everitta (1969)
kappa_spss <- function(a, b, k, wagi = c("brak", "liniowe", "kwadratowe")) {
  wagi <- match.arg(wagi)
  tab <- table(factor(a, 1:k), factor(b, 1:k)); N <- sum(tab); p <- tab / N
  w <- switch(wagi, brak = diag(k),
              liniowe = 1 - abs(outer(1:k, 1:k, "-")) / (k - 1),
              kwadratowe = 1 - (outer(1:k, 1:k, "-"))^2 / (k - 1)^2)
  pr <- rowSums(p); pc <- colSums(p)
  po <- sum(w * p); pe <- sum(w * outer(pr, pc)); kap <- (po - pe) / (1 - pe)
  wr <- as.vector(w %*% pc); wc <- as.vector(t(w) %*% pr)      # średnie wagi wiersza i kolumny
  sk <- outer(wr, wc, "+")
  var1 <- (sum(p * (w * (1 - pe) - sk * (1 - po))^2) - (po * pe - 2 * pe + po)^2) / (N * (1 - pe)^4)
  var0 <- (sum(outer(pr, pc) * (w - sk)^2) - pe^2) / (N * (1 - pe)^2)
  se <- sqrt(var1); se0 <- sqrt(var0); z <- kap / se0
  c(N = N, po = po, pe = pe, kappa = kap, ase = se, T_z = z, p = 2 * pnorm(-abs(z)),
    ci_dol = kap - qnorm(.975) * se, ci_gor = kap + qnorm(.975) * se)
}

# Krok 2. Crosstabs: kat_A × kat_B – tabela krzyżowa i Symmetric Measures (Kappa)
table(kat_A = dane$kat_A, kat_B = dane$kat_B)
k1 <- kappa_spss(dane$kat_A, dane$kat_B, 3)
round(k1, 3)
format.pval(k1[["p"]], digits = 3)
print(CohenKappa(dane$kat_A, dane$kat_B, conf.level = .95), digits = 4)   # kontrola: ten sam błąd standardowy

# Krok 3. Weighted Kappa: poz_A i poz_B, wagi liniowe (domyślne w SPSS) i kwadratowe
table(poz_A = dane$poz_A, poz_B = dane$poz_B)
kw <- rbind(nieważona = kappa_spss(dane$poz_A, dane$poz_B, 4),
            liniowe = kappa_spss(dane$poz_A, dane$poz_B, 4, "liniowe"),
            kwadratowe = kappa_spss(dane$poz_A, dane$poz_B, 4, "kwadratowe"))
round(kw, 3)
format.pval(kw[, "p"], digits = 3)
tab_poz <- table(factor(dane$poz_A, 1:4), factor(dane$poz_B, 1:4))
print(CohenKappa(tab_poz, weights = "Equal-Spacing", conf.level = .95), digits = 4)   # kontrola: wagi liniowe
print(CohenKappa(tab_poz, weights = "Fleiss-Cohen", conf.level = .95), digits = 4)    # kontrola: wagi kwadratowe

# Krok 4. ICC: średnie kwadraty z ANOVA raport × oceniający (pkt_A, pkt_B, pkt_C)
X <- as.matrix(dane[, c("pkt_A", "pkt_B", "pkt_C")]); k <- ncol(X)
round(colMeans(X), 2)
sr <- mean(X)
SSR <- k * sum((rowMeans(X) - sr)^2); SSC <- n * sum((colMeans(X) - sr)^2)
SST <- sum((X - sr)^2); SSE <- SST - SSR - SSC; SSW <- SST - SSR
MSR <- SSR / (n - 1); MSC <- SSC / (k - 1); MSE <- SSE / ((n - 1) * (k - 1)); MSW <- SSW / (n * (k - 1))
round(c(MSR = MSR, MSC = MSC, MSE = MSE, MSW = MSW), 3)

# Wzory McGrawa i Wonga (1996) – miara pojedyncza i średnia, przedziały 95%, test F przy wartości testowej 0
icc_spss <- function(model = c("ONEWAY", "RANDOM", "MIXED"), typ = c("CONSISTENCY", "ABSOLUTE"), alfa = .05) {
  model <- match.arg(model); typ <- match.arg(typ)
  if (model == "ONEWAY") {
    F0 <- MSR / MSW; df1 <- n - 1; df2 <- n * (k - 1)
    FL <- F0 / qf(1 - alfa / 2, df1, df2); FU <- F0 * qf(1 - alfa / 2, df2, df1)
    poj <- c((MSR - MSW) / (MSR + (k - 1) * MSW), (FL - 1) / (FL + k - 1), (FU - 1) / (FU + k - 1))
    sre <- c((MSR - MSW) / MSR, 1 - 1 / FL, 1 - 1 / FU)
  } else if (typ == "CONSISTENCY") {
    F0 <- MSR / MSE; df1 <- n - 1; df2 <- (n - 1) * (k - 1)
    FL <- F0 / qf(1 - alfa / 2, df1, df2); FU <- F0 * qf(1 - alfa / 2, df2, df1)
    poj <- c((MSR - MSE) / (MSR + (k - 1) * MSE), (FL - 1) / (FL + k - 1), (FU - 1) / (FU + k - 1))
    sre <- c((MSR - MSE) / MSR, 1 - 1 / FL, 1 - 1 / FU)
  } else {
    F0 <- MSR / MSE; df1 <- n - 1; df2 <- (n - 1) * (k - 1)
    r1 <- (MSR - MSE) / (MSR + (k - 1) * MSE + k * (MSC - MSE) / n)
    a <- k * r1 / (n * (1 - r1)); b <- 1 + k * r1 * (n - 1) / (n * (1 - r1))
    v <- (a * MSC + b * MSE)^2 / ((a * MSC)^2 / (k - 1) + (b * MSE)^2 / ((n - 1) * (k - 1)))
    Fs <- qf(1 - alfa / 2, n - 1, v); Fs2 <- qf(1 - alfa / 2, v, n - 1)
    lo <- n * (MSR - Fs * MSE) / (Fs * (k * MSC + (k * n - k - n) * MSE) + n * MSR)
    up <- n * (Fs2 * MSR - MSE) / (k * MSC + (k * n - k - n) * MSE + n * Fs2 * MSR)
    poj <- c(r1, lo, up)
    sre <- c((MSR - MSE) / (MSR + (MSC - MSE) / n), k * lo / (1 + (k - 1) * lo), k * up / (1 + (k - 1) * up))
  }
  out <- rbind(single = c(poj, F0, df1, df2), average = c(sre, F0, df1, df2))
  colnames(out) <- c("ICC", "lower", "upper", "F", "df1", "df2")
  cbind(round(out, 3), p = signif(pf(F0, df1, df2, lower.tail = FALSE), 3))
}
cat("\nTwo-Way Random, Absolute Agreement\n"); print(icc_spss("RANDOM", "ABSOLUTE"))
cat("\nTwo-Way Mixed, Consistency\n");        print(icc_spss("MIXED", "CONSISTENCY"))
cat("\nOne-Way Random\n");                    print(icc_spss("ONEWAY"))
# kontrola: psych::ICC (ICC1 = One-Way, ICC2 = absolutna zgodność, ICC3 = spójność; k = miara średnia)
print(psych::ICC(X, lmer = FALSE)$results[, c("type", "ICC", "F", "df1", "df2", "p", "lower bound", "upper bound")], digits = 3)

# Krok 5. Te same oceny bez oceniającego C: przesunięcie C obniżało zgodność bezwzględną
X <- X[, 1:2]; k <- 2; sr <- mean(X)
SSR <- k * sum((rowMeans(X) - sr)^2); SSC <- n * sum((colMeans(X) - sr)^2)
SSE <- sum((X - sr)^2) - SSR - SSC
MSR <- SSR / (n - 1); MSC <- SSC / (k - 1); MSE <- SSE / ((n - 1) * (k - 1))
cat("\nA i B: Two-Way Random, Absolute Agreement\n"); print(icc_spss("RANDOM", "ABSOLUTE"))

$id
[1] "label"       "format.spss"

$kat_A
[1] "label"       "format.spss" "class"       "labels"     

$kat_B
[1] "label"       "format.spss" "class"       "labels"     

$poz_A
[1] "label"       "format.spss" "class"       "labels"     

$poz_B
[1] "label"       "format.spss" "class"       "labels"     

$pkt_A
[1] "label"       "format.spss"

$pkt_B
[1] "label"       "format.spss"

$pkt_C
[1] "label"       "format.spss"

     kat_B
kat_A  1  2  3
    1 18  4  5
    2  4 21  4
    3  7  4 13
     N     po     pe  kappa    ase    T_z      p ci_dol ci_gor 
80.000  0.650  0.336  0.473  0.080  5.970  0.000  0.316  0.629 
[1] "2.37e-09"
 kappa lwr.ci upr.ci 
0.4727 0.3164 0.6289 
     poz_B
poz_A  1  2  3  4
    1 10  7  1  0
    2  2 16  5  0
    3  2  6  6  1
    4  0  1  8 15
            N    po    pe kappa   ase   T_z p ci_dol ci_gor
nieważona  80 0.588 0.254 0.447 0.072 7.011 0  0.305  0.589
liniowe    80 0.846 0.600 0.615 0.060 7.852 0  0.498  0.731
kwadratowe 80 0.937 0.744 0.756 0.051 6.840 0  0.657  0.855
[1] "2.36e-12" "4.09e-15" "7.94e-12"
 kappa lwr.ci upr.ci 
0.6146 0.4978 0.7313 
 kappa lwr.ci upr.ci 
0.7558 0.6567 0.8548 
pkt_A pkt_B pkt_C 
26.41 26.55 23.21 
    MSR     MSC     MSE     MSW 
 90.184 285.304   3.823   7.342 

Two-Way Random, Absolute Agreement
          ICC lower upper      F df1 df2        p
single  0.797 0.508 0.901 23.589  79 158 2.36e-58
average 0.922 0.756 0.965 23.589  79 158 2.36e-58

Two-Way Mixed, Consistency
          ICC lower upper      F df1 df2        p
single  0.883 0.836 0.919 23.589  79 158 2.36e-58
average 0.958 0.939 0.971 23.589  79 158 2.36e-58

One-Way Random
          ICC lower upper      F df1 df2        p
single  0.790 0.714 0.852 12.284  79 160 7.84e-40
average 0.919 0.882 0.945 12.284  79 160 7.84e-40
                         type   ICC    F df1 df2        p lower bound upper bound
Single_raters_absolute   ICC1 0.790 12.3  79 160 7.84e-40       0.714       0.852
Single_random_raters     ICC2 0.797 23.6  79 158 2.36e-58       0.508       0.901
Single_fixed_raters      ICC3 0.883 23.6  79 158 2.36e-58       0.836       0.919
Average_raters_absolute ICC1k 0.919 12.3  79 160 7.84e-40       0.882       0.945
Average_random_raters   ICC2k 0.922 23.6  79 158 2.36e-58       0.756       0.965
Average_fixed_raters    ICC3k 0.958 23.6  79 158 2.36e-58       0.939       0.971

A i B: Two-Way Random, Absolute Agreement
          ICC lower upper     F df1 df2        p
single  0.882 0.822 0.923 15.77  79  79 1.15e-27
average 0.937 0.902 0.960 15.77  79  79 1.15e-27