Adaptive comparative judgment: estimation, reliability, and evaluation

Comparative judgment (CJ) uses decisions about pairs of items to estimate their relative strengths. An adaptive run chooses later pairs using earlier results. Good pair selection, precise estimates, and accurate predictions are related, but each needs its own assessment:

Pair selection → parameter estimation → internal reliability → empirical score recovery → held-out predictive performance.

Step	Question	What to inspect
Pair selection	Which comparison should be collected next?	Strategy, allowed pairs, presentation order, and selection logs.
Parameter estimation	What strengths explain the collected outcomes?	Estimator, likelihood, adjustments, estimates, uncertainty, and convergence.
Internal reliability	How large is estimated item spread relative to estimated uncertainty?	SSR or Bayesian EAP reliability, with their definitions and provenance.
Empirical score recovery	How closely do scores reproduce an external reference?	Correlation, ranks, spread, and score differences against truth or a dense estimate.
Held-out predictive performance	How well do probabilities predict new judgments?	Log loss, Brier score, and pairwise accuracy on outcomes withheld from selection and fitting.

The arrows organize the workflow; they do not imply that passing one step proves the next. A dense score estimate is itself uncertain, and high internal reliability does not establish agreement with it, judge validity, or predictive calibration. Calibration means that predicted chances agree with observed win frequencies over comparable predictions.

This guide uses synthetic outcomes only. It makes no provider or network calls. The executable examples require the suggested withr package for local seeds; if it is absent, the code is displayed without execution. The conventional and Firth examples have additional optional-engine guards. Bayesian sampling is shown separately and never runs during document builds. Example sizes, penalties, and replicate counts illustrate interfaces, not recommended study settings.

Start with comparisons chosen in advance

Here theta is relative strength on a log-odds scale. For the simple BT model, the probability that A wins is plogis(theta[A] - theta[B]). Zero is a reference location, not a passing score. Separately fitted panels do not automatically share an origin or scale.

We choose all ordered pairs before drawing any winners. A fixed training set, a denser reference set, and an evaluation set use separate random streams. They share items and pair identities but contain independently generated judgments. No evaluation outcome is used to choose pairs, fit scores, or select settings. This assesses new judgments on known items, not generalization to unseen items or judges.

ids <- letters[1:6]
truth <- setNames(seq(-1, 1, length.out = length(ids)), ids)
pairs <- expand.grid(object1 = ids, object2 = ids, stringsAsFactors = FALSE)
pairs <- pairs[pairs$object1 != pairs$object2, ]
simulate_rows <- function(repeats, seed) {
  dat <- pairs[rep(seq_len(nrow(pairs)), each = repeats), ]
  p <- plogis(truth[dat$object1] - truth[dat$object2])
  dat$result <- withr::with_seed(seed, rbinom(nrow(dat), 1, p))
  rownames(dat) <- NULL
  dat
}
fixed_data <- simulate_rows(2L, 30601L)
reference_data <- simulate_rows(40L, 30602L)
heldout_data <- simulate_rows(20L, 30603L)
example_design <- data.frame(
  dataset = c("fixed training", "dense reference", "held-out evaluation"),
  outcome_seed = c(30601L, 30602L, 30603L),
  comparisons = c(nrow(fixed_data), nrow(reference_data), nrow(heldout_data)),
  selection = "all ordered pairs chosen before outcomes"
)
example_design
#>               dataset outcome_seed comparisons
#> 1      fixed training        30601          60
#> 2     dense reference        30602        1200
#> 3 held-out evaluation        30603         600
#>                                  selection
#> 1 all ordered pairs chosen before outcomes
#> 2 all ordered pairs chosen before outcomes
#> 3 all ordered pairs chosen before outcomes

fit_bt_model() accepts these three columns: object1, object2, and binary result, where 1 means the first item wins. For recorded winners in ID1/ID2/better_id or A_id/B_id/better_id columns, first use build_bt_data(). The getting-started guide demonstrates that conversion with bundled public data.

Conventional fitting and estimator provenance

The default engine = "auto" tries sirt, then BradleyTerry2 if sirt is missing or its execution fails. Explicit engines never switch to another estimator. Validation failures such as disconnected comparisons are not repaired by that automatic fallback. Use an explicit engine when an analysis requires a particular method, and retain the returned provenance, including effective settings, identification, convergence, engine version, and any fallback reason.

