Competitions and paired comparisons: Bradley-Terry and Elo

Module: Beyond the standard models
Tranche: extension (beyond a first course)
Before this: The Rasch model
Builds toward: none
IRW tables: nba_2012-2018, friedman2019_risk_harm, friedman2019_risk_disaster
From EDUC 252: n/a
Status: draft
Authorship: Written largely by Claude (Anthropic), from Ben Domingue’s EDUC 252 materials and under Ben’s direction. Ben reviews each lesson before it is marked done.
Comments: open an issue on GitHub about this lesson.

What this is for

Two basketball teams play; one wins. A survey respondent is shown two risks and picks the one that did more harm. Neither setting has items, only pairs set against each other, yet the model we need is one we already have. The Rasch model compares a person with an item. Replace the item with a second person and you have the Bradley-Terry model (Bradley & Terry, 1952), a model for games, tournaments and paired-comparison judgments alike. Its online cousin, the Elo rating, produced the chess ratings in The likelihood and logistic regression. This lesson builds both from the Rasch model, then asks a question item response data rarely force on us: what if the thing measured changes while we measure it?

Goals

By the end of this lesson you should be able to:

  • write down the Bradley-Terry model, fit it as a logistic regression, and say what fixes its origin and when its estimates exist;
  • explain why win percentage and Bradley-Terry strength agree when every agent faces the same opponents, and why they part ways when schedules differ;
  • derive the Elo update from the Bradley-Terry likelihood and say what its step size \(K\) trades off;
  • compare a static and a dynamic rating by how well each predicts games it hasn’t seen;
  • build a scale from paired-comparison judgments and say what the question asked does to it.

Core ideas

Competition data: no items, only opponents

In a competition each row pairs two agents, and the outcome says which came out ahead. The IRW stores such data in its competition standard: agent_a, agent_b, a winner (agent_a, agent_b or draw), a date, which side had the home advantage, and scores where there are scores; irw_fetch(name, source = "comp") loads them. Each outcome measures both agents at once, on one scale, and tells us only about the gap between them. The same shape fits judgments: a judge asked “which of these two is more harmful?” also produces a winner and a loser (Glickman & Jones, 2025, review models for both).

The Bradley-Terry model is a Rasch model in which the item is another agent

Recall from The Rasch model (Go deeper: sufficiency): take two items and condition on a respondent getting exactly one of them right. The probability that it was item \(i\) is \(1/(1 + e^{b_i - b_k})\), whatever the respondent’s \(\theta\). Swap the roles and compare two respondents on one item, and the comparison depends only on \(\theta_j - \theta_l\).

That comparison of two respondents is already a model for a competition. Write \(\theta_j\) for the strength of agent \(j\). The Bradley-Terry model says that when \(j\) meets \(k\), \[\Pr(j \text{ beats } k) = \frac{\exp(\theta_j - \theta_k)}{1 + \exp(\theta_j - \theta_k)}.\] Set it beside the Rasch model: the opponent’s strength \(\theta_k\) sits exactly where the item difficulty \(b_i\) sat. Zermelo (1929) proposed this model for tournament results; Bradley & Terry (1952) arrived at it for paired comparisons in experiments, and Luce (1959) reached the same form from a theory of choice behaviour. Andrich (1978) set out its relation to the Rasch model.

Anything else that tilts a game goes into the exponent. A home team has an advantage \(h\), so \(\Pr(\text{home } j \text{ beats } k) = \text{logistic}(h + \theta_j - \theta_k)\). In a judgment task, the same term captures a preference for whichever option is shown first.

Try it: the opponent is the item

Move the visitors’ strength and the curve slides without changing shape, as the Rasch curve did when we moved \(b\). The home advantage slides it too, as if every visitor were \(h\) logits weaker.

Two facts carry over from the Rasch model unchanged. The first is identification.

Recall from The Rasch model: adding the same constant to every \(\theta\) and every \(b\) leaves every probability unchanged, so the scale has no origin until we fix one.

Only differences \(\theta_j - \theta_k\) appear, so we fix one agent at 0 or make the strengths average 0. The second fact is that the model is a logistic regression: one row per game, outcome 1 if the first-listed agent won, and one column per agent, \(+1\) for the first agent, \(-1\) for the second and 0 otherwise. The intercept is \(h\). glm() fits it, as does the BradleyTerry2 package (Turner & Firth, 2012).

When does the maximum exist? Zermelo (1929) and Ford (1957) give the condition. Split the agents into any two groups; somebody in each group must have beaten somebody in the other. An agent who never loses fails it, and so does a league split into two halves that never meet.

Recall from The likelihood and logistic regression (Go deeper): if a predictor perfectly separates successes from failures, the log likelihood keeps rising as the slope grows and the MLE doesn’t exist (Albert & Anderson, 1984).

