Analysing rater and response data

Many studies produce the same shape of data: several people (raters, annotators, respondents, workers) each answer some of a set of items, and the analyst wants to know how much they agree, what the “true” answer to each item is, how reliable each rater is, whether two groups differ, and which of many comparisons survive correction. symplex::stats answers those questions exactly: every coefficient that is a rational function of the counts is an exact rational, every test statistic is an exact expression, and every p-value is an exact expression (or an exact rational for the discrete exact tests) that is only rounded when you ask with .eval_f64(). The oracles the tests cite are statsmodels, scipy, the krippendorff package and Python’s statistics module on Fractions.

The walk-through below is tests/v13/v13_walkthrough.rs; every number in it is produced by the library.

The data

Eight items rated by five raters into three categories 0, 1, 2, with two missing ratings:

use symplex::stats::agreement::*;
let table = RatingTable::from_i64_missing(&[
    &[Some(0), Some(0), Some(0), Some(0), Some(1)],
    &[Some(1), Some(1), Some(1), Some(2), Some(1)],
    &[Some(2), Some(2), Some(2), Some(2), Some(2)],
    &[Some(0), Some(1), Some(0), None,    Some(0)],
    &[Some(1), Some(1), Some(2), Some(1), Some(1)],
    &[Some(2), Some(1), Some(2), Some(2), None   ],
    &[Some(0), Some(0), Some(1), Some(0), Some(0)],
    &[Some(1), Some(2), Some(1), Some(1), Some(1)],
])?;

RatingTable is items × raters with Option<Q> cells (Q is the exact rational Ratio<BigInt>; symplex::linprog::{q, qi} build literals). from_i64, from_rows, and the raters-first from_raters_i64 (the krippendorff package’s layout) are the other constructors.

How much do the raters agree?

Every one of these is an exact rational: the krippendorff package reports 0.4524312896405921 for the nominal α, which is 214/473. The coefficients that need complete tables (Fleiss, ICC, Kendall’s W) say so in their errors; Krippendorff’s α is the one to reach for when cells are missing. confusion_matrix, category_frequencies and RatingTable::count_table expose the counts the coefficients are built from.

What is the answer to each item?

use symplex::stats::aggregation::*;
let labels = LabelTable::from_rows(&rows, 3)?;      // items × raters, Option<usize>
let votes = majority_votes(&labels);                // Vote { winner, tied, counts }
let ds = dawid_skene(&labels, &DawidSkeneOpts::default())?;
ds.labels();          // one label per item (argmax posterior)
ds.confusion[j];      // rater j's estimated confusion matrix
ds.priors;            // class prevalences

majority_vote / majority_votes / plurality(labels, &threshold) / weighted_vote(labels, &weights) are exact. dawid_skene is the Dawid–Skene (1979) EM algorithm in f64 — deterministic, with max_iter, tol, smoothing and an initialisation choice — returning per-item posteriors, per-rater confusion matrices, class priors and convergence information; dawid_skene_counts takes replicate ratings. For pairwise judgements (“is A better than B?”) bradley_terry(&wins, &opts) fits the Bradley–Terry model by Hunter’s MM algorithm.

Priors on the raters: MAP Dawid–Skene and MACE

With few items per rater the maximum-likelihood confusion matrices degenerate — a rater who only ever answered 0 gets the rows (1, 0, 0) and then carries no information at all. Two Bayesian variants keep the estimates in the interior:

// Dirichlet priors: α on the class prevalences, β on every confusion row
// ([true][observed], shared by all raters).  The M-step is the posterior
// mode (count + α − 1) / Σ(count + α − 1), so every α ≥ 1; all ones is the MLE.
let priors = DawidSkenePriors::symmetric(3, 1.0, 3.0, 1.5);   // K, class α, diagonal β, off-diagonal β
let ds = dawid_skene_map(&labels, &priors, &DawidSkeneOpts::default())?;

// MACE (Hovy et al. 2013): rater r copies the true label with probability θ_r,
// otherwise "spams" a label from its own distribution ξ_r.
let m = mace(&labels, &MaceOpts::default())?;       // smoothing 0.1, majority-vote start
m.competence;         // θ_r per rater — who to trust
m.spam_distribution;  // ξ_r — what a spammer types
m.labels();           // argmax posterior per item