conventional_fit <- fit_bt_model(
  fixed_data, engine = "sirt", sirt_eps = 0.3, verbose = FALSE
)
summarize_bt_fit(conventional_fit)
#> # A tibble: 6 × 6
#>   ID      theta    se  rank engine reliability
#>   <chr>   <dbl> <dbl> <int> <chr>        <dbl>
#> 1 f      1.20   0.567     1 sirt         0.465
#> 2 e      0.360  0.491     2 sirt         0.465
#> 3 d     -0.0144 0.482     3 sirt         0.465
#> 4 b     -0.387  0.488     4 sirt         0.465
#> 5 a     -0.578  0.496     6 sirt         0.465
#> 6 c     -0.578  0.496     5 sirt         0.465
conventional_fit$ssr
#> $observed_variance
#> [1] 0.4758405
#> 
#> $mean_squared_se
#> [1] 0.2543473
#> 
#> $true_score_variance
#> [1] 0.2214933
#> 
#> $ssr
#> [1] 0.4654779
#> 
#> $n_items
#> [1] 6
#> 
#> $n_finite
#> [1] 6
#> 
#> $valid
#> [1] TRUE
#> 
#> $status
#> [1] "ok"
#> 
#> $engine_reliability
#> [1] 0.4654779
#> 
#> $agrees
#> [1] TRUE
#> 
#> $absolute_difference
#> [1] 0
#> 
#> $tolerance
#> [1] 1e-12
conventional_fit$provenance[c(
  "requested_engine", "engine", "engine_version", "package_version",
  "adjustment", "convergence", "fallback_reason"
)]
#> $requested_engine
#> [1] "sirt"
#> 
#> $engine
#> [1] "sirt"
#> 
#> $engine_version
#> [1] "4.2.133"
#> 
#> $package_version
#> [1] "1.6.0"
#> 
#> $adjustment
#> $adjustment$method
#> [1] "epsilon"
#> 
#> $adjustment$eps
#> [1] 0.3
#> 
#> 
#> $convergence
#> $convergence$status
#> [1] "iteration_limit_reached"
#> 
#> $convergence$converged
#> [1] NA
#> 
#> $convergence$iterations
#> [1] 100
#> 
#> 
#> $fallback_reason
#> NULL

The explicit sirt_eps above controls sirt’s conventional epsilon adjustment; omitting it preserves the installed engine’s default. It is not a lapse probability. sirt reports iterations rather than an explicit convergence flag; inspect its recorded stopping status as well as reliability.

What the adjustments change

Estimator	Point-estimate adjustment	Reported uncertainty
`engine = "sirt"`	Conventional epsilon adjustment depends on observed win proportions.	Engine SEs and independently checked SSR; retain resolved sirt settings.
`engine = "brglm2"`	Firth / mean bias reduction adds half the log determinant of expected information to the log likelihood (`AS_mean`).	Inverse original-data expected information at the adjusted estimate, transformed to centered item coordinates.
`engine = "alpha"`	Hamilton-style adjustment adds equal pseudo-wins in both directions for every possible pair, including unobserved pairs.	Inverse original-data expected information at the adjusted estimate, transformed to centered item coordinates.

Firth is an explicit option for random or nonadaptive schedules, including connected data with undefeated or winless items. Alpha adjustment is motivated by adaptive scheduling. The package does not infer a schedule type from outcomes, choose alpha, or establish that either adjustment is universally superior. Neither adjustment selects pairs.

For $N$ items, alpha adds $c=\alpha/(N-1)$ pseudo-wins in each direction to every unordered pair. With $p_{ij}=\operatorname{logit}^{-1}(\theta_i-\theta_j)$ , the objective is

$\ell_\alpha(\theta)=\ell(\theta)+ \frac{\alpha}{N-1}\sum_{i<j}\log\{p_{ij}(1-p_{ij})\}.$

The Firth and alpha covariance matrices use the original observed comparisons, not the inverse penalized Hessian. Their SEs are conditional on the realized comparison graph: adjusting point estimates does not automatically account for the variability of an outcome-dependent schedule. See fit_bt_model() for the exact objectives, numerical checks, and the Hamilton and Tawn reference.

firth_fit <- fit_bt_model(fixed_data, engine = "brglm2", verbose = FALSE)
summarize_bt_fit(firth_fit)
#> # A tibble: 6 × 6
#>   ID      theta    se  rank engine reliability
#>   <chr>   <dbl> <dbl> <int> <chr>        <dbl>
#> 1 a     -0.520  0.402     5 brglm2       0.570
#> 2 b     -0.348  0.395     4 brglm2       0.570
#> 3 c     -0.520  0.402     6 brglm2       0.570
#> 4 d     -0.0134 0.390     3 brglm2       0.570
#> 5 e      0.323  0.396     2 brglm2       0.570
#> 6 f      1.08   0.455     1 brglm2       0.570
firth_fit$provenance[c("engine_version", "adjustment", "identification")]
#> $engine_version
#> [1] "1.1.1"
#> 
#> $adjustment
#> $adjustment$method
#> [1] "firth"
#> 
#> $adjustment$type
#> [1] "AS_mean"
#> 
#> $adjustment$log_determinant_multiplier
#> [1] 0.5
#> 
#> 
#> $identification
#> $identification$convention
#> [1] "sum_to_zero"
#> 
#> $identification$internal_reference
#> [1] "f"
#> 
#> $identification$transformation
#>    contrast1  contrast2  contrast3  contrast4  contrast5
#> a  0.8333333 -0.1666667 -0.1666667 -0.1666667 -0.1666667
#> b -0.1666667  0.8333333 -0.1666667 -0.1666667 -0.1666667
#> c -0.1666667 -0.1666667  0.8333333 -0.1666667 -0.1666667
#> d -0.1666667 -0.1666667 -0.1666667  0.8333333 -0.1666667
#> e -0.1666667 -0.1666667 -0.1666667 -0.1666667  0.8333333
#> f -0.1666667 -0.1666667 -0.1666667 -0.1666667 -0.1666667