Try it: the undefeated team

Four teams play each other twice; team A won five games and lost one, to D. The plot follows A’s estimated strength as the fitting algorithm runs (Zermelo’s iteration, which climbs the likelihood at every step; Hunter, 2004).

With one loss, the estimate settles within a few dozen iterations. Without it, every round of the algorithm pushes A higher: the data are best explained by an A infinitely stronger than everyone else. A single draw, counted as half a win as chess scores it, brings the estimate back. The data must contain some evidence that A can be beaten.

Wins are to Bradley-Terry what the sum score is to Rasch

Recall from The likelihood and logistic regression: in a logistic regression, the derivative of the log likelihood with respect to a coefficient is \(\sum_i (y_i - p_i)x_i\). For the intercept this says that at the maximum, the predicted number of successes equals the observed number.

Agent \(j\)’s column is \(+1\) in games where \(j\) is listed first and \(-1\) where it is listed second, so its derivative is zero exactly when \(j\)’s expected wins equal its actual wins. Given who played whom, the win total is sufficient for strength, as the sum score is for \(\theta\) under the Rasch model; in a round robin, ranking by strength and ranking by wins coincide. “Given who played whom” is where the two part ways.

Recall from Where CTT breaks: the sum score depends on which items were taken. A student scored more on an easy half of a test than on a hard half.

Win percentage ignores who you played; Bradley-Terry adjusts for strength of schedule, as the Rasch model adjusts for the items taken.

Try it: strength of schedule

Six teams in two divisions. Lions, Tigers and Bears are the three strongest; Wolves, Foxes and Hares the three weakest. Each team plays 40 games, a share of them within its own division. The widget uses the expected number of wins, with no chance in it, so that only the schedule differs.

At 0.4 the schedule is a balanced round robin, and win percentage orders the teams correctly. Push the slider up: the Bears, weakest of the strong division, lose more; the Wolves, best of the weak division, win more. From about 0.6 the Wolves have the better record, and at 0.9 they win 61% of their games to the Bears’ 39%. The Bradley-Terry estimates recover the true strengths at every setting, because the model knows whom each team played.

Win percentage assumes, in effect, that every agent faced the same opponents. When the schedule is balanced that holds, and it needs no model; in the simulation below it does as well as Bradley-Terry on a near-random schedule. When the schedule is lopsided, as in leagues with divisions or a comparative-judgment study where some essays meet only weak ones, the assumption fails and the adjustment is what separates the two.

Elo is Bradley-Terry estimated one game at a time

Bradley-Terry fits a season at once. A chess federation has no season: games arrive daily, and a rating is wanted after each. Elo (1978) proposed a rule. Before \(j\) meets \(k\), compute from the current ratings \(\hat p = \text{logistic}(\hat\theta_j - \hat\theta_k)\). After the game, with \(y = 1\) if \(j\) won, \(\tfrac12\) for a draw and 0 for a loss, \[\hat\theta_j \leftarrow \hat\theta_j + K\,(y - \hat p), \qquad \hat\theta_k \leftarrow \hat\theta_k - K\,(y - \hat p).\] A player who wins when expected to gains a little; a player who wins as an underdog gains a lot. One rating goes up by exactly what the other goes down, so the average rating never moves, which sets the origin.

Why this rule? \(y - \hat p\) is the derivative of one game’s Bradley-Terry log likelihood with respect to \(\theta_j\) (see the Go deeper callout below), so Elo takes a step of size \(K\) uphill on the likelihood of the latest game. A small \(K\) averages over many games: stable, but slow to notice a player who improves. A large \(K\) follows change quickly, and follows luck too.

Chess federations report ratings in points rather than logits. In the conversion table of the FIDE Rating Regulations (§8.1.1), an expected score of 0.75, odds of 3 to 1 or \(\log 3 = 1.10\) logits, corresponds to a rating difference of 193 points, so a logit is about 175 points there. FIDE uses \(K\) = 40, 20 or 10 depending on a player’s experience and rating (§8.3.3); \(K = 20\) points is about 0.11 logits.

Try it: Elo following a change

One team plays 200 games against opponents of known strength; its true strength is 0, then jumps at game 100. Elo starts at 0 and sees only results.

With \(K = 0.02\) the rating creeps up and is still well short of the new strength 100 games after the jump. With \(K = 0.6\) it passes 0.8 within three games of the jump but then swings by more than a logit either way. In between (here around \(K = 0.15\)) the mean squared distance is smallest. Set the jump to 0: now there is no change to follow, and the smallest \(K\) does best, because all a larger \(K\) adds is noise.

The same rule, with the item as the opponent, is used to track learners in adaptive practice systems (Pelánek, 2016). Glickman (1999) replaces the fixed \(K\) with a step that depends on how uncertain each rating is, so a rating we know little about moves further. That is the Glicko system.