// Rank items by how undecided the raters left them, and check the models
// against gold labels where you have them.
posterior_entropy(&m.posteriors);                   // bits per item
rater_confusion_from_gold(&labels, &gold)?;         // exact [gold][given] per rater

dawid_skene_map with symmetric(K, 1, 1 + s, 1 + s) is exactly dawid_skene with smoothing = s; the priors matter when a rater’s rows would otherwise be estimated from two or three items. MACE has one parameter per rater instead of a K × K matrix, so it is the model to reach for when raters are many and their labels few — it finds the constant-answer spammer that a majority vote is fooled by. Both are deterministic EM iterations in f64 (MaceInit::Posteriors gives a start of your choosing in place of random restarts).

How good is each rater?

use symplex::stats::estimation::{proportion_interval, IntervalMethod};
let acc = worker_accuracy(&rater_labels, &gold)?;   // Accuracy { correct: 5, answered: 8, accuracy: Some(5/8) }
category_metrics(&rater_labels, &gold, 3)?;         // precision / recall / F₁ per category, exact
let ci = proportion_interval(5, 8, 0.95, IntervalMethod::ClopperPearson)?;   // Interval<f64>
(ci.lower, ci.upper);                                // (0.2449, 0.9148)
proportion_interval(5, 8, 0.95, IntervalMethod::Wilson)?;
gold_screening(&labels, &gold, &q(2, 3))?;          // pass / fail per rater on the gold items

The Clopper–Pearson interval is the exact one (statsmodels proportion_confint(method='beta')); Wilson, Agresti–Coull and Wald are the usual approximations. The proportion intervals are interval estimates, so since 0.18 they live in stats::estimation beside the mean intervals (the stats::aggregation paths still re-export them for one release). A rater’s accuracy against chance is an exact binomial test (below).

Exact intervals

The f64 function rounds; the same intervals exist with exact endpoints. Wald, Wilson and Agresti–Coull are closed algebraic forms in the normal quantile z, so they take any expression for z — a symbol for the textbook formula, or z_for_confidence(&ctx, &q(95, 100)) for the exact √2·erfinv(19/20) (1.959963984540054). Clopper–Pearson’s endpoints are the roots in (0, 1) of the two binomial-tail polynomials Σ_{j≥k} C(n,j) pʲ(1−p)ⁿ⁻ʲ − α/2 and Σ_{j≤k} … − α/2, which proportion_interval_exact returns as RootOf algebraic numbers:

let z = ctx.symbol("z");
let ci = proportion_interval_symbolic(&ctx, 5, 8, &z, IntervalMethod::Wilson)?;   // Interval<Ex> in z, contains z^2

let ci = proportion_interval_exact(&ctx, 5, 8, &q(95, 100), IntervalMethod::ClopperPearson)?;
ci.lower;                          // RootOf(…degree-8 polynomial in _p…, k)
ci.lower.eval_f64()?;              // 0.2448632163665516   = scipy beta.ppf(0.025, 5, 4)
ci.upper.eval_f64()?;              // 0.9147665858627464   = scipy beta.ppf(0.975, 6, 3)
ci.lower.eval_decimal(30)?;        // as many digits as you like

proportion_interval_exact(&ctx, 1, 1, &q(9, 10), IntervalMethod::ClopperPearson)?;   // [1/20, 1]: linear tail, rational root

The exact closed forms are not clipped to [0, 1] (Agresti–Coull at k = 0 has a negative lower end that the f64 function clamps), and the Clopper–Pearson roots need a Sturm isolation and a factorisation of a degree-n polynomial — fine for tens of trials, not thousands. The known-σ mean interval has the same pair, confidence_interval_mean_z_symbolic / _exact, in the same module.

Do two groups differ?

Response times of two groups, in seconds:

use symplex::stats::{data, hypothesis::*};
let fast = data::from_i64(&[12, 15, 11, 14, 13, 16, 10, 17]);
let slow = data::from_i64(&[18, 22, 19, 25, 20, 21, 23, 24]);
data::mean(&fast)?;                                          // 27/2
data::variance(&fast, data::Ddof::Sample)?;                  // 6
data::quantile(&slow, &q(3, 4), data::QuantileMethod::Inclusive)?;   // 93/4