fixed_fit <- fit_bt_model(fixed_data, engine = "alpha", alpha = 0.5, verbose = FALSE)
summarize_bt_fit(fixed_fit)
#> # A tibble: 6 × 6
#>   ID      theta    se  rank engine reliability
#>   <chr>   <dbl> <dbl> <int> <chr>        <dbl>
#> 1 a     -0.535  0.403     6 alpha        0.590
#> 2 b     -0.358  0.396     4 alpha        0.590
#> 3 c     -0.535  0.403     5 alpha        0.590
#> 4 d     -0.0136 0.391     3 alpha        0.590
#> 5 e      0.333  0.398     2 alpha        0.590
#> 6 f      1.11   0.460     1 alpha        0.590
fixed_fit$provenance[c("engine_version", "package_version", "adjustment")]
#> $engine_version
#> [1] "4.6.1"
#> 
#> $package_version
#> [1] "1.6.0"
#> 
#> $adjustment
#> $adjustment$method
#> [1] "hamilton_alpha"
#> 
#> $adjustment$alpha
#> [1] 0.5
#> 
#> $adjustment$pseudo_count
#> [1] 0.1
#> 
#> $adjustment$pairs
#> [1] "all_unordered_pairs"
#> 
#> $adjustment$objective
#> [1] "log_likelihood + alpha/(N-1) * sum(log(p*(1-p)))"
#> 
#> $adjustment$reference
#> [1] "doi:10.1111/jedm.70022, equation (3)"

This alpha fit provides a common example that needs no additional fitting engine. The penalty is chosen in advance for demonstration; it is not selected using the reference or held-out outcomes.

Internal reliability: SSR and Bayesian EAP

Scale Separation Reliability (SSR)

Let $\hat\theta_i$ and $SE_i$ be estimates and standard errors on the same scale for the same $N$ items. Conventional SSR uses

$V_{\mathrm{observed}}=\frac{1}{N-1}\sum_i(\hat\theta_i-\bar{\hat\theta})^2, \qquad V_{\mathrm{error}}=\frac{1}{N}\sum_i SE_i^2,$ $\widehat V_{\mathrm{true}}=V_{\mathrm{observed}}-V_{\mathrm{error}}, \qquad SSR=\frac{\widehat V_{\mathrm{true}}}{V_{\mathrm{observed}}} =1-\frac{V_{\mathrm{error}}}{V_{\mathrm{observed}}}.$

In words, SSR subtracts average squared uncertainty from observed score spread, then expresses the remainder as a proportion of that spread. “Estimated true-score variance” is the name of this subtraction; it does not use known truth or prove that the remainder equals empirical recoverable variance.

ssr_components <- scale_separation_reliability(fixed_fit$theta$theta, fixed_fit$theta$se)
ssr_components
#> $observed_variance
#> [1] 0.4087629
#> 
#> $mean_squared_se
#> [1] 0.1675214
#> 
#> $true_score_variance
#> [1] 0.2412415
#> 
#> $ssr
#> [1] 0.5901746
#> 
#> $n_items
#> [1] 6
#> 
#> $n_finite
#> [1] 6
#> 
#> $valid
#> [1] TRUE
#> 
#> $status
#> [1] "ok"

The helper uses sample variance across items and the mean of squared SEs, not the square of their mean. It does not align names for you. Negative SSR and negative estimated true-score variance are retained: estimated uncertainty then exceeds observed spread. Invalid inputs are not repaired by dropping items.

Fit or calculation	SSR availability
sirt	Independently computed and checked against the engine’s `mle.rel`; missing/invalid SEs or undefined arithmetic cause an error.
Firth or alpha	Available with valid uncertainty and positive score variance; a valid equal-score fit is retained with `zero_score_variance` and missing reliability.
BradleyTerry2	Unavailable (`unavailable_se_convention`): its reference-based SEs are not used for centered-score SSR by this wrapper.
Experimental lapse model	Always unavailable (`prototype_ssr_unavailable`), for both interior and boundary fits.
Standalone SSR helper	Requires at least two finite scores with positive variance and matching finite, nonnegative SEs; otherwise errors.

Do not combine adjusted scores with another estimator’s SEs, or substitute bootstrap SDs to manufacture a “corrected SSR.” Report the components, estimator, SE convention, item set, and status along with the coefficient.

Bayesian BTL and EAP reliability

An expected a posteriori (EAP) score is an item’s posterior mean. If

$m_i=E(\theta_i\mid\text{data}),\qquad v_i=\operatorname{Var}(\theta_i\mid\text{data}),$

the implemented Bayesian coefficient is

$R_{\mathrm{EAP}}= \frac{\operatorname{var}_i(m_i)} {\operatorname{var}_i(m_i)+N^{-1}\sum_i v_i}.$

It describes posterior mean spread relative to that spread plus average posterior uncertainty. Unlike SSR, it adds uncertainty to the denominator instead of subtracting it from observed variance. The variances across items and across posterior draws use R’s sample-variance convention. Invalid draws, fewer than two items or draws, or nonpositive between-item variance give missing reliability; otherwise the result is clamped to $[0,1]$ .

The posterior model and priors affect both means and variances. Inspect sampling diagnostics before interpreting this coefficient. It does not measure empirical score recovery, held-out accuracy, or all sources of judgment error.

The following optional example fits the same fixed training outcomes with the simple btl likelihood. It requires installed CmdStan, cmdstanr, parallelly, and a compiler; it is not evaluated during builds. The seed and short sampler settings demonstrate reproducibility, not an adequate sampling budget.