Assumptions: games are independent given the strengths; game \(g\) between \(j\) and \(k\) follows Bradley-Terry with \(p_g = \text{logistic}(\theta_j - \theta_k)\); \(y_g \in \{0, \tfrac12, 1\}\) (a draw enters as \(\tfrac12\), which is the chess convention rather than a model for draws).

One game’s gradient. The log likelihood of game \(g\) is \(\ell_g = y_g \log p_g + (1 - y_g)\log(1 - p_g)\). As in The likelihood and logistic regression, with \(\partial p_g/\partial\theta_j = p_g(1-p_g)\), \[\frac{\partial \ell_g}{\partial \theta_j} = y_g - p_g, \qquad \frac{\partial \ell_g}{\partial \theta_k} = -(y_g - p_g),\] and every other agent’s derivative is 0. The Elo update adds \(K\) times this gradient, evaluated at the current ratings. The full log likelihood is \(\sum_g \ell_g\), so the Elo update is stochastic gradient ascent on the Bradley-Terry log likelihood: one game at a time, with learning rate \(K\).

What the fixed step does. Robbins & Monro (1951) gave a procedure of this kind, with steps that shrink as observations accumulate, whose estimates converge in probability to the value sought. Elo keeps \(K\) fixed, so its ratings never converge; they keep fluctuating around the strengths, and the fluctuations grow with \(K\). If strengths are fixed, that is a cost. If they drift, it is the point: a fixed step weights recent games more than old ones, so the ratings forget.

Two consequences. The ratings’ sum is fixed (each update adds and subtracts the same amount), so Elo sets the origin by keeping the mean at its starting value. And a home advantage, or any other term in the exponent, gets its own gradient, \(y_g - p_g\) for \(h\); problem 1 asks what updating it would do. Seen this way, Elo is one-step estimation of \(\theta\) with the opponent, or the item, treated as known. Aldous (2017) studies by simulation how accurate such ratings are.

Static or dynamic? Let prediction decide

Bradley-Terry assumes each agent has one strength for the whole period the data cover: plausible for one tournament or one survey, a strong assumption for six basketball seasons of trades and injuries. Elo assumes strength can drift, and pays for it with noise. Which assumption holds is a question the data can answer by prediction: fit on the past, predict games not yet seen, and compare the log likelihood of what happened (Glickman & Jones, 2025). In-sample fit can’t settle it, since a model with a strength per team per season fits better by construction whether or not teams changed. The NBA analysis below makes the comparison.

Paired comparisons as measurement

Thurstone (1927) used paired comparisons to scale things with no natural unit. His law of comparative judgment, in its simplest form (Case V), is Bradley-Terry with the normal curve in place of the logistic. Andrich (1978) showed that with the logistic it gives the same equations for the scale values as the Rasch model, except that the judge’s own level drops out by design, because both options go to the same judge.

What is measured here is the stimuli, not the people; the judges are instruments, and the question they were asked defines the scale. Ask “which is more harmful?” and “which could cause a bigger disaster?” of the same pairs and nothing obliges the two scales to agree. The same holds when teachers judge which of two essays is better (Pollitt, 2012): the question put to the judges plays the part of the rubric.

Simulate

The code below simulates 30 agents with fixed strengths playing 1,200 games on a lopsided schedule (three divisions formed by strength, 80% of games inside the division), fits Bradley-Terry with glm() and compares the estimates, and win percentage, with the truth. The second part lets strengths drift and runs Elo at seven values of \(K\), scoring each by the mean log likelihood of its predictions. It is base R and runs in about a second.

With the seed as written, win percentage correlates 0.65 with the true strengths and Bradley-Terry 0.94. Change within to 0.1, close to a random schedule, and win percentage (0.94) edges out the noisier Bradley-Terry estimates (0.90). Then try drift_sd at 0, 0.05 and 0.15: the best \(K\) grows with the drift, from 0.1 to 0.1–0.2 to 0.4.

Download this code to run it in your own R session. The irw package’s irw_simdata_comp() simulates competition data too.

With real data

Six NBA seasons

Our main example is nba_2012-2018, every regular-season NBA game from 2012-13 to 2017-18, from box scores compiled by Rossotti (2018). Its job is to show strengths that change over time. agent_a is always the home team, so the outcome is keyed as 1 when the home team wins. Download the code; it needs only base R.