let t = t_test_two_sample(&ctx, &fast, &slow, false, Alternative::TwoSided)?;   // Welch
t.statistic;          // exact: −(…)·√(…)   → −6.531972647421809
t.p_value;            // exact expression through betainc_regularized → 1.3298737271301488e-05
cohens_d(&ctx, &fast, &slow, true)?;

let u = mann_whitney_u(&ctx, &fast, &slow, Alternative::TwoSided, RankMethod::Exact)?;
u.statistic;          // 0
u.p_value;            // 1/6435 — the exact null distribution of U, as a rational

TestResult { statistic, p_value, df, alternative } carries exact expressions; p_value_f64() / statistic_f64() round them, and p_value_exact() returns the rational when there is one (the discrete exact tests). The family:

Approval counts by group, [[30, 10], [18, 22]]: fisher_exact gives p = 0.01150621201656047 exactly as a rational, and chi_square_independence(…, true) p = 0.01205961617749023 as a χ²-tail expression — both matching scipy to the last digit. Fisher’s p-value is an exact rational while the hypergeometric support has at most 2 000 points; for larger tables (cells in the 10⁴–10¹¹ range) it is numeric — the pmf walked from its mode in floating point, returned as ctx.from_f64(p) and agreeing with scipy to 1e-9 — and a [[10⁶, 10⁶ + 7], [10⁶ − 3, 10⁶]] table takes milliseconds. Below 1e-308 the expression is exp(ln p), so p_value_log10 still reads the tail.

Screening many raters at once

A rater who answered 14 of 20 gold questions correctly, against a chance rate of ⅓: binomial_test(&ctx, 14, 20, &q(1, 3), Alternative::Greater) has the exact p-value 1021403/1162261467 (≈ 8.79·10⁻⁴). Testing many raters multiplies the false positives; bonferroni, holm, benjamini_hochberg and benjamini_yekutieli return adjusted p-values and reject flags (matching statsmodels.multipletests):

let adj = benjamini_hochberg(&[0.001, 0.02, 0.03, 0.2, 0.8], 0.05)?;
adj.reject;   // [true, true, true, false, false]

bootstrap_ci and permutation_test (seeded through stats::Rng, so reproducible) cover statistics without a known null distribution; sample_size_for_proportion, sample_size_two_proportions, power_t_test_two_sample and sample_size_t_test_two_sample plan a study (the t-power uses the noncentral t by quadrature and matches statsmodels).

Estimating and modelling

stats::estimation fits distributions to data — fit_normal, fit_exponential, fit_poisson, fit_bernoulli, fit_geometric, fit_uniform, fit_log_normal by maximum likelihood (closed forms, exact where the estimate is rational), fit_gamma_moments, fit_beta_moments, fit_negative_binomial_moments by moments — and returns Distributions, so everything from the previous chapter (probabilities, quantiles, sampling) applies to the fitted model. log_likelihood, aic, bic compare fits. Conjugate Bayesian updating returns distributions too: beta_binomial_posterior(&ctx, &α, &β, successes, failures) is the posterior for a rater’s accuracy after s right and f wrong answers, credible_interval(&post, 0.95) its equal-tailed Interval<f64> (fields lower, upper), and posterior_predictive_beta_binomial the exact predictive table for the next n answers.

Sequences of states — a respondent moving between “engaged”, “guessing” and “gone”, say — are stats::markov::MarkovChain on an exact QMatrix: stationary_distribution() is one exact linear solve ([2/7, 3/7, 2/7] for the textbook 3-state chain), absorption_probabilities, expected_steps_to_absorption, expected_hitting_time, communication_classes, periods, and sample_path.

Descriptive statistics, exactly

stats::data has the everyday summaries on exact observations: mean, variance(Ddof::{Population, Sample}), std, median, quantile / quantiles (QuantileMethod::{Exclusive, Inclusive} — Python’s statistics.quantiles default and method='inclusive'), iqr, modes, frequencies, ranks (average ranks, as scipy.stats.rankdata), covariance, pearson, spearman, kendall_tau, skewness, kurtosis, median_abs_deviation, zscores, geometric_mean, harmonic_mean, trimmed_mean, and the outlier screens iqr_outliers (Tukey’s fences) and mad_outliers (modified z-scores) for response times. The ordinal association measures live here too: goodman_kruskal_gamma, somers_d(x, y, Dependent::{Y, X, Symmetric}), kendall_tau_c and the concordance_counts they are built from. from_f64 converts floats exactly (every f64 is a dyadic rational), to_f64 rounds back.