bayes_input <- build_btl_results_data(data.frame(
  ID1 = fixed_data$object1, ID2 = fixed_data$object2,
  better_id = ifelse(fixed_data$result == 1, fixed_data$object1, fixed_data$object2)
), backend = "synthetic", model = "fixed_binary_bt")
sampler_settings <- list(
  chains = 2L, parallel_chains = 1L, iter_warmup = 250L,
  iter_sampling = 250L, seed = 30604L
)
bayes_fit <- fit_bayes_btl_mcmc(
  bayes_input, ids = ids, model_variant = "btl", cmdstan = sampler_settings
)
list(model_variant = "btl", selection = example_design$selection[1L],
     sampler = sampler_settings, package_version = as.character(packageVersion("pairwiseLLM")))
summarize_refits(bayes_fit)[, c(
  "round_id", "total_pairs", "diagnostics_pass", "reliability_EAP",
  "max_rhat", "min_ess_bulk", "divergences"
)]

The Bayesian BTL guide explains setup, other likelihood variants, item summaries, and diagnostics. The adaptive design guide explains how EAP reliability participates in stopping. Passing that operational rule is not independent validation of the scores.

Collect an adaptive schedule

Every strategy begins with seeded comparisons connecting the items. Later selection uses the current state:

`pairing_strategy`	Later comparisons
`random`	Choose a minimum-degree focal item, then a random allowed partner.
`trueskill_p50`	Choose a minimum-degree focal item, then the partner closest to a TrueSkill win probability of one half (maximum ambiguity).
`trueskill_pollitt`	Choose a minimum-degree focal item, then the partner closest to a TrueSkill win probability of one third or two thirds. This is Pollitt-inspired, not the original BTL-probability procedure.
`hybrid` (default)	Mix anchor, long, mid, and local comparisons with exploration, coverage, and exposure constraints.

Degree means the number of completed comparisons involving an item. Direct strategies (the first three rows) balance presentation and reverse a repeated pair; they allow at most two judgments per pair. They can stop when the chosen focal item has no legal partner. Random here is therefore a specific constrained policy, not uniform sampling from all pairs. A sparse replay reservoir adds its own restrictions; this example uses ordinary within-set selection.

TrueSkill provides fast step-by-step selection estimates. Bayesian refits supply posterior scores, diagnostics, reliability, and stopping checks; their global_identified signal also affects hybrid routing. A later frequentist fit does not replace the scheduling model or retroactively alter collected pairs.

Hybrid recovery stays within hybrid

When an attempt cannot select a pair, the canonical ladder is:

Attempt	What changes
`base`	Use normal stage bounds and duplicate policy.
`expand_locality`	Widen mid/local bounds, retaining long-link semantics.
`uncertainty_pool`	Keep base bounds and increase exploration.
`dup_relax`	Allow the conditional third judgment when the configured ceiling permits it; ambiguity and utility conditions still apply.
`global_safe`	Use the broadest stage-compatible domain with the ordinary two-judgment limit.

Stage identity and hard constraints remain in force. If exploration misses every pair but a previous attempt had eligible candidates, the usual exploitation rule can select from the last viable pool. These are recovery attempts within hybrid, not substitution of random, p50, or Pollitt strategies. Exhaustion can stop collection before the requested budget. Inspect fallback_used, fallback_path, and candidate_starved in adaptive_step_log() and the optional report from summarize_adaptive(state, include_starvation = TRUE). Direct strategies do not invoke this ladder. See the design guide for the exact rules.

A small provider-free p50 run

Save the pristine initial state before collecting outcomes. The synthetic judge below supplies fresh winners. Its function name does not make a provider call: adaptive_rank_run_live() also accepts a local judge like this one.

initial <- adaptive_rank_start(ids, seed = 17L,
  adaptive_config = list(pairing_strategy = "trueskill_p50"))
judge <- function(A, B, state, ...) {
  p <- plogis(truth[[A$item_id]] - truth[[B$item_id]])
  list(is_valid = TRUE, Y = as.integer(runif(1) < p))
}
schedule_config <- list(refit_pairs_target = 100L)
observed <- withr::with_seed(7304L, adaptive_rank_run_live(
  initial, judge, n_steps = 12L, btl_config = schedule_config, progress = "none"
))
adaptive_data <- adaptive_results_history(observed)
adaptive_fit <- fit_bt_model(adaptive_data, engine = "alpha", alpha = 0.5, verbose = FALSE)
selection_log <- adaptive_step_log(observed)
list(strategy = initial$controller$pairing_strategy, selector_seed = 17L,
     outcome_seed = 7304L, scheduling = schedule_config,
     committed_comparisons = nrow(adaptive_data), estimator = adaptive_fit$provenance$adjustment)