Code 
# Competition tables have no items. Each row is one game: agent_a, agent_b, the
# winner ("agent_a", "agent_b" or "draw"), the date (seconds since 1970) and the
# scores. In this table agent_a is always the home team (hometeam == "agent_a").
nba_url <- "https://redivis.com/api/v1/tables/datapages.irw_competitions:v3_2.nba_2012-2018/rows?format=csv"
nba <- read.csv(nba_url)
head(nba)
agent_a score_a agent_b score_b date hometeam winner
SA 102 ATL 100 1385971200 agent_a agent_a
WAS 104 ATL 100 1490166000 agent_a agent_a
WAS 104 ATL 100 1490166000 agent_a agent_a
GS 105 ATL 100 1480320000 agent_a agent_a
NO 105 ATL 100 1391587200 agent_a agent_a
GS 105 ATL 100 1480320000 agent_a agent_a
Code 
nrow(nba)
[1] 9839

The 2016-17 and 2017-18 seasons appear twice, because two source files overlap; we keep one copy of each game.

Code 
# The 2016-17 and 2017-18 seasons appear twice (the source files overlap), so we
# keep one copy of each game. A season runs from October to April.
nba <- nba[!duplicated(nba), ]
nba$day <- as.Date(as.POSIXct(nba$date, origin = "1970-01-01", tz = "UTC"))
nba$season <- as.integer(format(nba$day, "%Y")) - (as.integer(format(nba$day, "%m")) < 8)
nba <- nba[order(nba$date), ]
nba$y <- as.integer(nba$winner == "agent_a")     # 1 if the home team won; no draws
table(season = nba$season)
season
2012 2013 2014 2015 2016 2017 
1229 1230 1230 1230 1230 1230
Code 
teams <- sort(unique(c(nba$agent_a, nba$agent_b)))
c(games = nrow(nba), teams = length(teams), home_win_rate = round(mean(nba$y), 3))
        games         teams home_win_rate 
     7379.000        30.000         0.586

That leaves 7,379 games among 30 teams, 1,230 in each season but one. Home teams won 58.6% of them.

Code 
# Bradley-Terry as a logistic regression: one row per game, a column per team,
# +1 for the home team, -1 for the visitors. The intercept is the home advantage h.
# Only differences are identified, so we drop the first team's column (its strength
# is fixed at 0) and then re-centre the strengths to mean 0.
design <- function(d) {
  X <- matrix(0, nrow(d), length(teams), dimnames = list(NULL, teams))
  X[cbind(seq_len(nrow(d)), match(d$agent_a, teams))] <- 1
  X[cbind(seq_len(nrow(d)), match(d$agent_b, teams))] <- -1
  X
}
fit_bt <- function(d) {
  X <- design(d)
  m <- glm(d$y ~ X[, -1], family = binomial)
  theta <- c(0, coef(m)[-1])
  names(theta) <- teams
  list(h = unname(coef(m)[1]), theta = theta - mean(theta), model = m)
}
bt_all <- fit_bt(nba)
round(c(home_advantage = bt_all$h, se = summary(bt_all$model)$coefficients[1, 2]), 3)
home_advantage             se 
         0.387          0.025
Code 
round(sort(bt_all$theta, decreasing = TRUE)[1:5], 2)
  GS   SA  HOU  LAC  OKC 
1.11 0.98 0.65 0.64 0.63

The home advantage is 0.39 logits (standard error 0.025): between equal teams, the home team wins with probability \(\text{logistic}(0.39) \approx 0.60\). Golden State (GS) and San Antonio (SA) are strongest over the six seasons. Does the schedule adjustment matter?

Code 
# Win percentage against the Bradley-Terry strength, over all six seasons.
win_pct <- sapply(teams, function(t) mean(c(nba$y[nba$agent_a == t], 1 - nba$y[nba$agent_b == t])))
round(c(correlation = cor(win_pct, bt_all$theta),
        rank_correlation = cor(win_pct, bt_all$theta, method = "spearman")), 4)
     correlation rank_correlation 
          0.9953           0.9947
Code 
# At the maximum, each team's expected wins equal its actual wins:
p_hat <- fitted(bt_all$model)
exp_wins <- sapply(teams, function(t) sum(p_hat[nba$agent_a == t]) + sum(1 - p_hat[nba$agent_b == t]))
act_wins <- sapply(teams, function(t) sum(nba$y[nba$agent_a == t]) + sum(1 - nba$y[nba$agent_b == t]))
round(head(cbind(actual = act_wins, expected = exp_wins)), 2)
    actual expected
ATL    257      257
BKN    200      200
BOS    262      262
CHA    217      217
CHI    253      253
CLE    268      268

Hardly: win percentage and strength correlate 0.995. Each team’s expected wins equal its actual wins (the first six are shown), and the NBA schedule is close enough to balanced that the two rankings nearly coincide. Now fit each season on its own.

Code 
# One Bradley-Terry fit per season: who is strongest, and how far apart are teams?
by_season <- t(sapply(sort(unique(nba$season)), function(s) {
  b <- fit_bt(nba[nba$season == s, ])
  c(season = s, home_adv = round(b$h, 2), sd_theta = round(sd(b$theta), 2),
    top_theta = round(max(b$theta), 2))
}))
data.frame(by_season, top_team = sapply(sort(unique(nba$season)), function(s)
  names(which.max(fit_bt(nba[nba$season == s, ])$theta))))
