--- title: "grassr: rater reliability on binary outcomes, from rating matrix to Report Card" author: - Austin Semmel - Rachel Gidaro date: "2026-07-12" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{grassr: rater reliability on binary outcomes, from rating matrix to Report Card} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ## The data: a rater panel A **rater panel** is `k` raters who each score the same `N` subjects on one binary question. The subjects are whatever is being rated. The raters are whoever rates them, human or machine. Each rater sees every subject exactly once and answers the same yes-or-no question about it, so the study produces `N x k` binary calls and nothing else. Chest films read for a finding, adverse events judged preventable or not, transcripts coded for a symptom mention: different fields, same object. That object is the package's entire input, laid out as a matrix. One row per subject. One column per rater. Each cell is that rater's call on that subject: 1 (or `TRUE`, or the positive level of a two-level factor) means present, 0 means absent. The same shape carries two different questions. - **Inter-rater agreement** asks whether different raters agree with each other. Columns are different people; this is the default and everything below is this case. - **Intra-rater agreement** asks whether one rater agrees with herself on re-reading. Columns are one person's repeated viewings of the same subjects, and the call is `grass_report(Y, axis = "intra")`; the printed card has the same fields. What the functions accept, completely: - a `matrix` or `data.frame`, subjects in rows, raters in columns, at least 2 raters and one row per subject; - rows are independent subjects. Every row must be a different subject, never a repetition of one (repeated-measures designs are addressed just below); - cells as numeric 0/1, `logical`, or a two-level `factor` (character columns are rejected; convert to factor first); - column names are optional and become rater labels on the card; - no `NA`s: the functions stop rather than silently drop, so filter first with `Y <- Y[complete.cases(Y), ]` and report the exclusions in your Methods; - binary outcomes only. Ordinal scales and continuous measurements are outside this package's calibration. The calibration assumes every row is a new subject. Two common clinical designs break that assumption. - **Repetitions within a session.** A patient performs the task three times and the raters score every repetition. Do not stack the repetitions as rows: repetitions of one patient are correlated, the stacked matrix overstates the number of subjects, and the band comes out narrower than the design earns. Build one matrix per repetition and compare the cards, or fix one repetition per patient by protocol. - **The same patients across sessions over time.** Within a session the rows are distinct patients, so each session is a valid panel. Score one card per session; the sequence of cards is the stability record. grassr does not pool sessions into one number and does not split disagreement into patient, session, and rater components. Those are variance-component analyses, outside this package's calibration. A 0/1 matrix does not say which rows share a patient, so grassr cannot detect either design. You enforce this part of the contract, not an error message. In practice you load your scores from a file. Here we simulate the panel instead, so the truth is known and every claim below can be checked against it: 150 subjects at 30% prevalence, three raters who each call 87% of cases correctly. ```r library(grassr) set.seed(7) truth <- rbinom(150, 1, 0.30) # the latent true classes Y0 <- sapply(1:3, function(j) # three raters, 87% accurate ifelse(truth == 1, rbinom(150, 1, 0.87), rbinom(150, 1, 0.13))) colnames(Y0) <- c("R1", "R2", "R3") head(Y0, 8) #> R1 R2 R3 #> [1,] 1 0 1 #> [2,] 0 0 0 #> [3,] 0 0 1 #> [4,] 0 0 0 #> [5,] 0 0 0 #> [6,] 1 1 1 #> [7,] 0 0 0 #> [8,] 1 1 0 ``` Row 6 is unanimous. Rows 1, 3, and 8 split. The question every reliability study asks is whether these raters agree well enough to trust the scores. The conventional answer computes one coefficient and stamps a label on it from a fixed scale: 0.21 to 0.40 is *fair*, 0.41 to 0.60 *moderate*, 0.61 to 0.80 *substantial*. ## The problem: the label tracks the study, not the raters Every chance-corrected coefficient is the same fraction, written in plain terms as ``` coefficient = (observed agreement − chance agreement) / (1 − chance agreement) ``` The catch is the chance term. It depends on how often the finding appears in front of the raters, so the coefficient moves when prevalence moves, with the raters unchanged. At two raters the expected value has a closed form, and four lines of code show the problem: ```r q <- 0.87 # every rater calls 87% of cases correctly prev <- seq(0.02, 0.98, 0.005) Pa <- q^2 + (1 - q)^2 # P(two such raters agree); prevalence-free pip <- prev * q + (1 - prev) * (1 - q) # observed positive rate kappa <- (Pa - (pip^2 + (1 - pip)^2)) / (1 - (pip^2 + (1 - pip)^2)) op <- par(mar = c(4.2, 4.2, 2.4, 1)) plot(prev, kappa, type = "l", lwd = 2, ylim = c(0, 0.8), xlab = "prevalence of the finding", ylab = "expected Cohen's kappa", main = "Same raters (q = 0.87) at every prevalence") abline(h = c(0.20, 0.40, 0.60), lty = 3, col = "grey55") text(rep(0.99, 4), c(0.10, 0.30, 0.50, 0.70), adj = 1, col = "grey40", cex = 0.8, labels = c("slight", "fair", "moderate", "substantial")) pts <- c(0.05, 0.50, 0.90) ptk <- sapply(pts, function(p) { pp <- p * q + (1 - p) * (1 - q) (Pa - (pp^2 + (1 - pp)^2)) / (1 - (pp^2 + (1 - pp)^2)) }) points(pts, ptk, pch = 19) text(pts, ptk + 0.05, sprintf("%.2f", ptk), cex = 0.85) ``` ![plot of chunk hook-figure](grassr-hook-figure-1.png) ```r par(op) ``` One rater panel, three published verdicts. A screening study at 5% prevalence reports *slight* agreement, a balanced study reports *moderate*, a high-prevalence study reports *fair*, and the raters never changed. Rater count and sample size move the value too. The label tracks the study design, not the people. The three coefficients in common use split the same way on a single dataset, because each computes the chance term differently. Fleiss' kappa computes it from the panel's marginal positive rate, PABAK fixes it at one half, and AC1 interpolates between the two. Same data, three numbers, three labels. ## What grassr reports instead grassr asks where the observed coefficient sits among the values a panel could produce at this design. Four terms carry the whole method. - **Panel quality `q`** is the rate at which a rater calls a case correctly (sensitivity = specificity = `q` under the calibration model). This is the thing we actually want to know. - The **reference surface** holds the coefficient values a quality-`q` panel produces at the design (prevalence, raters `k`, subjects `N`). It is derived from a data-generating process, not stipulated. - The **percentile** is where the observed coefficient falls within everything the design can produce. This is the headline number. - The **consistency band** holds the qualities `q` consistent with the observation (95%). Its width is the design's resolving power. A small study cannot produce a narrow band, and the reader sees that directly. A fifth term, the spread `delta_hat`, measures how much the family's coefficients disagree about `q`. When the spread is small against chance, one number is a fair summary. When it is large, no single number is honest and the card decomposes the panel rater by rater. The whole workflow is one call on the matrix from the opening section: ```r grass_report(ratings = Y0) #> GRASS Report Card #> #> sample = 3 raters, N = 150, pi_hat = 0.36 #> PABAK = 0.52 -> 70th percentile | quality 0.80-0.90 <- primary #> AC1 = 0.55 -> 70th percentile | quality 0.80-0.90 #> Fleiss kappa = 0.48 -> 71st percentile | quality 0.80-0.90 #> ICC = 0.62 -> 72nd percentile | quality 0.81-0.90 [distribution-sensitive] #> read: this panel agreed more tightly than 70% of what panels at this design can produce; the data are consistent with panel quality 0.80-0.90. #> delta = 0 pp implied-quality spread (aligned) #> matched null = (k=3, N=150, q=0.85): delta_hat at the 10.1 percentile [design snapped] #> #> Notes: #> - N=150 clamped to nearest sim-grid N=100. #> - delta_hat is the implied-quality spread over the agreement family (PABAK, mean AC1, Fleiss kappa). ICC is reported on the panel but does not enter delta_hat (v0.5.0 scope: ICC's reference depends on the full subject-prevalence distribution F and does not share the (q, pi_+) sufficient statistic the agreement family does). #> - flag from delta_hat's percentile on the matched null (k=3, N=150, q=0.85; 50,000 draws); design snapped to nearest calibrated cell. #> #> See `summary(...)` for full panel and CI details. #> See `plot(...)` for a surface-position visualization. ``` Reading the card: - **`sample`** is the study: rater count `k`, subject count `N`, and the observed positive rate. - **The coefficient lines** show each observed value, its percentile, and its consistency band. The `read:` line restates the headline in plain language. There is no bootstrap qualifier and no labeled band between the percentile and the reader; the band width is the uncertainty statement. - **`delta`** is the implied-quality spread across the agreement family (PABAK, mean AC1, Fleiss kappa; ICC excluded by construction). Each coefficient implies a panel quality, and the spread is judged by its percentile on a bundled null distribution at the matched design: **caution at the 95th percentile, divergent at the 99th**, with mid-p at ties. A panel just under a cut should be read by its printed percentile, not the flag word alone. At two raters the spread is `not_applicable`, because PABAK and AC1 imply identical quality by construction. Every panel below is synthetic, so the truth is planted and every claim the card makes can be checked against it. The cards print exactly as `grass_report()` returns them. ## Case studies ### Three labels, one position Five raters score 1,000 subjects. Every rater operates at quality 0.87 and the finding has true prevalence 0.15. ```r gen_flat <- function(N, prev, Se, Sp, k, seed) { set.seed(seed) C <- rbinom(N, 1L, prev) sapply(seq_len(k), function(j) ifelse(C == 1L, rbinom(N, 1L, Se), rbinom(N, 1L, 1 - Sp))) } Y1 <- gen_flat(1000, prev = 0.15, Se = 0.87, Sp = 0.87, k = 5, seed = 61) ``` The three coefficients draw three different fixed labels from this one panel: Fleiss' kappa 0.36 reads *fair*, PABAK 0.54 reads *moderate*, AC1 0.64 reads *substantial*. A reader who saw only the labels would conclude they describe three different panels. ```r grass_report(Y1, verbose = FALSE) #> GRASS Report Card #> #> sample = 5 raters, N = 1000, pi_hat = 0.24 #> PABAK = 0.54 -> 74th percentile | quality 0.85-0.90 <- primary #> AC1 = 0.64 -> 74th percentile | quality 0.85-0.90 #> Fleiss kappa = 0.36 -> 74th percentile | quality 0.85-0.90 #> ICC = 0.48 -> 63rd percentile | quality 0.80-0.85 [distribution-sensitive] #> read: this panel agreed more tightly than 74% of what panels at this design can produce; the data are consistent with panel quality 0.85-0.90. #> delta = 0 pp implied-quality spread (aligned) #> matched null = (k=5, N=1000, q=0.85): delta_hat at the 42.2 percentile [design snapped] #> #> Notes: #> - Consistency band narrower than the calibrated q-grid spacing; endpoints interpolated within one grid gap. #> - delta_hat is the implied-quality spread over the agreement family (PABAK, mean AC1, Fleiss kappa). ICC is reported on the panel but does not enter delta_hat (v0.5.0 scope: ICC's reference depends on the full subject-prevalence distribution F and does not share the (q, pi_+) sufficient statistic the agreement family does). #> - flag from delta_hat's percentile on the matched null (k=5, N=1000, q=0.85; 50,000 draws); design snapped to nearest calibrated cell. #> #> See `summary(...)` for full panel and CI details. #> See `plot(...)` for a surface-position visualization. ``` The card gives all three one reading. Each sits at the 74th percentile of what this design can produce, and each carries the same band, quality 0.85 to 0.90, which contains the true 0.87. The label disagreement was an artifact of reading three prevalence-sensitive statistics against one fixed scale. ### Two panels the fixed bands call identical Two five-rater panels, both at true prevalence 0.30 with 1,000 subjects, both returning Fleiss kappa 0.36. The fixed scale gives one verdict, *fair*, for both. The panels could hardly be more different: in the first, three raters operate near-ideally and two have a severe specificity problem; in the second, all five raters share one operating point. ```r gen_logitnormal <- function(seed, k, N, Se, Sp, pi, F_sigma2 = 0.25) { set.seed(seed) p_i <- plogis(rnorm(N, qlogis(pi), sqrt(F_sigma2))) C <- rbinom(N, 1L, p_i) Y <- matrix(0L, N, k) for (j in seq_len(k)) Y[, j] <- rbinom(N, 1L, ifelse(C == 1L, Se[j], 1 - Sp[j])) Y } Y_differ <- gen_logitnormal(51L, 5L, 1000L, Se = rep(0.92, 5), Sp = c(0.92, 0.92, 0.92, 0.51, 0.51), pi = 0.30) Y_agree <- gen_logitnormal(151L, 5L, 1000L, Se = rep(0.95, 5), Sp = rep(0.70, 5), pi = 0.30) ``` ![plot of chunk case2-scatter](grassr-case2-scatter-1.png) The first panel splits across the low-specificity zone, and the card routes to per-rater output. ```r grass_report(Y_differ, bootstrap_B = 200, verbose = FALSE) #> GRASS Report Card #> #> sample = 5 raters, N = 1000, pi_hat = 0.47 #> PABAK = 0.36 -> 58th percentile <- primary #> AC1 = 0.38 -> 60th percentile #> Fleiss kappa = 0.36 -> 58th percentile #> ICC = 0.46 -> 58th percentile [distribution-sensitive] #> panel-agg. = suppressed (divergent) #> delta = 0.42 pp (divergent) #> matched null = (k=5, N=1000, q=0.85): delta_hat at the 99.5 percentile [design snapped] #> #> pairwise PABAK / surface percentile (lower / upper): #> R1 R2 R3 R4 R5 #> R1 -- 84% 84% 48% 40% #> R2 0.71 -- 81% 45% 41% #> R3 0.70 0.67 -- 45% 41% #> R4 0.27 0.24 0.25 -- 35% #> R5 0.21 0.22 0.21 0.17 -- #> #> per-rater vs panel-majority of OTHER raters (pooled-reference): #> R1 Se_tilde = 0.84 Sp_tilde = 0.93 (n_pos = 360, n_neg = 457, excl = 183) #> R2 Se_tilde = 0.85 Sp_tilde = 0.91 (n_pos = 352, n_neg = 460, excl = 188) #> R3 Se_tilde = 0.83 Sp_tilde = 0.92 (n_pos = 351, n_neg = 459, excl = 190) #> R4 Se_tilde = 0.92 Sp_tilde = 0.51 (n_pos = 327, n_neg = 581, excl = 92) #> R5 Se_tilde = 0.88 Sp_tilde = 0.50 (n_pos = 330, n_neg = 579, excl = 91) #> #> per-rater (latent-class fit; alongside pairwise): #> R1 Se = 0.93 (0.89, 0.96) Sp = 0.93 (0.92, 0.96) #> R2 Se = 0.92 (0.89, 0.95) Sp = 0.91 (0.88, 0.93) #> R3 Se = 0.89 (0.86, 0.93) Sp = 0.91 (0.88, 0.94) #> R4 Se = 0.93 (0.90, 0.96) Sp = 0.51 (0.47, 0.55) #> R5 Se = 0.89 (0.86, 0.92) Sp = 0.50 (0.46, 0.54) #> #> Notes: #> - delta_hat is the implied-quality spread over the agreement family (PABAK, mean AC1, Fleiss kappa). ICC is reported on the panel but does not enter delta_hat (v0.5.0 scope: ICC's reference depends on the full subject-prevalence distribution F and does not share the (q, pi_+) sufficient statistic the agreement family does). #> - flag from delta_hat's percentile on the matched null (k=5, N=1000, q=0.85; 50,000 draws); design snapped to nearest calibrated cell. #> #> See `summary(...)` for full panel and CI details. #> See `plot(...)` for a surface-position visualization. ``` The observed coefficients are unremarkable. The implied-quality spread is not: 0.42 quality points, at the 99.5th percentile of what chance produces at this design, so the card declines to average the panel into one number. The latent-class fit names the structure, R4 and R5 near specificity 0.51 against 0.91 to 0.93 for the rest. This panel does not need more raters or more subjects. It needs R4 and R5 retrained on the negative-call criterion. The second panel is the honest limit of any internal diagnostic. ```r grass_report(Y_agree, verbose = FALSE) #> GRASS Report Card #> #> sample = 5 raters, N = 1000, pi_hat = 0.49 #> PABAK = 0.36 -> 56th percentile | quality 0.75-0.85 <- primary #> AC1 = 0.36 -> 56th percentile | quality 0.75-0.85 #> Fleiss kappa = 0.36 -> 56th percentile | quality 0.75-0.85 #> ICC = 0.47 -> 61st percentile | quality 0.76-0.85 [distribution-sensitive] #> read: this panel agreed more tightly than 56% of what panels at this design can produce; the data are consistent with panel quality 0.75-0.85. #> delta = 0 pp implied-quality spread (aligned) #> matched null = (k=5, N=1000, q=0.75): delta_hat at the 35.6 percentile [design snapped] #> #> Notes: #> - delta_hat is the implied-quality spread over the agreement family (PABAK, mean AC1, Fleiss kappa). ICC is reported on the panel but does not enter delta_hat (v0.5.0 scope: ICC's reference depends on the full subject-prevalence distribution F and does not share the (q, pi_+) sufficient statistic the agreement family does). #> - flag from delta_hat's percentile on the matched null (k=5, N=1000, q=0.75; 50,000 draws); design snapped to nearest calibrated cell. #> #> See `summary(...)` for full panel and CI details. #> See `plot(...)` for a surface-position visualization. ``` Nothing in the card objects, yet the panel calls about half its subjects positive against a true prevalence of 0.30. Every rater shares the same specificity deficit, so the shared error shifts all marginals in unison and no internal comparison can see it. Agreement diagnostics detect raters who disagree with each other, not a panel that is wrong together. Only an external reference standard reaches the second failure; the card's job is to make the first kind visible instead of averaging it away. ### A divergent panel and the per-rater route The flag also catches asymmetry that cancels at the panel level. Three raters lean toward over-calling (sensitivity 0.95, specificity 0.75); two lean the other way (0.75, 0.95). Each coefficient returns roughly the same observed value, yet the qualities they imply split. ```r Y_div <- gen_logitnormal(6L, 5L, 1000L, Se = c(0.95, 0.75, 0.95, 0.75, 0.95), Sp = c(0.75, 0.95, 0.75, 0.95, 0.75), pi = 0.50) card_div <- grass_report(ratings = Y_div, bootstrap_B = 200) card_div #> GRASS Report Card #> #> sample = 5 raters, N = 1000, pi_hat = 0.50 #> PABAK = 0.49 -> 67th percentile <- primary #> AC1 = 0.49 -> 69th percentile #> Fleiss kappa = 0.49 -> 67th percentile #> ICC = 0.62 -> 65th percentile [distribution-sensitive] #> panel-agg. = suppressed (divergent) #> delta = 0.25 pp (divergent) #> matched null = (k=5, N=1000, q=0.85): delta_hat at the 99.5 percentile [design snapped] #> #> pairwise PABAK / surface percentile (lower / upper): #> R1 R2 R3 R4 R5 #> R1 -- 65% 71% 67% 74% #> R2 0.47 -- 62% 73% 64% #> R3 0.51 0.41 -- 62% 73% #> R4 0.48 0.52 0.42 -- 70% #> R5 0.55 0.46 0.54 0.50 -- #> #> per-rater vs panel-majority of OTHER raters (pooled-reference): #> R1 Se_tilde = 0.94 Sp_tilde = 0.77 (n_pos = 434, n_neg = 471, excl = 95) #> R2 Se_tilde = 0.71 Sp_tilde = 0.95 (n_pos = 466, n_neg = 439, excl = 95) #> R3 Se_tilde = 0.95 Sp_tilde = 0.71 (n_pos = 433, n_neg = 485, excl = 82) #> R4 Se_tilde = 0.74 Sp_tilde = 0.95 (n_pos = 468, n_neg = 441, excl = 91) #> R5 Se_tilde = 0.96 Sp_tilde = 0.78 (n_pos = 428, n_neg = 476, excl = 96) #> #> per-rater (latent-class fit; alongside pairwise): #> R1 Se = 0.95 (0.21, 0.97) Sp = 0.78 (0.05, 0.82) #> R2 Se = 0.73 (0.04, 0.78) Sp = 0.95 (0.28, 0.97) #> R3 Se = 0.96 (0.26, 0.98) Sp = 0.72 (0.04, 0.75) #> R4 Se = 0.75 (0.05, 0.79) Sp = 0.95 (0.25, 0.97) #> R5 Se = 0.96 (0.21, 0.98) Sp = 0.78 (0.04, 0.82) #> #> Notes: #> - delta_hat is the implied-quality spread over the agreement family (PABAK, mean AC1, Fleiss kappa). ICC is reported on the panel but does not enter delta_hat (v0.5.0 scope: ICC's reference depends on the full subject-prevalence distribution F and does not share the (q, pi_+) sufficient statistic the agreement family does). #> - flag from delta_hat's percentile on the matched null (k=5, N=1000, q=0.85; 50,000 draws); design snapped to nearest calibrated cell. #> #> See `summary(...)` for full panel and CI details. #> See `plot(...)` for a surface-position visualization. ``` Three things change in the printed output. The panel-aggregate summary is marked `suppressed (divergent)`. The pairwise PABAK matrix and the per-rater pooled-reference table appear, exposing the structure pair by pair and rater by rater. And the latent-class table appears alongside, pinning R1, R3, R5 near (0.95, 0.75) and R2, R4 near (0.75, 0.95); its bootstrap intervals are wide near balanced prevalence because the resamples visit both the correct mode and its label-switched mirror. Under the divergent flag the recommended primary deliverable is `pairwise_agreement()`: ```r pw <- pairwise_agreement(Y_div) pw #> GRASS Pairwise Reliability #> #> sample = 5 raters, N = 1000, pi_hat = 0.50, tau2_hat = 0.122, axis = inter #> #> Pairwise PABAK (lower triangle = PABAK_ij; upper triangle = surface percentile): #> #> R1 R2 R3 R4 R5 #> R1 -- 65% 71% 67% 74% #> R2 0.47 -- 62% 73% 64% #> R3 0.51 0.41 -- 62% 73% #> R4 0.48 0.52 0.42 -- 70% #> R5 0.55 0.46 0.54 0.50 -- #> #> Per-rater behavior against pooled panel-majority: #> #> rater Se_tilde Sp_tilde n_pos n_neg n_excl #> R1 0.94 0.77 434 471 95 #> R2 0.71 0.95 466 439 95 #> R3 0.95 0.71 433 485 82 #> R4 0.74 0.95 468 441 91 #> R5 0.96 0.78 428 476 96 #> (Se_tilde, Sp_tilde are calls vs panel-majority of OTHER raters, #> not against external truth. n_excl: subjects with tied majority, #> excluded from the per-rater pool.) ``` Read the per-rater table first. The reference is the observable panel majority rather than a latent class, so label-switching is moot: an over-caller shows high `Se_tilde` and lower `Sp_tilde`, an under-caller the mirror. In the matrix, the lower triangle holds the pairwise agreements and the upper triangle their percentiles at each pair's own marginal. The same PABAK value can land at different percentiles across pairs because each pair sees a different positive rate; the flag is exposing exactly the structure a panel average would have hidden. ### Two raters The headline call works at k = 2 with no signature change. ```r set.seed(42) N <- 200 truth <- rbinom(N, 1, 0.30) Y_k2 <- cbind( ifelse(truth == 1, rbinom(N, 1, 0.85), rbinom(N, 1, 0.15)), ifelse(truth == 1, rbinom(N, 1, 0.85), rbinom(N, 1, 0.15)) ) grass_report(ratings = Y_k2) #> GRASS Report Card #> #> sample = 2 raters, N = 200, pi_hat = 0.39 #> PABAK = 0.50 -> 69th percentile | quality 0.80-0.90 <- primary #> AC1 = 0.52 -> 69th percentile | quality 0.80-0.90 #> read: this panel agreed more tightly than 69% of what panels at this design can produce; the data are consistent with panel quality 0.80-0.90. #> delta = 0 pp implied-quality spread (not_applicable) #> matched null = n/a at k = 2 (coefficients cannot disagree; see pairwise/bounds path) #> #> Notes: #> - delta_hat is not applicable at k = 2: the two-coefficient agreement family (PABAK, AC1) implies identical panel quality by construction, so cross-coefficient discordance cannot be observed. Use the k = 2 identifiable bounds and pairwise path for asymmetry assessment. #> #> See `summary(...)` for full panel and CI details. #> See `plot(...)` for a surface-position visualization. ``` With two raters the agreement family reduces to PABAK and AC1, which imply identical quality by construction, so `delta_hat` reports `not_applicable` rather than a flag. Per-rater sensitivity and specificity are not point-identified from a two-rater matrix, so the latent-class fit returns Hui-Walter bounds instead of point estimates: ```r latent_class_fit(ratings = Y_k2, B = 200) #> grass latent-class fit #> method : hui_walter #> raters (k) : 2 #> bootstrap B : 200 #> per-rater #> R1 Se in [0.395, 1.000] Sp in [0.605, 1.000] (bounds) #> R2 Se in [0.385, 1.000] Sp in [0.615, 1.000] (bounds) #> #> Note: at k = 2, per-rater Se/Sp are not point-identified #> without external information. The intervals are #> Hui-Walter (1980) inequality bounds, not point #> estimates with sampling uncertainty. ``` The bounds carry the same meaning as the k >= 3 point estimates, an honest statement of what the data say about each rater. A study that needs per-rater points should add a third rater. ## Layered access Every field beneath the printed card is available on demand. The examples reuse the aligned case. ```r card1 <- grass_report(Y1, verbose = FALSE) summary(card1) #> GRASS Report Card -- summary #> #> sample : k = 5 raters, N = 1000, pi_hat = 0.237, axis = inter #> tau2_hat : 0.065 #> #> primary coefficient #> name : pabak #> observed : 0.537 #> percentile : 73.86 pp (pooled; position in the design's achievable range) #> q_hat : 0.866 #> band : consistent with panel quality 0.85-0.90 (95%) #> basis : pooled-achievable-range #> #> delta (cross-coefficient asymmetry) #> implied-q spread (pp) : 0.00 #> flag : aligned #> delta percentile : 42.2 pctile on matched null (k=5, N=1000, q=0.85) #> #> panel (full table) #> pabak observed = 0.537 q_hat = 0.866 pct = 73.86 pp quality 0.85-0.90 ref = closed-form #> mean_ac1 observed = 0.637 q_hat = 0.866 pct = 73.86 pp quality 0.85-0.90 ref = closed-form #> fleiss_kappa observed = 0.361 q_hat = 0.866 pct = 73.86 pp quality 0.85-0.90 ref = closed-form #> icc observed = 0.484 q_hat = 0.822 pct = 62.50 pp quality 0.80-0.85 ref = fitted-icc #> #> notes #> - Consistency band narrower than the calibrated q-grid spacing; endpoints interpolated within one grid gap. #> - Fitted-ICC F_key picked via glmer: mu_hat=-1.775, tau2_hat=3.089 -> F_key tau2=4.0000, mu=-1.386. #> - Fitted-ICC reference (GLMM-gap corrected) at F_key=LN_mu=-1.386_tau2=4.0000, k=5, N=1000 (family=logit_normal, M1=0.300). #> - delta_hat is the implied-quality spread over the agreement family (PABAK, mean AC1, Fleiss kappa). ICC is reported on the panel but does not enter delta_hat (v0.5.0 scope: ICC's reference depends on the full subject-prevalence distribution F and does not share the (q, pi_+) sufficient statistic the agreement family does). #> - flag from delta_hat's percentile on the matched null (k=5, N=1000, q=0.85; 50,000 draws); design snapped to nearest calibrated cell. #> #> grass version : 0.7.4 #> timestamp : 2026-07-12 08:27:47 ``` `summary()` is the audit view: the full panel with percentiles and bands, every caveat from the surface lookup, the package version, the timestamp. ```r head(as.data.frame(card1), 5) #> coefficient observed_value surface_percentile band_lo band_hi band_open_low #> 1 pabak 0.5372000 73.86364 0.85125 0.89875 FALSE #> 2 mean_ac1 0.6372062 73.86364 0.85125 0.89875 FALSE #> 3 fleiss_kappa 0.3609191 73.86364 0.85125 0.89875 FALSE #> 4 icc 0.4842582 62.50000 0.80125 0.84875 FALSE #> band_open_high q_hat se_q_hat clamped reference_used in_delta_hat #> 1 FALSE 0.8664693 0.005014575 FALSE closed-form TRUE #> 2 FALSE 0.8664471 0.005014821 FALSE closed-form TRUE #> 3 FALSE 0.8664693 0.005014575 FALSE closed-form TRUE #> 4 FALSE 0.8217253 0.002468171 FALSE fitted-icc FALSE #> is_primary delta_hat delta_percentile delta_flag matched_null_k #> 1 TRUE 0.00221665 42.21665 aligned 5 #> 2 FALSE 0.00221665 42.21665 aligned 5 #> 3 FALSE 0.00221665 42.21665 aligned 5 #> 4 FALSE 0.00221665 42.21665 aligned 5 #> matched_null_N matched_null_q #> 1 1000 0.85 #> 2 1000 0.85 #> 3 1000 0.85 #> 4 1000 0.85 ``` `as.data.frame()` is the tidy view, one row per coefficient, for binding many studies' cards into one table. ```r plot(card1) ```
Surface position plot