#> $strategy
#> [1] "trueskill_p50"
#> 
#> $selector_seed
#> [1] 17
#> 
#> $outcome_seed
#> [1] 7304
#> 
#> $scheduling
#> $scheduling$refit_pairs_target
#> [1] 100
#> 
#> 
#> $committed_comparisons
#> [1] 12
#> 
#> $estimator
#> $estimator$method
#> [1] "hamilton_alpha"
#> 
#> $estimator$alpha
#> [1] 0.5
#> 
#> $estimator$pseudo_count
#> [1] 0.1
#> 
#> $estimator$pairs
#> [1] "all_unordered_pairs"
#> 
#> $estimator$objective
#> [1] "log_likelihood + alpha/(N-1) * sum(log(p*(1-p)))"
#> 
#> $estimator$reference
#> [1] "doi:10.1111/jedm.70022, equation (3)"
selection_log[, c("A_id", "B_id", "Y", "pairing_strategy", "fallback_used")]
#> # A tibble: 12 × 5
#>    A_id  B_id      Y pairing_strategy fallback_used
#>    <chr> <chr> <int> <chr>            <chr>        
#>  1 a     b         0 trueskill_p50    warm_start   
#>  2 f     a         1 trueskill_p50    warm_start   
#>  3 c     f         0 trueskill_p50    warm_start   
#>  4 e     c         1 trueskill_p50    warm_start   
#>  5 d     e         0 trueskill_p50    warm_start   
#>  6 c     d         0 trueskill_p50    base         
#>  7 b     f         0 trueskill_p50    base         
#>  8 a     c         0 trueskill_p50    base         
#>  9 b     d         0 trueskill_p50    base         
#> 10 f     e         0 trueskill_p50    base         
#> 11 c     a         1 trueskill_p50    base         
#> 12 d     f         0 trueskill_p50    base

adaptive_results_history() already returns the three columns needed by fit_bt_model() and defaults to committed comparisons only. It needs no further build_bt_data() conversion. This demonstration specifies a refit interval larger than its entire budget, so it has no Bayesian refit or EAP coefficient. It does not demonstrate statistical stopping. In a design with scheduled refits, preserve those refits, including their optional CmdStan requirements. The practical adaptive guide covers that full lifecycle.

Empirical recovery and scale dispersion

In a simulation we know the generating scores. With real data, a dense or high-information fit can instead serve as a reference, but it is not known truth. Its uncertainty, estimator, overlap with training outcomes, and model assumptions affect apparent recovery. Here the dense reference uses independent judgments. Align IDs and center scores before comparing differences; do not rescale each vector to unit SD if scale expansion is the question.

reference_fit <- fit_bt_model(reference_data, engine = "alpha", alpha = 0.5, verbose = FALSE)
theta_adaptive <- setNames(adaptive_fit$theta$theta, adaptive_fit$theta$ID)[ids]
theta_reference <- setNames(reference_fit$theta$theta, reference_fit$theta$ID)[ids]
score_recovery <- function(estimate, reference) {
  data.frame(
    pearson_r_squared = cor(estimate, reference)^2,
    rank_correlation = cor(estimate, reference, method = "spearman"),
    score_sd = sd(estimate), reference_sd = sd(reference),
    sd_ratio = sd(estimate) / sd(reference),
    centered_rmse = sqrt(mean(((estimate - mean(estimate)) -
      (reference - mean(reference)))^2))
  )
}
recovery <- rbind(
  against_truth = score_recovery(theta_adaptive, truth[ids]),
  against_dense_estimate = score_recovery(theta_adaptive, theta_reference)
)
recovery
#>                        pearson_r_squared rank_correlation score_sd reference_sd
#> against_truth                  0.8500247        0.9428571 1.790407    0.7483315
#> against_dense_estimate         0.8309166        0.9428571 1.790407    0.8672850
#>                        sd_ratio centered_rmse
#> against_truth          2.392532      1.038837
#> against_dense_estimate 2.064382      0.969044
reference_fit$provenance$adjustment
#> $method
#> [1] "hamilton_alpha"
#> 
#> $alpha
#> [1] 0.5
#> 
#> $pseudo_count
#> [1] 0.1
#> 
#> $pairs
#> [1] "all_unordered_pairs"
#> 
#> $objective
#> [1] "log_likelihood + alpha/(N-1) * sum(log(p*(1-p)))"
#> 
#> $reference
#> [1] "doi:10.1111/jedm.70022, equation (3)"

Correlation describes association, not equality. Multiplying every score by a positive constant leaves ranks and Pearson correlation unchanged, so it also leaves squared Pearson correlation unchanged. This is the meaning of “R-squared” here; a prediction-error statistic such as $1-SSE/SST$ need not be invariant. The score SD and variance do change, and BT differences become larger log odds. Thus an adaptive schedule could stretch the scale while preserving strong rank recovery. Examine spread separately instead of using correlation as a test for the absence of scale distortion.

scale_comparison <- rbind(
  reference = score_recovery(theta_reference, theta_reference),
  doubled = score_recovery(2 * theta_reference, theta_reference)
)
scale_comparison
#>           pearson_r_squared rank_correlation score_sd reference_sd sd_ratio
#> reference                 1                1 0.867285     0.867285        1
#> doubled                   1                1 1.734570     0.867285        2
#>           centered_rmse
#> reference     0.0000000
#> doubled       0.7917192

Doubling here is an algebraic illustration, not a new estimator. This small example is not a benchmark establishing which selector or estimator performs best. Neither recovery table estimates calibration of SSR or EAP reliability; assessing that requires an appropriate repeated-sampling design and reference target.

Held-out predictive evaluation

Score recovery compares item scores. Predictive evaluation compares probabilities with new binary judgments. For first-item outcomes $y_k$ and probabilities $p_k$ , use the following summaries (lower loss is better, higher accuracy is better):

$\mathrm{log\ loss}=-\frac1M\sum_k[y_k\log p_k+(1-y_k)\log(1-p_k)], \qquad \mathrm{Brier}=\frac1M\sum_k(p_k-y_k)^2.$

Pairwise accuracy counts correct winners at a stated probability threshold. It discards the strength of the probability: 0.51 and 0.99 choose the same winner but express very different confidence. Log loss particularly penalizes confident mistakes. Brier score also uses probabilities. These are overall predictive scores, not complete calibration diagnostics.