season home_adv sd_theta top_theta top_team
2012 0.57 0.77 1.50 MIA
2013 0.41 0.79 1.36 SA
2014 0.39 0.82 1.69 GS
2015 0.46 0.89 2.30 GS
2016 0.40 0.63 1.59 GS
2017 0.38 0.69 1.43 HOU

The strongest team changes: Miami, then San Antonio, then Golden State for three seasons, then Houston. The spread of strengths (SD 0.63 to 0.89 logits) and the home advantage (0.38 to 0.57) move less. A pooled fit gives each team one strength for all six seasons; how much does that cost when we predict?

Code 
# Fit on 2012-13 to 2016-17, predict every game of 2017-18. The score is the mean
# log likelihood per game (closer to 0 is better; always guessing 0.5 gives -0.693).
train <- nba[nba$season < 2017, ]
test  <- nba[nba$season == 2017, ]
mean_ll <- function(p, y) mean(y * log(p) + (1 - y) * log(1 - p))
bt5 <- fit_bt(train)                                  # all five seasons pooled
bt1 <- fit_bt(nba[nba$season == 2016, ])              # the last season only
p_const <- rep(mean(train$y), nrow(test))             # home advantage, nothing else
p_bt5 <- plogis(bt5$h + design(test) %*% bt5$theta)
p_bt1 <- plogis(bt1$h + design(test) %*% bt1$theta)

# Elo on the logit scale: before each game, predict from the current ratings; after
# it, move both ratings by K times the surprise. Ratings start at 0 and carry over
# from season to season. The home advantage is the five-season estimate.
elo <- function(d, K, h) {
  r <- setNames(rep(0, length(teams)), teams)
  p <- numeric(nrow(d))
  for (g in seq_len(nrow(d))) {
    a <- d$agent_a[g]; b <- d$agent_b[g]
    p[g] <- plogis(h + r[a] - r[b])
    r[a] <- r[a] + K * (d$y[g] - p[g])
    r[b] <- r[b] - K * (d$y[g] - p[g])
  }
  list(p = p, ratings = r)
}
# Choose K on the training seasons only. Elo's predictions are made before each
# game, so they are already out of sample; 2012-13 is a burn-in.
Ks <- c(0.01, 0.02, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.5)
tune <- sapply(Ks, function(K) {
  e <- elo(train, K, bt5$h)
  mean_ll(e$p[train$season > 2012], train$y[train$season > 2012])
})
round(setNames(tune, Ks), 4)
   0.01    0.02    0.05   0.075     0.1    0.15     0.2     0.3     0.5 
-0.6422 -0.6279 -0.6141 -0.6107 -0.6096 -0.6104 -0.6131 -0.6213 -0.6438
Code 
K <- Ks[which.max(tune)]
p_elo <- elo(nba, K, bt5$h)$p[nba$season == 2017]

scores <- c(constant = mean_ll(p_const, test$y), bt_five_seasons = mean_ll(p_bt5, test$y),
            bt_last_season = mean_ll(p_bt1, test$y), elo = mean_ll(p_elo, test$y))
round(scores, 3)
       constant bt_five_seasons  bt_last_season             elo 
         -0.681          -0.697          -0.658          -0.622
Code 
# Is Elo's edge over the pooled fit more than noise? A paired standard error over games.
d_ll <- (test$y * log(p_elo) + (1 - test$y) * log(1 - p_elo)) -
        (test$y * log(p_bt5) + (1 - test$y) * log(1 - p_bt5))
round(c(K = K, difference = mean(d_ll), se = sd(d_ll) / sqrt(length(d_ll))), 3)
         K difference         se 
     0.100      0.075      0.010
Code 
# How much did the teams change? Strengths from the pooled fit against 2017-18 alone.
round(cor(bt5$theta, fit_bt(test)$theta), 2)
[1] 0.33

The step size was chosen on the training seasons only: \(K = 0.1\) logits (about 17 FIDE points) predicted the 2013-14 to 2016-17 games best. On 2017-18, Elo’s mean log likelihood is −0.622 per game, against −0.658 for Bradley-Terry on 2016-17 alone, −0.681 for the constant and −0.697 for the pooled fit. Elo’s edge over the pooled fit is 0.075 per game (standard error 0.010). The last line says why: the pooled strengths correlate only 0.33 with the strengths 2017-18 alone produces. Glickman & Jones (2025) analyse NBA outcomes with dynamic models in more depth.