Is the questionnaire itself reliable?

When each item is meant to measure the same thing, stats::reliability asks whether they do — exactly. With respondents as rows and items as columns of a RatingTable:

Every coefficient that is a rational function of the scores is an exact rational (Cronbach’s α, KR-20); the ones with roots are exact expressions. Since 0.18 stats::reliability holds only scale reliability and item analysis; its former neighbours moved to the module their rule names (the old paths re-export them for one release):

Explaining accuracy or time by features

stats::regression fits models to the data you have about respondents:

use symplex::stats::regression::{ols, logit, LogitOpts, Design};
// Response time explained by two features, with an intercept — exactly.
let fit = ols(&y, &rows, true)?;                     // or Design::new().intercept().column(&x1).column(&x2).fit(&y)?
fit.coefficients;                                     // Vec<Q>, exact (XᵀX)⁻¹Xᵀy
fit.r_squared; fit.adjusted_r_squared;                // exact rationals
fit.standard_errors(&ctx)?;                           // exact expressions (√ of σ̂²(XᵀX)⁻¹)
fit.coefficient_tests(&ctx)?;                         // TestResult per coefficient, p through StudentT
fit.f_test(&ctx)?; fit.anova_table();                 // overall F, exact
fit.conf_int(0.95)?; fit.prediction_interval(&x_new, 0.95)?;   // Vec<Interval<f64>>, Interval<f64> — no ctx: the limits are f64
fit.leverage(); fit.cooks_distance(); fit.durbin_watson(); vif(&rows)?;
// Correct / incorrect explained by features: logistic regression (IRLS, f64).
let lg = logit(&correct, &features, true, &LogitOpts::default())?;
lg.coefficients; lg.odds_ratios(); lg.p_values; lg.pseudo_r_squared; lg.predict_proba(&x_new);

ols, wls, simple_linear_regression and polyfit are exact and match statsmodels OLS attribute for attribute; logit matches Logit to 1e-6 and reports perfect separation as an error rather than as enormous coefficients.

Several answers, or ordered answers

When the response has more than two categories, mnlogit fits a multinomial logit (statsmodels MNLogit: category 0 is the reference, one equation per other category) and ologit a proportional-odds model (statsmodels OrderedModel(distr='logit'): P(y ≤ j | x) = σ(θ_j − xβ), thresholds instead of an intercept). Both are Newton–Raphson in f64 with the same LogitOpts, and both refuse separated data with an error that names the diverging coefficient.

use symplex::stats::regression::{mnlogit, ologit, LogitOpts};
// Which of three answers a respondent picks, explained by one feature x = 1..24.
let y = [0, 0, 1, 0, 0, 1, 2, 0, 1, 1, 2, 0, 1, 2, 1, 2, 1, 2, 2, 1, 2, 0, 2, 2];
let x: Vec<Vec<f64>> = (1..=24).map(|i| vec![f64::from(i)]).collect();
let mn = mnlogit(&y, &x, true, &LogitOpts::default())?;
println!("{:.4} {:.4}", mn.coefficients[1][0], mn.coefficients[1][1]);   // -2.8750 0.2536  (answer 2 vs 0: intercept, slope)
println!("{:.4}", mn.relative_risk_ratios()[1][1]);                     // 1.2887  (× per unit of x for P(2)/P(0))
let p = mn.predict_proba(&[12.0])?;
println!("{:.3} {:.3} {:.3}", p[0], p[1], p[2]);                         // 0.272 0.406 0.322
println!("{}", mn.predict(&[12.0])?);                                    // 1  (the modal answer)
println!("{:.4}", mn.pseudo_r_squared);                                  // 0.1572