p_adaptive <- predict(adaptive_fit, newdata = heldout_data[c("object1", "object2")])
p_reference <- predict(reference_fit, newdata = heldout_data[c("object1", "object2")])
p_expanded <- plogis(2 * qlogis(p_reference))
prediction_metrics <- function(p, y) {
  data.frame(
    log_loss = -mean(ifelse(y == 1, log(p), log1p(-p))),
    brier = mean((p - y)^2),
    accuracy = mean((p >= 0.5) == y)
  )
}
predictive_results <- rbind(
  adaptive = prediction_metrics(p_adaptive, heldout_data$result),
  dense_reference = prediction_metrics(p_reference, heldout_data$result),
  doubled_reference = prediction_metrics(p_expanded, heldout_data$result)
)
predictive_results
#>                    log_loss     brier  accuracy
#> adaptive          0.8177778 0.2607401 0.6583333
#> dense_reference   0.6097954 0.2101448 0.6750000
#> doubled_reference 0.7423706 0.2387389 0.6750000

The example classifies exactly 0.5 as an A win. It does not clip probabilities: a probability-zero realized event gives infinite log loss. Predictions here are plug-in probabilities from the alpha fits. The doubled reference retains the same predicted winners while moving probabilities farther from one half. Held-out losses can expose overconfidence that rank correlation or accuracy misses; no single finite test set guarantees this ordering.

Keep evaluation outcomes out of adaptive state updates, reference fitting, penalty selection, and stopping decisions. If you tune a method using evaluation results, reserve another independent test set for the final assessment. Sharing item identities is intentional here, but copying the same judgment into training and evaluation would not be independent evaluation.

Schedule-aware parametric bootstrap

bootstrap_bt_model() simulates new binary BT outcomes from the initial fitted scores, refits the chosen estimator, aligns item IDs, and centers scores. For each item it computes

$\widehat{\mathrm{bias}}_i=\overline{\theta_i^*}-\hat\theta_i, \qquad \hat\theta_{i,\mathrm{corrected}}=\hat\theta_i-\widehat{\mathrm{bias}}_i.$

This estimates point-estimate bias under the fitted generating model and the specified design. Corrections can change ranks and need not improve recovery.

Fixed versus adaptive scheduling

Fixed mode retains an outcome-independent schedule’s ordered rows and repeats, and draws fresh outcomes. It is appropriate for the training pairs chosen in advance above:

fixed_boot <- bootstrap_bt_model(fixed_fit, mode = "fixed", n_rep = 20L,
                                 seed = 1304L, keep = "theta")
fixed_boot$theta
#> # A tibble: 6 × 7
#>   ID    theta_initial bootstrap_mean    bias theta_corrected bootstrap_sd
#>   <chr>         <dbl>          <dbl>   <dbl>           <dbl>        <dbl>
#> 1 a           -0.535          -0.622 -0.0862         -0.449         0.523
#> 2 b           -0.358          -0.479 -0.120          -0.238         0.384
#> 3 c           -0.535          -0.454  0.0813         -0.617         0.515
#> 4 d           -0.0136         -0.124 -0.110           0.0966        0.392
#> 5 e            0.333           0.425  0.0915          0.242         0.484
#> 6 f            1.11            1.25   0.144           0.966         0.573
#> # ℹ 1 more variable: mcse_bias <dbl>
fixed_boot$provenance[c("mode", "seed", "pairing_strategy", "refit", "versions")]
#> $mode
#> [1] "fixed"
#> 
#> $seed
#> [1] 1304
#> 
#> $pairing_strategy
#> [1] "fixed"
#> 
#> $refit
#> $refit$estimator
#> [1] "alpha"
#> 
#> $refit$args
#> $refit$args$alpha
#> [1] 0.5
#> 
#> 
#> 
#> $versions
#>           R pairwiseLLM       stats 
#>     "4.6.1"     "1.6.0"     "4.6.1"

An adaptive schedule depends on earlier outcomes. Its bootstrap must repeat the selection decisions after simulating those outcomes. Freezing the observed adaptive pairs conditions on a realized schedule and answers a different question. The package uses the canonical runner, not a separate approximation to p50, Pollitt-inspired, random, or hybrid selection:

adaptive_boot <- bootstrap_bt_model(adaptive_fit, mode = "adaptive", n_rep = 4L,
  seed = 2304L, initial_state = initial, budget = 12L,
  btl_config = schedule_config, keep = "theta")
