---
title: "Rasch analysis of comparative judgements"
author: "Joshua A. McGrane"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Rasch analysis of comparative judgements}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 4.2)
options(digits = 4)
source("precomputed.R")
```

```{r library}
library(rasch)
```

## Fit the Bradley--Terry--Luce model

Paired-comparison data record two objects and an observed preference. In the
Bradley--Terry--Luce model (Bradley and Terry 1952; Luce 1959), the log odds
of choosing object A over object B are their location difference:

$$
P(A\succ B)=\frac{\exp(\beta_A)}
{\exp(\beta_A)+\exp(\beta_B)}.
$$

This is the conditional form of the dichotomous Rasch model (Rasch 1960;
Andrich 1978). The comparison graph must connect all objects;
otherwise their relative locations are not identified.

For ordered comparisons, `btl` fits the adjacent-category extension

$$
\log\frac{P(Y=r)}{P(Y=r-1)}=\beta_A-\beta_B-\tau_r,
$$

with thresholds symmetric under reversal of presentation order (Tutz 1986).

```{r fit}
d <- simulate_btl(
  n_objects = 8,
  n_judges = 48,
  reps_per_pair = 84,
  erratic_judges = 2 / 48,
  dependence = list(exposure = 0.7, carry_over = 0),
  seed = 2
)
truth <- attr(d, "truth")
judge_number <- as.integer(sub("^J", "", d$judge))
d$panel <- factor(ifelse(judge_number %% 2L,
                         "panel A", "panel B"))
d$experience <- factor(ifelse(judge_number <= 24,
                              "experienced", "novice"))

fit <- btl(d, object_a = "object_a", object_b = "object_b",
           winner = "winner", judge = "judge", order = "order")
fit
fit$objects
```

Panel and experience are crossed judge factors with no planted object
difference. Two judges answer at random, an exposure effect is added in
judgment order, and carry-over remains zero. Residual dimensionality is
examined separately below because combining a second attribute with this
history model would make the carry-over estimate conditional on a deliberately
misspecified common scale.

Judge residuals describe agreement with the common object scale; they are not
person measures. Object fit, judge fit, targeting, and comparison information
address different parts of the design and should be considered together.

```{r judges}
judge_order <- order(abs(fit$judges$fit_resid), decreasing = TRUE)
head(fit$judges[judge_order, ], 6)
erratic_fit <- fit$judges[match(truth$erratic, fit$judges$judge), ]
erratic_judge <- erratic_fit$judge[
  which.max(abs(erratic_fit$fit_resid))
]
judge_surprise(fit, erratic_judge)
btl_information(fit)
```

## Bootstrap goodness of fit

The pairwise chi-square and the object and judge residuals use approximate
references. `fit_bootstrap()` generates outcomes on the observed comparison
design and refits the model. It reports the whole-model probability and joint
adjusted probabilities for pairs, objects and judges. Ordered thresholds and
history effects are retained; history-dependent outcomes are generated in
sequence.
Report the usable, non-converged, and other-failure counts. A maxT adjustment
is unavailable for the complete family if one testable member lacks a usable
joint null.
Each statistic is adjusted under the fitted global null. The adjustment does
not guarantee familywise error among fitting pairs, objects or judges when
another member departs from the model.

With exposure and carry-over fitted jointly at 30 judges, a 1,000-dataset
dichotomous study gave 4.8% total error and 3.9%, 4.8% and 5.5% adjusted
familywise error for pairs, objects and judges. A 200-dataset four-category
study gave 3.0%, 3.5%, 6.0% and 6.0%. All 238,800 refits were usable.

```{r btl-fit-bootstrap, eval = FALSE}
boot <- fit_bootstrap(fit, B = 999, seed = 2026)
boot$total
head(boot$pairs[order(boot$pairs$chisq_p_boot_adj), c(
  "object_a", "object_b", "chisq", "chisq_p_boot_adj"
)])
head(boot$objects[order(boot$objects$fit_resid_p_boot_adj), c(
  "object", "fit_resid", "fit_resid_p_boot_adj"
)])
head(boot$judges[order(boot$judges$fit_resid_p_boot_adj), c(
  "judge", "fit_resid", "fit_resid_p_boot_adj"
)])
```

## Check transitivity and residual structure

Circular triads (Kendall and Babington Smith 1940) identify local
contradictions in the observed ordering.
Residual dimensionality asks whether comparisons contain a structured second
attribute after the primary scale is fitted.

When a judgment-order column is supplied, `fit$dependence` reports the
exposure and carry-over effects; adding `position = TRUE` reports the position
effect alongside them. Use `p_adj`, which applies Holm's correction across the
declared effects even when one probability is withheld; `p` is retained as
the raw probability.

```{r structure}
tr <- btl_transitivity(fit)
tr