// A rating on an ordered scale 0 < 1 < 2, explained by x = 0..19: no intercept column.
let y = [0, 0, 1, 0, 1, 1, 2, 1, 2, 2, 0, 1, 2, 2, 1, 0, 0, 1, 2, 2];
let x: Vec<Vec<f64>> = (0..20).map(|i| vec![f64::from(i)]).collect();
let ol = ologit(&y, &x, &LogitOpts::default())?;
println!("{:.4}", ol.coefficients[0]);                                   // 0.1140  (β)
println!("{:.4} {:.4}", ol.thresholds[0], ol.thresholds[1]);             // 0.1212 1.7317  (θ₀ < θ₁, as thresholds)
println!("{:.4}", ol.odds_ratios()[0]);                                  // 1.1208  (odds of a higher rating, per unit of x)
let p = ol.predict_proba(&[7.0])?;
println!("{:.3} {:.3} {:.3}", p[0], p[1], p[2]);                         // 0.337 0.381 0.282
println!("{:.3}", ol.llr_test(&ctx)?.p_value_f64()?);                    // 0.126  (H₀: β = 0, χ²₁)

statsmodels reports the ordinal thresholds as θ₀ followed by the log-differences ln(θ_j − θ_{j−1}); ologit reports the thresholds themselves (statsmodels’ transform_threshold_params), and its standard errors are for those. Every printed number above is asserted in tests/v19/v19_mnlogit.rs.

Screening while the answers arrive

stats::sequential::Sprt::bernoulli(&p0, &p1, α, β) is Wald’s sequential probability ratio test: feed each gold-question outcome to update, and the moment the exact log-likelihood ratio (log_likelihood_ratio(&ctx)) crosses a boundary the decision is AcceptH0 (the rater performs at the chance rate p0) or AcceptH1 (at the competent rate p1); until then Continue. The decision is sticky — once reached, every later update returns it, is_decided() is true and stopped_at() gives the sample size at which the test stopped; reset() starts over. expected_sample_size_bernoulli says how many questions that takes on average.

Comparing label distributions

stats::information measures, exactly, how two raters’ (or two populations’) label distributions differ: kl_divergence, js_divergence, total_variation, hellinger, bhattacharyya_distance, cross_entropy; and from a joint table how much one variable tells about another: mutual_information, conditional_entropy, normalized_mutual_information (Norm::{Arithmetic, Geometric, Min, Max}), joint_from_counts. Logs are kept as exact ln expressions over prime factors, so H(½, ¼, ¼) is exactly 3/2 bits.

Time to completion, time to attrition

stats::survival handles right-censored durations — how long until a respondent finishes (or is still working when the study ends), how long a worker stays active:

use symplex::stats::survival::{KaplanMeier, Observation, CiMethod, log_rank_test};
let obs = Observation::from_i64(&[3, 5, 6, 7, 8, 10, 12, 12], &[true, false, true, true, false, true, true, false]);
let km = KaplanMeier::fit(&obs)?;
km.survival_at(&q(6, 1));               // 35/48 — an exact step function
km.variance_at(&q(6, 1));               // Greenwood, 1505/55296
km.cumulative_hazard_at(&q(6, 1));      // Nelson–Aalen, 7/24
km.median();                            // Some(10)
km.confidence_interval(&q(7, 1), 0.95, CiMethod::LogLog)?;   // Interval<f64>
km.restricted_mean(&q(12, 1));          // 1279/144
let r = log_rank_test(&ctx, &all_obs, &groups)?;   // exact χ² statistic 149059681/48496587 ≈ 3.0736, p ≈ 0.0796

The Kaplan–Meier steps, Greenwood variances and Nelson–Aalen hazards are exact rationals matching statsmodels’ SurvfuncRight; the log-rank statistic is an exact rational matching scipy’s logrank and statsmodels’ survdiff. exponential_rate is the censored MLE of a constant hazard, and survival_function / hazard_function give S(t) and h(t) of any Distribution as expressions.

Several measures at once, and extremes

stats::multivariate::MultivariateNormal (symbolic mean vector and covariance Matrix) has an exact density, marginals, conditionals by the Schur complement, Mahalanobis distance and sampling; covariance_matrix and correlation_matrix summarise several columns of ratings at once (exact), and pca finds their principal components exactly (pca_f64 for larger matrices). stats::order::{order_statistic, minimum_of, maximum_of} give the distribution of the k-th smallest of n independent draws as a Distribution in its own right — the fastest of n responses, the worst of n ratings — exactly for both continuous families and finite tables.