adaptive_boot$theta
#> # A tibble: 6 × 7
#>   ID    theta_initial bootstrap_mean    bias theta_corrected bootstrap_sd
#>   <chr>         <dbl>          <dbl>   <dbl>           <dbl>        <dbl>
#> 1 a            -2.61          -2.14   0.469           -3.08         0.601
#> 2 b            -0.870         -0.197  0.674           -1.54         0.571
#> 3 c            -0.715         -0.648  0.0667          -0.782        0.681
#> 4 d             0.397         -0.453 -0.851            1.25         1.17 
#> 5 e             2.44           1.71  -0.724            3.16         0.489
#> 6 f             1.36           1.73   0.366            0.993        0.395
#> # ℹ 1 more variable: mcse_bias <dbl>
adaptive_boot$replicates[, c(
  "replicate", "selector", "outcome", "scheduling_refit", "estimator",
  "success", "n_comparisons", "schedule_digest"
)]
#> # A tibble: 4 × 8
#>   replicate  selector   outcome scheduling_refit estimator success n_comparisons
#>       <int>     <int>     <int>            <int>     <int> <lgl>           <int>
#> 1         1 156533677 156533778        156533879 156533980 TRUE               12
#> 2         2 156543684 156543785        156543886 156543987 TRUE               12
#> 3         3 156553691 156553792        156553893 156553994 TRUE               12
#> 4         4 156563698 156563799        156563900 156564001 TRUE               12
#> # ℹ 1 more variable: schedule_digest <chr>
adaptive_boot$provenance[c("mode", "seed", "pairing_strategy", "refit", "versions")]
#> $mode
#> [1] "adaptive"
#> 
#> $seed
#> [1] 2304
#> 
#> $pairing_strategy
#> [1] "trueskill_p50"
#> 
#> $refit
#> $refit$estimator
#> [1] "alpha"
#> 
#> $refit$args
#> $refit$args$alpha
#> [1] 0.5
#> 
#> 
#> 
#> $versions
#>           R pairwiseLLM       stats 
#>     "4.6.1"     "1.6.0"     "4.6.1"

The pristine initial connected pairs and presentations remain fixed; subsequent selector randomness and simulated outcomes use deterministic per-replicate seeds. Separate seeds cover selection, outcomes, scheduling refits, and the final estimator. Compatible software and deterministic callbacks are needed for exact reproduction. Serial and optional parallel execution use the same replicate seeds and aggregation order.

Keep the original scheduling configuration. This example deliberately has no scheduled refits within its budget; a design that uses Bayesian refits must rerun them. In hybrid, those refits can affect later routing and tapering. The final bootstrap estimator is separate from the scheduling fitter. The bootstrap budget counts committed comparisons including initialization; statistical stopping does not shorten that fixed-budget experiment, but terminal selection exhaustion fails a replicate. Hybrid recovery remains hybrid, including inside bootstrap.

What the Monte Carlo summaries mean

Output	Meaning
`bias`, `theta_corrected`	Estimated bias and corrected point estimate under the fitted model and specified schedule semantics.
`bootstrap_sd`	Dispersion of the successfully refitted replicate estimates.
`mcse_bias`	Simulation error in the estimated bias, calculated as `bootstrap_sd / sqrt(n_success)`.
`n_success`, `n_failed`, `replicates`	Whether the correction has sufficient successful simulations and what happened in each replicate.

Neither bootstrap SD nor Monte Carlo error is a corrected-score standard error. The API does not provide corrected-score confidence intervals or corrected SSR. These small replicate counts keep the example fast; their Monte Carlo error can be substantial and they are not study recommendations. The default requires all replicates to succeed. A deliberately lower min_success excludes failures and therefore changes the summary to successful refits; inspect every failure. An unmet threshold raises an error carrying condition$result diagnostics and missing corrected scores.

The built-in final estimators are alpha and Firth; other supported estimators require an explicit callback. The generating model is simple binary BT, without lapse, position, or ties. Lapse fits are rejected. Keep full fit/bootstrap objects to retain settings and provenance. The bootstrap guide explains callback requirements, optional parallel workers, and artifact retention.

Experimental lapse and positional-bias model matching

fit_bt_model(..., engine = "lapse") uses the Bayesian btl_e_b likelihood with frequentist unpenalized maximum likelihood:

$P(A\text{ wins})=(1-\epsilon) \operatorname{logit}^{-1}(\theta_A-\theta_B+\beta)+\epsilon/2.$

The item strengths sum to zero; positive positional beta favors the first presented item; lapse epsilon is the probability of guessing. Unlike conventional sirt epsilon, lapse epsilon is an estimated likelihood parameter. The estimator jointly fits all three parameter groups without Bayesian priors. It helps separate differences in probability models from differences in estimation; it is experimental, not the default BT estimator or a Bayesian replacement.

The following separate four-item example uses both presentation directions and many inexpensive simulated binary rows to identify the extra parameters. The counts are demonstration choices, not guidance on the data needed in practice.