Where does that leave the choice? It turns on whether strength stayed put over the period the data cover, and held-out games answer that directly: over these six seasons, a rating that forgets predicted 2017-18 better than the pooled fit, by 0.075 per game. Over a short, fixed window, such as one tournament or one survey, the same comparison may come out the other way, and the Bradley-Terry fit brings standard errors and explicit terms for home advantage or order. Fitting on the past and scoring predictions of what came next is a comparison you can run on any competition table, and its result is a sound basis for choosing.

One hundred risks

friedman2019_risk_harm comes from a survey in which Americans compared 100 public risks in pairs (Friedman, 2019; replication data). Its job is to show paired comparisons measuring the things compared. Per the study’s codebook, respondents picked the risk that caused more harm in the past year, so the scale is keyed toward more harm; the picked risk is agent_a.

Code 
# 100 public risks. Each rater saw pairs and picked one: for harm, the risk that
# caused more harm in the past year; for disaster, the one with more potential for
# disaster (Friedman, 2019, codebook). agent_a is always the risk the rater picked.
harm_url <- "https://redivis.com/api/v1/tables/datapages.irw_competitions:v3_2.friedman2019_risk_harm/rows?format=csv"
dis_url  <- "https://redivis.com/api/v1/tables/datapages.irw_competitions:v3_2.friedman2019_risk_disaster/rows?format=csv"
harm <- read.csv(harm_url)
disaster <- read.csv(dis_url)
head(harm, 3)
agent_a agent_b rater winner
Air pollution Air warfare 2125 agent_a
Air pollution Air warfare 889 agent_a
Air pollution Alzheimer’s disease 940 agent_a
Code 
c(judgments = nrow(harm), raters = length(unique(harm$rater)),
  risks = length(unique(c(harm$agent_a, harm$agent_b))), picked_a = mean(harm$winner == "agent_a"),
  median_per_rater = median(table(harm$rater)))
       judgments           raters            risks         picked_a 
           21390             1194              100                1 
median_per_rater 
              15

There are 21,390 judgments by 1,194 raters, 15 for the median rater. Since the chosen risk is always agent_a, the table can’t say which risk was shown first, so there is no order effect to estimate: we fit with no intercept, flipping a coin for which risk is written first so the outcome varies.

Code 
# The table records the chosen risk as agent_a, so we can't estimate an order
# effect. We flip a coin for which risk is written first, so that the outcome
# varies; with no intercept the fit is the same either way.
fit_judgments <- function(d, seed = 1) {
  risks <- sort(unique(c(d$agent_a, d$agent_b)))
  set.seed(seed)
  flip <- rbinom(nrow(d), 1, 0.5) == 1
  first <- ifelse(flip, d$agent_b, d$agent_a)
  second <- ifelse(flip, d$agent_a, d$agent_b)
  y <- as.integer(!flip)                              # 1 if the first risk was picked
  X <- matrix(0, nrow(d), length(risks), dimnames = list(NULL, risks))
  X[cbind(seq_len(nrow(d)), match(first, risks))] <- 1
  X[cbind(seq_len(nrow(d)), match(second, risks))] <- -1
  m <- glm(y ~ 0 + X[, -1], family = binomial)
  theta <- c(0, coef(m))
  names(theta) <- risks
  theta - mean(theta)
}
harm_theta <- fit_judgments(harm)
round(sort(harm_theta, decreasing = TRUE)[1:6], 2)
                Alcohol use                      Cancer 
                       2.05                        1.96 
                    Obesity               Drunk driving 
                       1.85                        1.81 
Prescription drug overdoses                     Smoking 
                       1.75                        1.74
Code 
round(sort(harm_theta)[1:6], 2)
             Asteroid collision               Extraterrestrials 
                          -3.09                           -2.99 
             Volcanic eruptions                    Solar flares 
                          -2.30                           -1.97 
                    Nuclear war High-energy physics experiments 
                          -1.97                           -1.83

Alcohol use (2.05), cancer (1.96) and obesity (1.85) are judged most harmful; asteroid collisions (−3.09) and extraterrestrials (−2.99) least. Nuclear war (−1.97) is near the bottom, consistent with a question about harm in the past year. How repeatable is the scale? Split the raters at random and fit each half.

Code 
# Split the raters at random into two halves and fit each half separately.
set.seed(2)
raters <- unique(harm$rater)
half <- sample(raters, length(raters) / 2)
th1 <- fit_judgments(harm[harm$rater %in% half, ])
th2 <- fit_judgments(harm[!harm$rater %in% half, ])
round(cor(th1, th2), 2)
[1] 0.97

The two halves’ scale values correlate 0.97: each rater made about 15 judgments, but each half pools nearly 600 raters.

Our contrast is friedman2019_risk_disaster, from the same survey and so the same data source: respondents compared the same risks for their potential for disaster.