dim_data <- simulate_btl(
  n_objects = 8, n_judges = 48, reps_per_pair = 84,
  second_attribute = list(rho = 0.3), seed = 47
)
dim_fit <- btl(dim_data, object_a = "object_a", object_b = "object_b",
               winner = "winner", judge = "judge")
```

The dimensionality reference is precomputed in the distributed vignette.
The following call reproduces it.

```{r dimensions, eval = recompute}
dimensions <- btl_dimensionality(dim_fit, reps = 20, seed = 2026,
                                independent_comparisons = TRUE)
dimensions
```

```{r dimensions-precomputed, echo = FALSE}
if (!recompute) {
  dimensions <- vignette_result("paired-comparisons")$dimensions
  dimensions
}
```

```{r plot-transitivity, fig.alt = "Per-judge consistency of the paired comparisons."}
plot_btl_transitivity(tr)
```

```{r plot-scree, fig.alt = "Residual bimension strengths against the simulated noise reference."}
plot_btl_scree(dimensions)
```

The decomposition compares observed and fitted expected points within each
object pair. For ordered responses this includes the category thresholds;
fitted position and history effects also enter the expectation.

The simulated comparisons in this example are independent, so the call
explicitly requests that reference. For judge-clustered data, the default
shows only the residual decomposition: clustered repetition can look like
dimensional structure against an independent-comparison reference.
The shaded area runs from the simulated mean to the finite-simulation 5%
upper limit. Twenty replicates keep this example quick; a final sensitivity
analysis needs enough replicates to stabilise the reference distribution.
An ordered analysis needs one row per comparison:
a replication count does not retain the sequence needed to construct exposure
and carry-over histories, so `btl()` refuses that combination.

## Examine DIF across judge groups

`btl_dif` tests whether object locations differ across nominated judge factors.
The omnibus analysis uses judges as the independent units. A significant term
is resolved into factor-specific object locations and pairwise logit
differences. HC3 covariance allows the precision of judge means to vary with
their comparison workloads. Omnibus and pairwise inference require at least
eight judges and eight effective judges in every relevant factor cell;
estimates remain descriptive below that boundary. The tables report both
counts. Panel and experience are null in this example, so an adjusted flag
would be a false positive rather than a planted result.

```{r btl-dif-example}
bd <- btl_dif(fit, d[c("panel", "experience")])
bd_order <- order(bd$summary$p_uniform_adj)
head(bd$summary[bd_order, c(
  "object", "term", "F_uniform", "p_uniform_adj", "eta2_uniform",
  "min_judges", "min_effective_judges"
)], 6)
```

An optional fitted-outcome bootstrap repeats the comparison fit and the
complete object-by-term family. It is a sensitivity analysis beside the
primary HC3 residual analysis. Its familywise probabilities refer to the
fitted global invariant null.

```{r btl-dif-bootstrap, eval = FALSE}
bd_boot <- dif_bootstrap(fit, bd, B = 999, seed = 2026)
bd_boot$summary
```

## Equate panels through common objects

`btl_equate` aligns two calibrations that share at least two objects. Two
common objects identify a descriptive origin shift; at least three with
usable covariance information are needed for object-drift tests. For two
fitted calibrations, drift inference is withheld until independent judges and
comparisons are stated explicitly. The shift is precision-weighted when at
least two common objects have usable variances. If they do not, the function
uses the unweighted mean location difference, labels it in `shift_method`, and
keeps the link descriptive. An exact common anchor fixes the shift even when
other common objects have unavailable SEs.

```{r equating-example, eval = FALSE}
eq <- btl_equate(current_panel, reference_panel, independent = TRUE)
eq$table
eq$equated                 # reference panel on the current panel's origin
```

The SEs in `eq$equated` include uncertainty in the shift and its covariance
with the reference locations. They are withheld when that joint uncertainty
is unavailable. The table retains its covariance and finite-sample degrees
of freedom as attributes; it is not independent of either input calibration.

A bank table may be used in place of `reference_panel`. Marginal object
standard errors are not enough for drift tests because they omit the
covariance created by the bank's fitted origin. Attach the joint matrix as
`attr(bank, "cov_location")`, ordered like the bank rows or named by object;
otherwise the alignment is descriptive. A bank treated as fixed may instead
carry zero standard errors. Dependent panels require a joint or paired
bootstrap outside this function.

## Linked frames for paired comparisons

`btl_efrm` combines the comparative judgement model with Humphry's extended
frame of reference structure (Humphry and Andrich 2008). Judges belong to
panels, and objects belong to linked sets. For object $k$ in set $s$,

$$
v_k=\alpha_s\beta_k+\kappa_s.
$$

A same-set comparison in panel $g$ has logit
$\phi_g(\beta_A-\beta_B)$; a cross-set comparison has logit
$\phi_g(v_A-v_B)$. As in Humphry and Andrich (2008, eq. 15), $\phi_g$ and
$\alpha_s$ are unit ratios: the common reference unit over the frame's own
unit. The panel units $\phi_g$ have geometric mean one. For the linked item
sets, the first set fixes $\alpha_1=1$ and $\kappa_1=0$, defining the unit
and origin of the common scale; the remaining set transformations are
estimated relative to it. A value above one denotes a finer unit and steeper
comparisons. Cross-set comparisons identify the set units and origins.
The cross-set likelihood holds the within-set locations and panel units fixed.
It estimates the set transformations directly from the comparison outcomes
and does not use the finite-grid person-distribution link in `rasch_efrm()`.

```{r efrm}
de <- simulate_btl_efrm(
  n_objects_per_set = 5, n_sets = 2,
  n_judges_per_panel = 6, n_panels = 2,
  reps_within = 15, reps_cross = 15,
  set_units = c(1, 1.3), set_origins = c(0, 0.6), seed = 9
)
ef <- btl_efrm(
  de, "object_a", "object_b", winner = "winner", judge = "judge",
  panels = "panel", object_sets = attr(de, "truth")$object_sets,
  se_method = "conditional"
)
ef$phi_table
ef$alpha_table
ef$kappa_table
```

The default judge bootstrap resamples judges within panels and refits both
stages. A set's panel-ratio fit must pass an exact likelihood-curvature check
before entering the reconciliation. Failed sets are omitted from that
reconciliation; estimation stops if the remaining sets cannot link all panels.
Cross-set outcomes are checked for complete and quasi-complete separation;
neither supports a finite link. Reaching the iteration limit is not convergence.
Saved frame fits without a current likelihood-check record require refitting
before reopening in the app; their source files remain unchanged.
The parametric bootstrap (`se_method = "bootstrap"`) draws independent
outcomes from the fitted model. The conditional option used above reports
estimates and conditional standard errors but withholds probabilities because
it does not propagate stage-one uncertainty into the set link. With either
bootstrap, omnibus probabilities are Holm-adjusted across the three unit
families. Individual panel-unit, set-unit and set-origin contrasts form one
separate Holm-adjusted follow-up family. Judge-bootstrap tests require at least six judges
and 5.5 effective judges in every panel, and eight of each on a set link.
The pooled pairwise chi-square is retained as a descriptive fit summary, but
its row-based probability is not reported because judges contribute repeated
comparisons.
Judge resamples are distributed over four workers by default, or fewer where
the system imposes a lower limit. Set `seed` to reproduce the resamples. The
parametric bootstrap remains serial because its refits are inexpensive. In the
application, frame estimation runs in the background and may be cancelled.

With six judges per panel and 20 repetitions per pair, a 1,000-fit null study
gave raw marginal rejection of 3.3 per cent for panel units, 6.7 per cent for
set units and 6.0 per cent for set origins. Holm familywise rejection was 3.9
per cent across the omnibus decisions and 3.0 per cent across their follow-ups;
set-unit coverage was 0.900. This design lies in the caution band reported by
the fit. With 12 judges per panel and the default 200 resamples, raw set-unit
rejection was 4.6 per cent, the empirical-to-reported SE ratio was 0.992 and
coverage was 0.934 over 500 null fits.
With more within- and cross-set information, finite-sample attenuation of the
set-unit estimate declined: log-unit bias was -0.041 at 20 repetitions per
pair, -0.016 at 50 and -0.006 at 100. Reported decisions use Holm adjustment
across the three omnibus tests; the individual unit contrasts form a separate
Holm-adjusted family.

## The same analysis in the application

`rasch::run_app()` fits comparative judgement designs alongside the other
models, and the sequence above maps onto its panels. Choosing **Comparative
Judgement** as the model changes what the **Data** panel asks for: instead of
item columns it wants the two object columns and the observed preference,
with the judge column and any judge factors beside them. A judgment-order
column is assigned here too, which is what makes the exposure and carry-over
terms available later.

```{r app-cj-data, echo = FALSE, out.width = "100%", fig.alt = "The Data panel with the comparative judgement roles assigned: the two object columns, the observed preference, and the judge column."}
knitr::include_graphics("figures/app-cj-data.png")
```

**Summary** reports what `btl()` prints: the design counts, the pairwise
chi-square, the object separation index, and the within-judge exposure and
carry-over effects when a judgment-order column was assigned --- with
inference withheld on designs too small to support it, stated as a note
rather than silently. A design whose comparison graph does not connect is
refused for the same reason it is refused in the code --- the relative
locations are not identified. The fit-bootstrap button runs the fitted-design
bootstrap in the background and adds adjusted probabilities to the pair,
object and judge tables. It can be cancelled without losing the fitted model.

```{r app-cj-summary, echo = FALSE, out.width = "100%", fig.alt = "The Summary panel for a comparative judgement fit: the design counts, object separation and pairwise fit tiles, and the test-of-fit table."}
knitr::include_graphics("figures/app-cj-summary.png")
```

**Items** becomes the object panel. The caterpillar plot orders the objects
by location with their confidence intervals, which is the graphical form of
`fit$objects`; the panels beside it hold the symmetric thresholds, the
threshold components and the category probability curves for an ordered
design, and the pairwise fit that `btl_transitivity()` summarises.

```{r app-cj-items, echo = FALSE, out.width = "100%", fig.alt = "The object panel, showing the caterpillar plot of object locations with confidence intervals."}
knitr::include_graphics("figures/app-cj-items.png")
```

**Persons** becomes the judge panel: judge fit, the surprise index that
`judge_surprise()` returns, and the consistency of each judge's comparisons.
These describe agreement with the common object scale, and are not person
measures --- the same caution the code section makes.

```{r app-cj-judges, echo = FALSE, out.width = "100%", fig.alt = "The judge panel, showing judge fit and consistency against the common object scale."}
knitr::include_graphics("figures/app-cj-judges.png")
```

Under **Invariance**, the DIF panel runs `btl_dif()` over the nominated judge
factors. It reports the omnibus test with judges as the independent units,
the factor-specific object locations and their pairwise differences in
logits, and the raw and effective judge counts in the least-supported cells
--- so a design below the eight-judge boundary shows why its inference is
withheld rather than simply returning nothing.

```{r app-cj-dif, echo = FALSE, out.width = "100%", fig.alt = "The DIF panel for a comparative judgement fit, showing the omnibus test across judge groups with the judge counts per level."}
knitr::include_graphics("figures/app-cj-dif.png")
```

Equating and extended frames sit in the same menu. Frame estimation runs in
the background and can be cancelled, because the judge bootstrap is the
expensive part of the analysis. As in the Rasch panels, every table and
figure carries an **R code** disclosure holding the call that produced it.

## A worked analysis on real data

The [party blocs case
study](https://drjoshmcgrane.github.io/rasch/articles/articles/party-blocs-crisis.html)
fits the comparative judgement frame model to real paired comparisons
between political parties, with judge panels and ideological blocs as
frames. Its script ships with the package under `casestudies`.

## References

Andrich, D. (1978). Relationships between the Thurstone and Rasch approaches
to item scaling. *Applied Psychological Measurement*, 2, 451--462.

Bradley, R. A., and Terry, M. E. (1952). Rank analysis of incomplete block
designs: I. The method of paired comparisons. *Biometrika*, 39, 324--345.

Humphry, S. M., and Andrich, D. (2008). Understanding the unit in the
Rasch model. *Journal of Applied Measurement*, 9(3), 249--264.

Kendall, M. G., and Babington Smith, B. (1940). On the method of paired
comparisons. *Biometrika*, 31(3/4), 324--345.

Tutz, G. (1986). Bradley-Terry-Luce models with an ordered response.
*Journal of Mathematical Psychology*, 30(3), 306--316.

Luce, R. D. (1959). *Individual Choice Behavior: A Theoretical Analysis*.
Wiley.

Rasch, G. (1960). *Probabilistic Models for Some Intelligence and Attainment
Tests*. Copenhagen: Danish Institute for Educational Research. (Expanded
edition, 1980, Chicago: University of Chicago Press.)