lapse_ids <- letters[1:4]
lapse_truth <- setNames(seq(-2, 2, length.out = 4), lapse_ids)
lapse_pairs <- expand.grid(object1 = lapse_ids, object2 = lapse_ids, stringsAsFactors = FALSE)
lapse_pairs <- lapse_pairs[lapse_pairs$object1 != lapse_pairs$object2, ]
lapse_pairs <- lapse_pairs[order(lapse_pairs$object1, lapse_pairs$object2), ]
fit_lapse_example <- function(beta, epsilon, seed) {
  p <- (1 - epsilon) * plogis(lapse_truth[lapse_pairs$object1] -
    lapse_truth[lapse_pairs$object2] + beta) + epsilon / 2
  wins <- withr::with_seed(seed, rbinom(nrow(lapse_pairs), 2000L, p))
  dat <- lapse_pairs[rep(seq_len(nrow(lapse_pairs)), each = 2000L), ]
  dat$result <- unlist(lapply(wins, function(w) c(rep(1L, w), rep(0L, 2000L - w))))
  fit_bt_model(dat, engine = "lapse", verbose = FALSE)
}
lapse_interior <- fit_lapse_example(beta = 0.3, epsilon = 0.2, seed = 305L)
lapse_boundary <- fit_lapse_example(beta = 0, epsilon = 0, seed = 30501L)
lapse_examples <- list(interior = lapse_interior, boundary = lapse_boundary)
lapply(lapse_examples, function(fit) {
  list(engine = fit$engine, beta = fit$beta, epsilon = fit$epsilon,
       convergence = fit$provenance$convergence$status,
       uncertainty = fit$provenance$uncertainty, ssr = fit$ssr)
})
#> $interior
#> $interior$engine
#> [1] "lapse"
#> 
#> $interior$beta
#> [1] 0.3010381
#> 
#> $interior$epsilon
#> [1] 0.2117678
#> 
#> $interior$convergence
#> [1] "converged"
#> 
#> $interior$uncertainty
#> $interior$uncertainty$method
#> [1] "inverse_joint_observed_information"
#> 
#> $interior$uncertainty$valid
#> [1] TRUE
#> 
#> $interior$uncertainty$status
#> [1] "valid"
#> 
#> $interior$uncertainty$coordinates
#> [1] "sum_to_zero_theta_beta_epsilon"
#> 
#> $interior$uncertainty$scope
#> [1] "conditional_on_realized_comparison_graph"
#> 
#> $interior$uncertainty$schedule_aware
#> [1] FALSE
#> 
#> 
#> $interior$ssr
#> $interior$ssr$ssr
#> [1] NA
#> 
#> $interior$ssr$valid
#> [1] FALSE
#> 
#> $interior$ssr$status
#> [1] "prototype_ssr_unavailable"
#> 
#> 
#> 
#> $boundary
#> $boundary$engine
#> [1] "lapse"
#> 
#> $boundary$beta
#> [1] -0.02964059
#> 
#> $boundary$epsilon
#> [1] 0
#> 
#> $boundary$convergence
#> [1] "converged_boundary"
#> 
#> $boundary$uncertainty
#> $boundary$uncertainty$method
#> [1] "none"
#> 
#> $boundary$uncertainty$valid
#> [1] FALSE
#> 
#> $boundary$uncertainty$status
#> [1] "nonregular_boundary"
#> 
#> $boundary$uncertainty$coordinates
#> [1] "sum_to_zero_theta_beta_epsilon"
#> 
#> $boundary$uncertainty$scope
#> [1] "conditional_on_realized_comparison_graph"
#> 
#> $boundary$uncertainty$schedule_aware
#> [1] FALSE
#> 
#> 
#> $boundary$ssr
#> $boundary$ssr$ssr
#> [1] NA
#> 
#> $boundary$ssr$valid
#> [1] FALSE
#> 
#> $boundary$ssr$status
#> [1] "prototype_ssr_unavailable"
list(selection = "all ordered pairs chosen before outcomes",
     outcome_seeds = c(interior = 305L, boundary = 30501L),
     model_variant = lapse_interior$model_variant,
     adjustment = lapse_interior$provenance$adjustment,
     package_version = lapse_interior$provenance$package_version)
#> $selection
#> [1] "all ordered pairs chosen before outcomes"
#> 
#> $outcome_seeds
#> interior boundary 
#>      305    30501 
#> 
#> $model_variant
#> [1] "btl_e_b"
#> 
#> $adjustment
#> $adjustment$method
#> [1] "none"
#> 
#> $adjustment$log_penalty
#> [1] 0
#> 
#> 
#> $package_version
#> [1] "1.6.0"
predict(lapse_boundary, newdata = lapse_pairs)
#>  [1] 0.21486917 0.06826458 0.01805636 0.77496115 0.20628310 0.06123446
#>  [7] 0.92786699 0.78384172 0.19591582 0.98086198 0.93526782 0.79457685

A valid interior MLE has $0<\epsilon<1$ and passes identification, stationarity, and information/Hessian checks. Its joint covariance includes estimation of beta and epsilon; item SEs are local model-based approximations conditional on the realized graph, not a claim of repeated-sample calibration.

A valid epsilon-zero boundary MLE needs no guessing component. It retains theta, beta, likelihoods, and predictions, with convergence status converged_boundary. Ordinary joint interior uncertainty is unavailable: theta$se is missing, vcov and parameter_vcov are NULL, and uncertainty status is nonregular_boundary. A small positive epsilon alone does not establish a boundary optimum. Epsilon one leaves strengths and beta unidentified and is rejected. Connectedness alone is insufficient for this model’s identification, and genuine numerical/identification failures still error.

Conventional SSR is unavailable for both interior and boundary lapse fits. Do not compute it manually from the interior SEs, or feed a lapse fit into the simple-BT bootstrap. See fit_bt_model() and predict.pairwiseLLM_bt_lapse() for the full experimental contract and failure diagnostics.

Reporting the complete workflow

Report which strategy collected the comparisons, which estimator fit them, its settings and uncertainty convention, and which internal reliability formula was used. Then report recovery against the stated reference separately from held-out predictive performance. Include score spread, evaluation-data provenance, seeds, and software versions. Retain selection logs and complete fit/bootstrap objects.

The examples illustrate those distinctions. They do not establish universal estimator or strategy superiority, universal reliability calibration, suitability for a particular study, or equivalence of internal reliability and empirical recovery.

Citation

Mercer, S. H. (2026). Adaptive comparative judgment: estimation, reliability, and evaluation [R package vignette]. Comprehensive R Archive Network. https://doi.org/10.32614/CRAN.package.pairwiseLLM

Five related questions