plot of chunk plot-surface-card

The default `plot()` view places the panel on the expected-value surface at its design, with the observation pinned. Other views (`type = "thermometer"` for the `delta_hat` gauge, `"panel"`, `"intervals"`, `"per_rater"`, `"diagnostic"`) are documented in `?plot.grass_card`. ## Building blocks Three of the internals `grass_report()` calls are exported. ```r position_on_surface(ratings = Y1, metric = "pabak") #> grass surface-position report (sweep convention, v0.7.1) #> metric : pabak #> observed value : 0.537 #> design (pi_hat,k,N) : (0.237, 5, 1000) #> implied quality q_hat: 0.866 +/- 0.005 #> pooled percentile : 73.9 (of the design's achievable range) #> consistency band : consistent with panel quality 0.85-0.90 (95%) #> sampling method : empirical #> notes : #> - Consistency band narrower than the calibrated q-grid spacing; endpoints interpolated within one grid gap. ``` `position_on_surface()` does the lookup for one coefficient and returns the percentile, the band, the implied quality, and the full sweep profile. ```r check_asymmetry(ratings = Y1) #> GRASS panel asymmetry diagnostic #> #> delta_hat = 0.0 pp (spread of the implied panel qualities) #> flag = aligned (42.2 percentile of matched null: k=5, N=1000, q=0.85) #> #> panel: #> coefficient observed implied q pooled pctile in delta_hat #> pabak 0.54 0.866 73.9 yes #> mean_ac1 0.64 0.866 73.9 yes #> fleiss_kappa 0.36 0.866 73.9 yes #> icc 0.48 0.822 62.5 no [distribution-sensitive] #> #> Surface caveats: #> - Consistency band narrower than the calibrated q-grid spacing; endpoints interpolated within one grid gap. #> - Fitted-ICC F_key picked via glmer: mu_hat=-1.775, tau2_hat=3.089 -> F_key tau2=4.0000, mu=-1.386. #> - Fitted-ICC reference (GLMM-gap corrected) at F_key=LN_mu=-1.386_tau2=4.0000, k=5, N=1000 (family=logit_normal, M1=0.300). #> - delta_hat is the implied-quality spread over the agreement family (PABAK, mean AC1, Fleiss kappa). ICC is reported on the panel but does not enter delta_hat (v0.5.0 scope: ICC's reference depends on the full subject-prevalence distribution F and does not share the (q, pi_+) sufficient statistic the agreement family does). #> - flag from delta_hat's percentile on the matched null (k=5, N=1000, q=0.85; 50,000 draws); design snapped to nearest calibrated cell. ``` `check_asymmetry()` returns the spread and its flag without per-rater detail. ```r latent_class_fit(ratings = Y_div, B = 200) #> grass latent-class fit #> method : dawid_skene_em #> raters (k) : 5 #> bootstrap B : 200 #> converged : TRUE #> iterations : 10 #> prevalence_hat: 0.475 #> log-likelihood: -2504.371 #> per-rater #> R1 Se = 0.950 (0.227, 0.968) Sp = 0.783 (0.048, 0.818) #> R2 Se = 0.729 (0.044, 0.767) Sp = 0.953 (0.267, 0.972) #> R3 Se = 0.956 (0.282, 0.977) Sp = 0.718 (0.047, 0.751) #> R4 Se = 0.750 (0.048, 0.785) Sp = 0.949 (0.252, 0.968) #> R5 Se = 0.962 (0.211, 0.981) Sp = 0.781 (0.035, 0.811) ``` `latent_class_fit()` produces the per-rater table directly. For prospective design, before any data exist, `plot_surface()` draws a coefficient's expected-value surface at a planned design and can pin a hypothetical observed value; see `?plot_surface`. ## Reference resolution and the clamp guard `delta_hat` inverts each coefficient to an implied panel quality. If an observed value falls outside the achievable range of its reference surface at the study's design, the inversion clamps to the boundary and flags `clamped = TRUE`, and `grass_report()` excludes clamped coefficients from `delta_hat` whenever at least two unclamped ones remain. Otherwise a single boundary-clamped coefficient could inflate the spread and fire a false divergent flag. Clamping is the exception, and the example below does not clamp. When the exact design is not a calibrated grid point, the reference resolves to the nearest calibrated cell and the resolution is disclosed in the `Notes`: ```r Y_edge <- gen_logitnormal(5L, 5L, 1000L, Se = rep(0.85, 5), Sp = rep(0.85, 5), pi = 0.08) card_edge <- grass_report(Y_edge, bootstrap_B = 50) card_edge$panel[, c("coefficient", "surface_percentile", "clamped", "reference_used")] #> coefficient surface_percentile clamped reference_used #> 1 pabak 67.75846 FALSE closed-form #> 2 mean_ac1 67.67121 FALSE closed-form #> 3 fleiss_kappa 67.78692 FALSE closed-form #> 4 icc 51.13636 FALSE fitted-icc ``` When a clamp does occur the card surfaces it twice, in the `Notes` and in the `panel` data frame's `clamped` and `reference_used` columns. The trade-off is deliberate: a coefficient with no calibrated reference at the user's design contributes nothing, and in exchange a fired flag reflects genuine disagreement rather than a boundary artifact. ## The machinery The sections above are everything a practitioner needs. This section records the definitions and calibration facts behind them; the accompanying paper carries the derivations and the full calibration appendices. **Notation.** `Y` is the N x k binary rating matrix. `pi_hat` is the observed marginal positive rate. `Se_j` and `Sp_j` are per-rater sensitivity and specificity inferred from the matrix, not against external truth. `q` is diagonal rater quality (`Se = Sp = q` for every rater under the calibration case). `q_hat` is the implied panel quality recovered from an observed coefficient. `delta_hat` is the spread of implied qualities, in quality percentage points, over PABAK, mean AC1, and Fleiss kappa. **Recovering `q_hat`.** PABAK has an algebraic inverse, `q_hat = (1 + sqrt(PABAK)) / 2`. Fleiss kappa, AC1, Krippendorff alpha, and ICC invert by a 501-point grid lookup with linear interpolation. Out-of-range observations clamp to the boundary with a flag (previous section). **ICC and the subject-prevalence distribution.** The agreement family's expectations depend only on the panel marginal and the rater quality. ICC depends additionally on the full distribution of subject-level positive probabilities: two studies can share a 30% positive rate and differ on whether it comes from obvious positives or borderline cases. The bundled reference covers a 48-point logit-normal grid plus four discrete-mixture profiles; a panel whose true distribution sits elsewhere carries a mismatch cost the agreement family does not. That is why ICC carries its `[distribution-sensitive]` marker and stays out of `delta_hat`. Variant selection (ICC(1,1), ICC(2,1), ICC(3,1)) follows the Koo and Li (2016) decision tree through the `metric` argument. **The matched-null flag convention.** The null distribution of `delta_hat` moves with panel quality and design, and at small N it includes a point mass at zero, so no fixed cut in quality points works everywhere. The package reads the observed spread's percentile on a bundled null at the matched `(k, N, q_hat)` cell, exactly how every coefficient on the card is already read. The nulls are generated through the production pipeline itself, pooled over prevalence: 440 cells at 50,000 draws each, with the 385 cells at k >= 3 shipped (the flag is `not_applicable` at k = 2). Monte Carlo error on a reported percentile stays under a quarter of a percentile point at the shipped draw counts, and cells with unstable extreme tails carry a transparency flag in the bundled object. The card prints the matched cell, its draw count, and any design snap. **Pairwise PABAK and pooled-reference identifiability.** PABAK is the pairwise coefficient because under symmetric raters its expectation reduces to `(2q - 1)^2`, a function of quality alone; every rater pair has its own observed marginal, and Cohen's kappa at the pair level is prevalence-attenuated (the kappa paradox of Byrt, Bishop, and Carlin, 1993). The pooled-reference estimates read each rater against the observable panel majority, so they are label-flip clean whenever the majority tracks the truth direction; under a panel whose minority is right, they read backwards. A 54-cell Monte Carlo validation recovers each rater's true bias direction in 85.8% of in-majority draws overall, with a per-cell median of 96.1%. Recovery exceeds 95% at the designs where the divergent flag is the intended deliverable (k of 5 or more at asymmetry 0.2 and above); small panels at mild asymmetry (k = 3, A = 0.1) are the structural floor. The package reports the per-rater pool size alongside each estimate. **Where the rest lives.** The accompanying paper carries the reference-surface construction, the closed-form derivations, the calibration grids (44,616 surface cells; 22-million-draw null and power programs), the flag's operating characteristics, and two published-panel reanalyses. The three case studies above are the paper's worked examples, regenerated from the same seeds; the quick-start and two-rater demonstration panels are this vignette's own. ## Extending the calibration Every number above is read from a simulated reference, and the reference is only as fine as the compute that built it. The bundled calibration holds 44,616 surface cells and 385 `delta_hat` null cells, built from 133 million simulated panels. Its bounds are compute limits, not method limits. Designs between grid points snap to the nearest calibrated cell. The ICC reference covers 52 subject-prevalence profiles. The `delta_hat` null pools five calibration prevalences, and 21 of its 385 cells flag unstable extreme tails. Denser grids, a wider profile compendium, a prevalence-stratified null, and deeper draws would each tighten the readings, and each needs replication rather than new theory. The package includes the tools to close these gaps. The bundled manifest lists every open calibration block. Each block is one seeded simulation cell, small enough to run on any machine with grassr installed: ```r manifest <- grass_calibration_manifest() table(manifest$program) #> #> lattice_k prev_strata tail_topup #> 220 1925 21 ``` `tail_topup` deepens the 21 flagged null cells, `prev_strata` rebuilds the `delta_hat` null with full draws in every prevalence stratum, and `lattice_k` adds rater counts between the shipped grid points. To see what a time budget buys on your machine, ask for a plan. Block selection is randomized so uncoordinated contributors rarely collide; we pin the selection here so the vignette reproduces: ```r set.seed(7) plan <- grass_contribute(dir = tempdir(), hours = 5, dry_run = TRUE) plan[, c("block_id", "program", "k", "N", "q", "prev", "draws")] #> block_id program k N q prev draws #> 1 33 prev_strata 3 15 0.85 0.20 50000 #> 2 149 prev_strata 3 100 0.65 0.50 50000 #> 3 161 prev_strata 3 100 0.85 0.95 50000 #> 4 197 prev_strata 3 200 0.65 0.05 50000 #> 5 232 prev_strata 3 300 0.85 0.05 50000 #> 6 244 prev_strata 3 300 0.97 0.50 50000 #> 7 253 prev_strata 3 500 0.75 0.20 50000 #> 8 327 prev_strata 5 20 0.75 0.05 50000 #> 9 335 prev_strata 5 20 0.85 0.80 50000 #> 10 337 prev_strata 5 20 0.92 0.05 50000 #> 11 342 prev_strata 5 20 0.97 0.05 50000 #> 12 440 prev_strata 5 100 0.92 0.80 50000 #> 13 464 prev_strata 5 150 0.92 0.50 50000 #> 14 489 prev_strata 5 200 0.92 0.50 50000 #> 15 550 prev_strata 5 1000 0.65 0.80 50000 #> 16 583 prev_strata 6 15 0.85 0.20 50000 #> 17 777 prev_strata 6 300 0.75 0.05 50000 #> 18 904 prev_strata 8 30 0.75 0.50 50000 #> 19 1024 prev_strata 8 200 0.65 0.50 50000 #> 20 1030 prev_strata 8 200 0.75 0.80 50000 #> 21 1068 prev_strata 8 300 0.97 0.20 50000 #> 22 1122 prev_strata 10 15 0.65 0.05 50000 #> 23 1151 prev_strata 10 20 0.65 0.95 50000 #> 24 1174 prev_strata 10 30 0.65 0.50 50000 #> 25 1215 prev_strata 10 50 0.92 0.80 50000 #> 26 1216 prev_strata 10 50 0.92 0.95 50000 #> 27 1427 prev_strata 15 20 0.75 0.05 50000 #> 28 1483 prev_strata 15 50 0.85 0.20 50000 #> 29 1500 prev_strata 15 75 0.65 0.80 50000 #> 30 1634 prev_strata 15 500 0.85 0.50 50000 #> 31 1706 prev_strata 25 20 0.75 0.95 50000 #> 32 1765 prev_strata 25 50 0.92 0.80 50000 #> 33 1869 prev_strata 25 200 0.97 0.50 50000 #> 34 1872 prev_strata 25 300 0.65 0.05 50000 #> 35 1928 prev_strata 25 1000 0.75 0.20 50000 #> 36 2062 lattice_k 12 20 0.65 NA 50000 ``` The same call without `dry_run = TRUE` runs the selected blocks from their declared seeds and writes a submission bundle to `dir`: one result file per block plus a manifest with checksums and session information. `grass_verify_contribution()` checks the bundle before you send it. Submission is a pull request against the `calibration-contrib` branch of the source repository (https://github.com/defense031/grassr). The maintainers verify a submission by re-executing its seeds, which reproduce bit for bit at the same package version, and merged cells fold into the bundled calibration in a future data release. The pipeline is already live. Outside contributors had run and verified 165 of the 2,166 blocks when this release was built. The repository manifest tracks which blocks remain open. Duplicate runs of a block are not wasted, they cross-verify. For a large budget, claim a block range through a repository issue first; the full protocol is in the CONTRIBUTING file. The bundled surfaces were themselves computed on two machines and merged by cell. The wider simulation programs (the reference surfaces and the ICC profile compendium) are in the repository under `simulation/` for contributors who want more than the packaged blocks. ## References - Semmel, A. & Gidaro, R. (2026). *A context-conditioned reporting convention for rater reliability on binary outcomes.* (Working paper.) The package is the computational companion to this work. - Byrt, T., Bishop, J. & Carlin, J. B. (1993). Bias, prevalence and kappa. *Journal of Clinical Epidemiology* 46(5), 423-429. - Cohen, J. (1960). A coefficient of agreement for nominal scales. *Educational and Psychological Measurement* 20(1), 37-46. - Dawid, A. P. & Skene, A. M. (1979). Maximum likelihood estimation of observer error-rates using the EM algorithm. *Applied Statistics* 28(1), 20-28. - Fleiss, J. L. (1971). Measuring nominal scale agreement among many raters. *Psychological Bulletin* 76(5), 378-382. - Gwet, K. L. (2008). Computing inter-rater reliability and its variance in the presence of high agreement. *British Journal of Mathematical and Statistical Psychology* 61(1), 29-48. - Hui, S. L. & Walter, S. D. (1980). Estimating the error rates of diagnostic tests. *Biometrics* 36(1), 167-171. - Koo, T. K. & Li, M. Y. (2016). A guideline of selecting and reporting intraclass correlation coefficients for reliability research. *Journal of Chiropractic Medicine* 15(2), 155-163.