Code 
# The same 100 risks, judged for their potential for disaster.
dis_theta <- fit_judgments(disaster)
round(cor(harm_theta, dis_theta), 2)
[1] 0.63
Code 
shift <- dis_theta - harm_theta
round(sort(shift, decreasing = TRUE)[1:5], 2)   # rise most on disaster
         Nuclear war   Asteroid collision Biological terrorism 
                2.97                 2.32                 2.01 
 Pandemics / plagues    Nuclear meltdowns 
                1.78                 1.76
Code 
round(sort(shift)[1:5], 2)                      # fall most
Alcohol use     Smoking     Obesity     Falling      Cancer 
      -1.25       -1.19       -1.08       -1.03       -1.00

The two scales correlate only 0.63. Nuclear war rises 2.97 logits, asteroid collisions 2.32 and biological terrorism 2.01; alcohol use falls 1.25, smoking 1.19 and obesity 1.08. The same hundred risks, ordered differently, because the scales answer different questions: a paired-comparison scale measures the stimuli with respect to a criterion, and the criterion belongs in its name. Friedman asked about eight criteria (problem 5 uses a third), and his argument concerns how Americans with different values judge the same risks differently, which one scale averages over (problem 6).

Problems

  1. Derivation. Add a home advantage to the Elo rule, \(\hat p = \text{logistic}(\hat h + \hat\theta_j - \hat\theta_k)\), and use the Go deeper callout to write an update for \(\hat h\) with its own step size. Show that the average rating still never moves. Run it on the NBA data: does \(\hat h\) settle near 0.39?
  2. Real data with a twist. Did the NBA home advantage change? Fit season-specific team strengths with one \(h\), then with a separate \(h\) per season, and compare by likelihood ratio. Would separate \(h\)s have improved the 2017-18 predictions?
  3. Judgment. A district scores teachers by the average Bradley-Terry strength of their students’ essays, from judges’ paired comparisons. What would you need to believe for these scores to be fair to teachers? Name three assumptions, and say which data like the Friedman tables could check.
  4. Design. You have 100 essays and budget for 1,000 comparisons. How would you choose the pairs? Consider when the estimates exist, how much each comparison tells you (information \(p(1-p)\)), and adaptive pairing (Pollitt, 2012). Check your design with the simulation.
  5. Real data. Friedman’s other criteria are IRW tables with the same layout, such as friedman2019_risk_worry. Fit one and compare it with the harm scale. Which risks move most, and why might they? Check the wording in the codebook first.
  6. Challenge (open). The Friedman raters disagree. Split them by a feature you can justify (the replication data hold respondent covariates), fit a mixture of Bradley-Terry models, or give raters random effects. Is there one scale of harm, or several? What would count as evidence either way?

Stuck, or curious about something this lesson doesn’t cover? Ask Claude about this lesson opens a new conversation that already knows what the lesson is about; add your question at the end. Claude shows a caution notice on any prompt that arrives through a link; that’s expected.

