grassr answers one question about a rating study: how good are these
raters? Give grass_report() a matrix of binary ratings, subjects in
rows and raters in columns, and it returns a Report Card with the
coefficient to report, where that value falls among everything a study
with this many raters, this many subjects, and this positive rate could
produce, and the range of rater quality the data are consistent with,
where quality is the share of subjects a rater labels correctly. If the
coefficients imply different qualities, the card flags the panel and
shows the raters one at a time. grass_power() sizes a rating study
before any data exist.
The package does not use fixed kappa benchmarks such as the Landis and Koch (1977) scale. Two raters who each label 87% of subjects correctly read as slight at 5% prevalence, moderate at 50%, and fair at 90%, even though the underlying rater quality never changed.
One matrix or data frame, subjects in rows and raters in columns, at
least two raters. Cells are numeric 0/1, logical, or a two-level
factor or character column, where 1, TRUE, or the positive level
means present. For a factor or character column the level named yes,
true, positive, present, or case, in any letter case, is taken as
positive; any other pair of labels is refused with a request for a 0/1
recode. Identifier columns such as a subject ID must be dropped before
the call. Column names, if given, become rater labels on the card.
Ratings must be complete (no missing data). The functions stop on NA
rather than dropping rows, so remove incomplete rows first with
Y <- Y[complete.cases(Y), ]; the card never sees the dropped rows.
Every row must be a different subject. A 0/1 matrix does not say which rows share a patient, and stacking repeated measures of one patient as rows overstates the number of subjects and narrows the quality band falsely. Build one matrix per repetition or session and compare the cards.
The example panel is simulated so the card can be checked against a known answer: 150 subjects at 30% prevalence, three raters who each label 87% of subjects correctly.
library(grassr) set.seed(7) truth <- rbinom(150, 1, 0.30) # the true class of each subject Y0 <- sapply(1:3, function(j) # three raters, each 87% correct 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
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 -> 71st percentile | quality 0.81-0.90 [distribution-sensitive] #> read: this panel's agreement exceeds 70% of what panels in this study #> context can produce; the data are consistent with panel quality 0.80-0.90. #> delta = 0.00 pp implied-quality spread (aligned) #> matched null = (k=3, N=150, q=0.86): delta_hat at percentile 11.6 of the #> null #> #> Notes: #> - ICC row: N=150 is not a calibrated ICC size, so its reference uses the #> nearest on the log scale, N=200. PABAK, AC1, Fleiss kappa, and #> delta_hat interpolate at N=150. #> #> See `summary(...)` for full panel and CI details. #> See `plot(...)` for a surface-position visualization.
Reading the card, line by line:
sample is the study context: rater count k, subject count N,
and the observed positive rate pi_hat, the share of all N x k
ratings that are positive. Here pi_hat is 0.36
against a true prevalence of 0.30, because false positives on the
many negatives outnumber misses on the few positives. Every reference
panel is simulated in this context.<- primary marks the coefficient to report: PABAK
(prevalence-adjusted bias-adjusted kappa; Byrt, Bishop, and Carlin,
1993) at moderate prevalence, and AC1 (Gwet, 2008) when the positive
rate is below 0.15 or above 0.85 (0.20 or 0.80 at two raters), where
its percentile is less noisy. The ICC is ICC(1,1) and is tagged
[distribution-sensitive] because its reference depends on how the
positives are spread across subjects, not only on how many there are,
so its percentile is less certain than the other three.read is the headline as a sentence, on the primary coefficient.
There is no verdict label; the width of the band is the uncertainty.delta is the spread of the qualities implied by PABAK, AC1, and
Fleiss' kappa, in percentage points (pp) of quality, followed by its
flag: aligned, caution, or divergent.In one sentence, this card says: "Three raters rated 150 subjects, observed positive rate 0.36. PABAK was 0.52, at the 70th percentile of what a three-rater, 150-subject study at this positive rate can produce, and the data are consistent with panel quality 0.80 to 0.90 (spread 0.00 pp, aligned)."
The band comes from a sweep across the calibrated qualities. At each quality the reference holds the panels that quality produces in this study context; the curve is the share of them with PABAK below the observed 0.52. The band is every quality where that share stays between 2.5% and 97.5%.
The consistency band is every quality at which the observed value is neither unusually low nor unusually high.
Two raters. The card is shorter at k = 2. Fleiss' kappa and the ICC
need three or more raters, Cohen's kappa has no calibrated reference,
and PABAK and AC1 always imply the same quality, so the flag is
not_applicable. Per-rater sensitivity and specificity are returned as
Hui and Walter (1980) bounds rather than point estimates.
Intra-rater. One rater scoring the same subjects on repeated
occasions uses the same layout with the occasions as columns and the
call grass_report(Y, axis = "intra"). ICC becomes the primary
coefficient; the three agreement coefficients still print and still
feed the flag.
The flag answers one question: do these raters differ from each other
enough that one number would hide it? aligned means no, and the card
reports the panel as one coefficient with its band. divergent means
yes, at least one rater is behaving differently from the rest, and the
card withholds the one-number summary and shows the raters
individually so the odd one out can be found. caution is the border.
The test behind it: each of PABAK, AC1, and Fleiss' kappa implies a
panel quality, and delta_hat is the gap between the highest and
lowest of the three. Equal-quality raters produce some gap by chance,
so the observed gap is compared with what chance produces at this
rater count, sample size, quality, and prevalence, the matched null the
card names. Below the null's 95th percentile the flag is aligned,
from the 95th caution, from the 99th divergent.
Here five raters score 1,000 subjects at true prevalence 0.30; three have sensitivity and specificity 0.92, and two keep that sensitivity but drop to specificity 0.51.
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_split <- 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) grass_report(Y_split, 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] #> summary = suppressed (divergent) #> delta = 0.42 pp (divergent) #> matched null = (k=5, N=1000, q=0.80): delta_hat at percentile 99.5 of the #> null #> #> 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 the majority of the OTHER raters: #> 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) #> (excl = subjects the other raters tied on, left out of that rater's pool) #> #> per-rater (latent-class fit): #> 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) #> #> See `summary(...)` for full panel and CI details. #> See `plot(...)` for a surface-position visualization.
The coefficients look ordinary, but the spread of implied quality lies
beyond the 99th percentile of what chance produces in this study (the
card prints 99.5, the top of its stored null), so the card drops the bands and the read: line, marks the summary
suppressed, and shows the raters one at a time in three views:
Se_tilde is the share of
subjects the others called positive that this rater also called
positive; Sp_tilde is the same share on the subjects the others
called negative. R4 and R5 are over-callers: Sp_tilde 0.51 and 0.50
against 0.91 to 0.93 for the rest.bootstrap_B resamples (1,000 by default; 200 here to
run quickly). Because nothing in the ratings says which hidden class
is the positive one, it can report every rater's two numbers
swapped. The majority table is the safe view and the latent-class
fit is the check.This panel does not need more raters or more subjects. It needs R4 and R5 retrained on when to call a subject negative. In one sentence: "Five raters rated 1,000 subjects, observed positive rate 0.47. The qualities implied by the three agreement coefficients were 0.42 percentage points apart, beyond the 99th percentile of what equal-quality raters produce in this study context, so the panel is not summarized by one coefficient. Against the majority of the other raters, R4 and R5 had specificity 0.51 and 0.50 versus 0.91 to 0.93 for the rest."
A caution card keeps its bands and is read the same way as an aligned
card, after a look at the pairwise table.
The flag detects raters who disagree with each other. A panel that is wrong together, every rater sharing the same specificity deficit, raises no flag, because no internal comparison can see it. Only an external reference standard catches that case.
grass_power()Assume a prevalence and a rough rater quality, and grass_power()
gives the probability that a study of a given size shows panel quality
above a value you choose, and the size at which that probability
reaches the level you want. Fix four of quality, prevalence, rater
count, sample size, and power, and it solves for the fifth. Give
prevalence = when planning; give pi_hat = only when a card is
already in hand, since the card reports the observed positive rate and
the function converts between the two. Here the panel is assumed to be
0.90 and a screening study with three raters at 10% prevalence is sized
to show it is above 0.80.
pw <- grass_power("fleiss_kappa", q = 0.90, q0 = 0.80, prevalence = 0.10, k = 3, power = 0.80) pw #> #> GRASS power analysis: show panel quality above 0.80 (panel assumed 0.90) #> coefficient: Fleiss kappa #> #> q = 0.90 #> prevalence = 0.10 #> pi_hat = 0.18 (implied) #> k = 3 #> N = 61 <- solved #> power = 0.80 #> #> Power is the probability that a study of this size shows panel quality #> above 0.80 when the panel is 0.90. #> See `plot()` for the curve over Number of subjects (N). plot(pw)
Power curve over sample size for the planned screening study, with the requested power as a dashed line and the solved N marked.
Solve for k instead and the function returns the smallest rater
count that reaches the power, or NA with the reason. Rater count
snaps to the calibrated grid (2, 3, 5, 8, 15, 25). Leave the prevalence
empty and it returns the range of prevalences over which a fixed design
keeps that power. target = sizes a study against a fixed coefficient
value instead, for a threshold imposed from outside. Plan on the
coefficient the card will mark primary.
grass_report() composes four functions, each callable on a rating
matrix:
position_on_surface() returns the percentile, band, and implied
quality for one coefficient. Given a hypothetical value and a study
context instead of a matrix, it answers a planning question.check_asymmetry() returns delta_hat and its flag without
per-rater detail.pairwise_agreement() returns the pairwise matrix and the majority
table for any panel, flagged or not.latent_class_fit() returns the per-rater sensitivity and
specificity table, and at two raters the Hui and Walter bounds.plot_surface() draws a coefficient's expected value over prevalence
and quality for a planned rater count and sample size, with an
observed value pinned on it.
A card has three other views. summary(card) prints the full panel,
every caveat from the surface lookup, the package version, and the
timestamp. as.data.frame(card) returns one row per coefficient, for
binding many studies' cards into one table; delta_hat is in quality
percentage points, the unit the card prints. plot(card) places the
panel on the primary coefficient's expected-value surface over
prevalence and quality.
Other plot() types (type = "thermometer" for the delta_hat gauge,
"panel", "intervals", "per_rater", "diagnostic") are documented
in ?plot.grass_card.
Every percentile and band comes from a bundled simulated reference.
Study contexts between calibrated cells interpolate in quality, sample
size, and prevalence; rater count snaps to the nearest calibrated
value, and the card says so. A coefficient outside the range its study
context can produce is clamped to the nearest edge, marked
clamped = TRUE in the panel data frame, and dropped from delta_hat.
The simulation programs are in the source repository under
simulation/.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.