Going further

  • Glickman, M. E., & Jones, A. C. (2025). Models and rating systems for head-to-head competition. Annual Review of Statistics and Its Application, 12, 259–282. https://doi.org/10.1146/annurev-statistics-040722-061813. A review of the whole lesson and much beyond it, from Bradley-Terry and ties to Glicko and systems for online games, with an NBA application.
  • Zermelo, E. (1929). Die Berechnung der Turnier-Ergebnisse als ein Maximumproblem der Wahrscheinlichkeitsrechnung. Mathematische Zeitschrift, 29, 436–460. https://doi.org/10.1007/BF01180541. The model, its fitting algorithm and the existence condition, for chess tournaments.
  • Bradley, R. A., & Terry, M. E. (1952). Rank analysis of incomplete block designs: I. The method of paired comparisons. Biometrika, 39(3–4), 324–345. https://doi.org/10.1093/biomet/39.3-4.324. Where the model gets its name.
  • Ford, L. R. (1957). Solution of a ranking problem from binary comparisons. The American Mathematical Monthly, 64(8), 28–33. https://doi.org/10.1080/00029890.1957.11989117. When the maximum likelihood estimates exist.
  • Luce, R. D. (1959). Individual choice behavior: A theoretical analysis. Wiley. https://openlibrary.org/works/OL4319279W. The choice axiom, from which the model follows.
  • Thurstone, L. L. (1927). A law of comparative judgment. Psychological Review, 34(4), 273–286. https://doi.org/10.1037/h0070288. Paired comparisons as psychological measurement.
  • Andrich, D. (1978). Relationships between the Thurstone and Rasch approaches to item scaling. Applied Psychological Measurement, 2(3), 451–462. https://doi.org/10.1177/014662167800200319. The Rasch model and paired comparisons side by side.
  • Davidson, R. R. (1970). On extending the Bradley-Terry model to accommodate ties in paired comparison experiments. Journal of the American Statistical Association, 65(329), 317–328. https://doi.org/10.1080/01621459.1970.10481082; and Rao, P. V., & Kupper, L. L. (1967). Ties in paired-comparison experiments: A generalization of the Bradley-Terry model. Journal of the American Statistical Association, 62(317), 194–204. https://doi.org/10.1080/01621459.1967.10482901. Two ways to model draws rather than score them as half a win; ordered response categories are the subject of Models for polytomous responses.
  • Hunter, D. R. (2004). MM algorithms for generalized Bradley-Terry models. The Annals of Statistics, 32(1). https://doi.org/10.1214/aos/1079120141. The fitting algorithm in the undefeated-team widget, and its extensions.
  • Turner, H., & Firth, D. (2012). Bradley-Terry models in R: The BradleyTerry2 package. Journal of Statistical Software, 48(9). https://doi.org/10.18637/jss.v048.i09. Fitting the model, with covariates, in R.
  • Elo, A. E. (1978). The rating of chessplayers, past and present. Arco. https://openlibrary.org/works/OL6807294W. The rating system, from its author.
  • Glickman, M. E. (1999). Parameter estimation in large dynamic paired comparison experiments. Journal of the Royal Statistical Society, Series C, 48(3), 377–394. https://doi.org/10.1111/1467-9876.00159. The Glicko system, which carries each rating’s uncertainty; Glickman’s Glicko-2 note adds a volatility for each player.
  • Herbrich, R., Minka, T., & Graepel, T. (2006). TrueSkill: A Bayesian skill rating system. Advances in Neural Information Processing Systems, 19. https://proceedings.neurips.cc/paper/2006/hash/f44ee263952e65b3610b8ba51229d1f9-Abstract.html. A Bayesian rating system from online gaming, which also rates players who compete in teams.
  • Glickman, M. E., & Stern, H. S. (1998). A state-space model for National Football League scores. Journal of the American Statistical Association, 93(441), 25–35. https://doi.org/10.1080/01621459.1998.10474084. A full dynamic model, with strengths that evolve over time.
  • Aldous, D. (2017). Elo ratings and the sports model: A neglected topic in applied probability? Statistical Science, 32(4). https://doi.org/10.1214/17-STS628. Open questions about the model, and a simulation study of how accurate Elo ratings are.
  • Hvattum, L. M., & Arntzen, H. (2010). Using ELO ratings for match result prediction in association football. International Journal of Forecasting, 26(3), 460–470. https://doi.org/10.1016/j.ijforecast.2009.10.002. Elo ratings used for prediction.
  • Pelánek, R. (2016). Applications of the Elo rating system in adaptive educational systems. Computers & Education, 98, 169–179. https://doi.org/10.1016/j.compedu.2016.03.017. Elo with the item as the opponent.
  • Pollitt, A. (2012). The method of Adaptive Comparative Judgement. Assessment in Education: Principles, Policy & Practice, 19(3), 281–300. https://doi.org/10.1080/0969594X.2012.665354. Paired comparisons for grading open-ended work.
  • Robbins, H., & Monro, S. (1951). A stochastic approximation method. The Annals of Mathematical Statistics, 22(3), 400–407. https://doi.org/10.1214/aoms/1177729586. Stochastic approximation, for the Go deeper callout.
  • Friedman, J. A. (2019). Priorities for preventive action: Explaining Americans’ divergent reactions to 100 public risks. American Journal of Political Science, 63(1), 181–196. https://doi.org/10.1111/ajps.12400. The source of the risk data, and an argument about why reactions to risks differ across people.
  • From the IRW: the competition data standard, with a worked Elo and Bradley-Terry example on NFL data.

Data sources

Data from the Item Response Warehouse (each table name links to its IRW page); references from IRW biblio.

For instructors

  • Source: new for this course (issue #86); there is no EDUC 252 class behind it. It is built on the IRW competition standard.
  • Timing: the core ideas and widgets take about 50 minutes; the NBA and risk analyses work as a 30-minute lab. The real-data code runs in about ten seconds, most of it the download.
  • Data notes: in nba_2012-2018, agent_a is always the home team (the column is hometeam in this table, homefield in the standard), and two seasons appear twice because of overlapping source files, so the code removes exact duplicates; the 2012-13 season has 1,229 games. In the Friedman tables agent_a is always the chosen risk, which rules out an order effect; the survey module and the position of each pair aren’t in the table. The two Friedman tables are one data source, used as main example and contrast.
  • Held back: draws modelled properly (Davidson; ordinal models), intransitivity and multidimensional strength, rater-level covariates and mixtures (problem 6), Glicko-2 and TrueSkill (Going further). Elo in adaptive learning systems is picked up in Item banks and adaptive testing.
  • Solutions: held back in solutions/competitions.qmd (not published). Problem 6 is open.