Dominance hierarchies from wins

behaviour

social structure

statistics

ecology tutorial

Proportion of wins ranks who an animal fought, not how good it is. David’s score, Bradley-Terry and Elo compared on a win matrix in base R, no packages.

Author

Tidy Ecology

Published

2026-08-09

Watch a group for a season, write down who displaced whom, and you have a square matrix of wins. Turning that into an order looks like counting: whoever wins the highest share of their contests goes on top.

That works when everybody fights everybody equally often, and no group does. Some pairs meet constantly and some almost never, and an animal’s win rate is then a statement about who it happened to face as much as about what it can do. The corrections for this are old, they are short, and none of the standard ones needs a package.

This post builds a group where the answer is known, shows the proportion of wins putting a middling animal on top, and works through what David’s score, Bradley-Terry and Elo each do about it. All four are written out in base R, which is also the only way to see what each one is weighting.

A group with a known order

Ten animals with abilities on a logistic scale, and an encounter design of the sort field data actually has: two animals contest a lot with a few specific partners and once with everybody else.

library(ggplot2)
library(patchwork)

te_paper  <- "#f5f4ee"
te_ink    <- "#16241d"
te_body   <- "#2c3a31"
te_forest <- "#275139"
te_rust   <- "#b5534e"
te_gold   <- "#c9b458"
te_line   <- "#dad9ca"

theme_datasheet <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = te_paper, colour = NA),
          panel.background = element_rect(fill = te_paper, colour = NA),
          panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
          panel.grid.minor = element_blank(),
          text             = element_text(colour = te_body),
          plot.title       = element_text(colour = te_ink, face = "bold"),
          plot.subtitle    = element_text(colour = te_body),
          axis.text        = element_text(colour = te_body))
}

n_ind   <- 10
ids     <- LETTERS[seq_len(n_ind)]
ability <- setNames(seq(2.4, -2.4, length.out = n_ind), ids)
n_many  <- 20
n_few   <- 1

encounters <- matrix(n_few, n_ind, n_ind, dimnames = list(ids, ids))
diag(encounters) <- 0
encounters["G", c("H", "I", "J")] <- n_many
encounters[c("H", "I", "J"), "G"] <- n_many
encounters["A", c("B", "C", "D")] <- n_many
encounters[c("B", "C", "D"), "A"] <- n_many

simulate_wins <- function(seed) {
  set.seed(seed)
  W <- matrix(0, n_ind, n_ind, dimnames = list(ids, ids))
  for (i in seq_len(n_ind - 1)) for (j in (i + 1):n_ind) {
    k <- encounters[i, j]
    if (k == 0) next
    p <- 1 / (1 + exp(-(ability[i] - ability[j])))
    w <- rbinom(1, k, p)
    W[i, j] <- w
    W[j, i] <- k - w
  }
  W
}
wins  <- simulate_wins(31)
total <- wins + t(wins)

A is the strongest animal and J the weakest. G sits seventh of ten, and G’s contests are almost all against the three animals below it. A’s are almost all against the three animals just below it, which are the hardest opponents in the group after A itself.

prop_win  <- rowSums(wins) / rowSums(total)
mean_opp  <- as.numeric(total %*% ability) / rowSums(total)
names(mean_opp) <- ids
n_events  <- sum(wins)
opp_gap   <- mean_opp["A"] - mean_opp["G"]

159 recorded contests. The average opponent A faced has an ability of +1.12 and the average opponent G faced has -1.60, a gap of 2.72 on the same scale the abilities are measured on.

The proportion of wins ranks the opponents

rank_of  <- function(x) rank(-x, ties.method = "min")
prop_top <- names(which.max(prop_win))
prop_rho <- cor(prop_win, ability, method = "spearman")
g_rank_p <- rank_of(prop_win)["G"]

G wins 74.2 per cent of its contests and A wins 71.2 per cent, so the win rate ranks G first, above A. G is seventh of ten. The Spearman correlation between the win rate and the true ability is 0.63, which is not nothing, and is also not an ordering anyone should report.

The mechanism is not observation effort. Effort inflates centrality scores in association networks, which is a separate problem covered elsewhere on the site, and it would show up here as a confounding between win rate and number of contests. What is happening instead is that a win rate has no way to express that one animal’s wins were harder to get than another’s.

Two corrections, and what each one weights

David’s score replaces the raw count with a sum of dyadic proportions, then adds a second round weighted by how strong each defeated opponent was, and subtracts the same construction for losses.

dyadic <- function(W, corrected) {
  tt <- W + t(W)
  P  <- ifelse(tt > 0, W / pmax(tt, 1), 0)
  if (corrected) P <- ifelse(tt > 0, P - (P - 0.5) / (tt + 1), 0)
  P
}
david_score <- function(W, corrected = FALSE) {
  P  <- dyadic(W, corrected)
  w1 <- rowSums(P)
  l1 <- colSums(P)
  setNames(w1 + as.numeric(P %*% w1) - l1 - as.numeric(t(P) %*% l1), ids)
}
ds_raw  <- david_score(wins)
ds_corr <- david_score(wins, corrected = TRUE)

Bradley-Terry goes further and writes a likelihood. Each animal gets a latent ability, the probability that i beats j is the logistic function of the difference, and the whole thing is an ordinary binomial regression on a design matrix of plus and minus ones.

bt_design <- function(W) {
  tt <- W + t(W)
  pr <- which(upper.tri(tt) & tt > 0, arr.ind = TRUE)
  X  <- matrix(0, nrow(pr), n_ind, dimnames = list(NULL, ids))
  for (r in seq_len(nrow(pr))) {
    X[r, pr[r, 1]] <-  1
    X[r, pr[r, 2]] <- -1
  }
  list(X = X[, -n_ind, drop = FALSE],
       y = cbind(W[cbind(pr[, 1], pr[, 2])], W[cbind(pr[, 2], pr[, 1])]))
}
bt_fit <- function(W) {
  d <- bt_design(W)
  suppressWarnings(glm(d$y ~ 0 + d$X, family = binomial))
}
bt_ability <- function(W) {
  a <- c(coef(bt_fit(W)), 0)
  names(a) <- ids
  a - mean(a)
}
bt_hat  <- bt_ability(wins)
bt_rho  <- cor(bt_hat, ability, method = "spearman")
dsr_rho <- cor(ds_raw, ability, method = "spearman")
dsc_rho <- cor(ds_corr, ability, method = "spearman")
g_rank_dsr <- rank_of(ds_raw)["G"]
g_rank_bt  <- rank_of(bt_hat)["G"]
bt_top     <- names(which.max(bt_hat))

Both corrections put G back where it belongs. David’s score moves it from first to 7, Bradley-Terry to 7, and the Bradley-Terry ordering here reproduces the true one exactly, with A on top and a Spearman of 0.99.

rank_tab <- data.frame(
  id    = rep(ids, 4),
  truth = rep(rank_of(ability), 4),
  est   = c(rank_of(prop_win), rank_of(ds_raw), rank_of(ds_corr), rank_of(bt_hat)),
  score = rep(c("proportion of wins", "David's score",
                "David's score, corrected", "Bradley-Terry"), each = n_ind))
rank_tab$score <- factor(rank_tab$score,
  levels = c("proportion of wins", "David's score",
             "David's score, corrected", "Bradley-Terry"))

rank_tab$focal <- ifelse(rank_tab$id == "G", "G", "the rest")

ggplot(rank_tab, aes(truth, est, colour = focal)) +
  geom_abline(slope = 1, intercept = 0, colour = te_line, linewidth = 0.8) +
  geom_point(size = 2.8) +
  scale_colour_manual(values = c(te_rust, te_forest), name = NULL) +
  scale_x_continuous(breaks = seq(1, n_ind, 3)) +
  scale_y_continuous(breaks = seq(1, n_ind, 3)) +
  facet_wrap(~score, nrow = 1) +
  labs(x = "true rank", y = "estimated rank",
       title = "One animal in the wrong place, and where each score puts it") +
  theme_datasheet() +
  theme(legend.position = "bottom",
        strip.text = element_text(colour = te_ink, face = "bold", size = 9))

Four small scatter panels, one per score, each with true rank on the horizontal axis and estimated rank on the vertical and a faint diagonal line. In the proportion of wins panel the highlighted animal, truly seventh, sits at estimated rank one, far below the diagonal. In the three other panels it sits on or beside the diagonal. — Figure 1: Rank under each score against true rank, for one simulated group.

One group is an anecdote, so here is the same design simulated three hundred times.

n_set  <- 300
sweep_scores <- t(vapply(seq_len(n_set), function(s) {
  W  <- simulate_wins(s)
  tt <- W + t(W)
  pw <- rowSums(W) / rowSums(tt)
  dr <- david_score(W)
  dc <- david_score(W, corrected = TRUE)
  bh <- bt_ability(W)
  c(cor(pw, ability, method = "spearman"), cor(dr, ability, method = "spearman"),
    cor(dc, ability, method = "spearman"), cor(bh, ability, method = "spearman"),
    rank_of(pw)["G"], rank_of(dr)["G"], rank_of(dc)["G"], rank_of(bh)["G"],
    names(which.max(pw)) != "A", names(which.max(dr)) != "A",
    names(which.max(dc)) != "A", names(which.max(bh)) != "A")
}, numeric(12)))
score_names <- c("proportion of wins", "David's score",
                 "David's score, corrected", "Bradley-Terry")
summ <- data.frame(
  score    = factor(score_names, levels = score_names),
  rho      = colMeans(sweep_scores[, 1:4]),
  g_rank   = colMeans(sweep_scores[, 5:8]),
  top_miss = colMeans(sweep_scores[, 9:12]),
  g_top3   = colMeans(sweep_scores[, 5:8] <= 3))
top3_cut  <- 3
g_top3_max <- max(summ$g_top3[2:4])
g_rank_mid <- mean(summ$g_rank[2:4])
mc_half    <- 0.5

The win rate places the seventh ranked animal in the top 3 in 98 per cent of the simulated groups, with a mean rank of 2.0. Both versions of David’s score and Bradley-Terry put it in the top 3 in at most 1.7 per cent, with mean ranks around 6.5. On the whole ordering the mean Spearman rises from 0.60 to 0.92 and 0.95.

The interesting column is the last one. Asked only who the top animal is, the raw David’s score is wrong 41 per cent of the time, which is worse than the win rate at 20 per cent. The chance corrected version drops to 15 per cent and Bradley-Terry to 12 per cent.

rho_long <- data.frame(
  rho   = as.numeric(sweep_scores[, 1:4]),
  score = factor(rep(score_names, each = n_set), levels = score_names))

p_rho <- ggplot(rho_long, aes(score, rho, fill = score)) +
  geom_boxplot(outlier.size = 0.6, colour = te_ink, linewidth = 0.35) +
  scale_fill_manual(values = c(te_rust, te_gold, te_ink, te_forest)) +
  labs(x = NULL, y = "Spearman with true ability",
       title = "Whole ordering") +
  theme_datasheet() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 25, hjust = 1, size = 9))

p_top <- ggplot(summ, aes(score, 100 * top_miss, fill = score)) +
  geom_col(width = 0.65) +
  scale_fill_manual(values = c(te_rust, te_gold, te_ink, te_forest)) +
  labs(x = NULL, y = "per cent of groups", title = "Wrong animal named first") +
  theme_datasheet() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 25, hjust = 1, size = 9))

p_rho + p_top + plot_annotation(theme = theme_datasheet())

Two panels. The left panel shows box plots of Spearman correlation with the truth for four scores, lowest and widest for proportion of wins and highest for Bradley-Terry. The right panel is a bar chart of how often each score names the wrong top animal, where the raw David's score bar is the tallest, taller than proportion of wins, and the corrected version and Bradley-Terry are the shortest. — Figure 2: Agreement with the true order, and error at the top, over 300 simulated groups.

The reason is visible in the code. David’s score builds from a matrix of dyadic proportions, and a dyad seen once contributes a proportion of zero or one exactly as firmly as a dyad seen twenty times. In this design most dyads were seen once, so most of the score is carried by the least informative cells. Bradley-Terry never forms a proportion: it counts wins and losses in a likelihood, so a dyad contributes in proportion to how often it was contested. The chance corrected dyadic index, which pulls a proportion towards one half by an amount that depends on the number of encounters, closes most of the gap for the price of one extra term.

When the likelihood has no answer

Bradley-Terry has one structural requirement, and glm() will not tell you when it fails. The maximum likelihood estimate exists only if the digraph of wins is strongly connected: for every way of splitting the group into two parts, somebody in each part must have beaten somebody in the other. Otherwise one group of animals can be pushed arbitrarily far above the other and the likelihood keeps improving.

mat_pow <- function(A, k) {
  B <- diag(nrow(A))
  for (i in seq_len(k)) B <- B %*% A
  B
}
strongly_connected <- function(W) all(mat_pow((W > 0) + diag(nrow(W)), nrow(W) - 1) > 0)

conn <- t(vapply(seq_len(n_set), function(s) {
  W <- simulate_wins(s)
  m <- bt_fit(W)
  c(strongly_connected(W),
    any(rowSums(W) == 0) || any(colSums(W) == 0),
    max(abs(coef(m))))
}, numeric(3)))

share_broken  <- mean(conn[, 1] == 0)
share_extreme <- mean(conn[, 2] == 1)
big_coef      <- median(conn[conn[, 1] == 0, 3])
ok_coef       <- median(conn[conn[, 1] == 1, 3])

Across the three hundred groups the digraph fails to be strongly connected in 17 per cent of them, while an animal that never won or never lost turns up in only 1.3 per cent. Scanning for an undefeated alpha or a hopeless omega therefore catches almost none of the cases, which is the argument for running the digraph test instead. In the affected groups the largest fitted ability has a median absolute value of 24 against 5.3 in the rest, and the standard errors go with them.

undefeated <- wins
undefeated["A", ] <- undefeated["A", ] + undefeated[, "A"]
undefeated[, "A"] <- 0
diag(undefeated) <- 0

m_sep  <- bt_fit(undefeated)
se_sep <- max(summary(m_sep)$coefficients[, 2])
ab_sep <- max(abs(coef(m_sep)))
conn_sep <- strongly_connected(undefeated)

half_add <- undefeated
seen     <- (undefeated + t(undefeated)) > 0
half_add[seen] <- half_add[seen] + 0.5
m_half  <- bt_fit(half_add)
se_half <- max(summary(m_half)$coefficients[, 2])
a_half  <- c(coef(m_half), 0); names(a_half) <- ids
a_half  <- a_half - mean(a_half)

Handing A every one of its contests makes the point loudly. The digraph test returns FALSE, glm() returns coefficients as large as 108 with standard errors up to 37698, and it does so after converging without an error. Adding half a win to each side of every observed dyad brings the largest standard error down to 0.69 and A’s ability to 2.90, at the cost of shrinking every other estimate as well. A penalty or a weakly informative prior is the same repair with a better justification, and the site covers what that looks like for logistic regression in general.

Elo depends on the order you feed it

Elo updates a rating after every contest, which makes it attractive for long observation records and gives it a property the other two do not have: the answer depends on the sequence.

elo_run <- function(seq_i, seq_j, k_factor, init_rating = 1000) {
  r <- setNames(rep(init_rating, n_ind), ids)
  for (t in seq_along(seq_i)) {
    i <- seq_i[t]; j <- seq_j[t]
    p <- 1 / (1 + 10^((r[j] - r[i]) / 400))
    r[i] <- r[i] + k_factor * (1 - p)
    r[j] <- r[j] - k_factor * (1 - p)
  }
  r
}

event_i <- rep(row(wins)[wins > 0], wins[wins > 0])
event_j <- rep(col(wins)[wins > 0], wins[wins > 0])

n_shuffle <- 500
k_grid    <- c(20, 50, 100, 200)
set.seed(4)
elo_sweep <- do.call(rbind, lapply(k_grid, function(kk) {
  out <- t(replicate(n_shuffle, {
    o <- sample(length(event_i))
    r <- elo_run(event_i[o], event_j[o], kk)
    c(cor(r, ability, method = "spearman"), which.max(r), r["A"] - r["B"])
  }))
  data.frame(k_factor = kk, rho = out[, 1],
             top_miss = out[, 2] != 1, gap = out[, 3])
}))
elo_summ <- aggregate(cbind(rho, top_miss, gap) ~ k_factor, elo_sweep, mean)
elo_sd   <- aggregate(cbind(rho, gap) ~ k_factor, elo_sweep, sd)
elo_neg  <- aggregate(gap ~ k_factor, elo_sweep, function(v) mean(v < 0))
k_mid    <- 100
row_mid  <- which(elo_summ$k_factor == k_mid)

The contests are the same contests in every one of the 500 runs; only the order changes. At an update constant of 100 the mean Spearman with the truth is 0.90, which is respectable, and the animal with the highest final rating is not A in 34 per cent of orderings. The rating gap between A and B averages 87 points with a standard deviation of 120, and it comes out negative on 24 per cent of the orderings.

p_k1 <- ggplot(elo_summ, aes(k_factor, rho)) +
  geom_line(colour = te_forest, linewidth = 0.9) +
  geom_point(colour = te_forest, size = 2.6) +
  labs(x = NULL, y = "mean Spearman", title = "A larger update learns faster") +
  theme_datasheet()

p_k2 <- ggplot(elo_summ, aes(k_factor, 100 * top_miss)) +
  geom_line(colour = te_rust, linewidth = 0.9) +
  geom_point(colour = te_rust, size = 2.6) +
  labs(x = "update constant", y = "per cent of orderings",
       title = "And forgets the past faster") +
  theme_datasheet()

p_k1 / p_k2 + plot_annotation(theme = theme_datasheet())

Two panels sharing a horizontal axis of update constant. The upper panel shows mean Spearman correlation with the truth rising from about 0.75 at the smallest constant and levelling off. The lower panel shows the percentage of orderings that name the wrong top animal rising steadily from under 10 per cent to nearly half. — Figure 3: Elo across 500 random orderings of the same contests, at four update constants.

Raising the update constant helps the ordering, because the ratings reach a sensible spread within a short record, and hurts the top of it, because the last few contests then carry most of the weight. There is no setting that removes the dependence, only settings that trade one symptom for the other. If the ratings are meant to describe a stable hierarchy rather than track a changing one, that dependence is pure noise, and the only honest use of a single Elo ordering is one that reports how much it moves under reordering.

What to report

Give the encounter matrix, not only the ranking. The number of contests per dyad is what determines whether any of this is estimable, and it is one table.

Do not rank by proportion of wins. It is the only method here that gets the order wrong for a structural reason rather than a statistical one, and the structure is present in every observational dataset.

Fit Bradley-Terry with glm() and report the standard errors on the abilities. It is three lines to build the design matrix, the model is the same likelihood the specialist packages fit, and the errors tell you which adjacent ranks are not separated.

Run the strong connectivity test before trusting the fit. It costs one matrix power and it is the only reliable warning that the estimates are on their way to infinity.

If Elo is used, report the spread over random reorderings of the same contests, and the update constant. A single Elo ranking without both is one draw from a distribution the method does not show you.

Honest limits

Everything here assumes a single fixed ability per animal for the whole observation period. Real hierarchies change: animals mature, alliances form, an alpha is deposed. Bradley-Terry and David’s score both average over that change and report something in the middle, while Elo tracks it, which is the one setting where Elo’s order dependence is a feature rather than a defect. Nothing in this post tests for a changing order, and the methods that do are a different exercise.

The simulation generates contests from exactly the Bradley-Terry model, so Bradley-Terry is being scored on its home ground. When real dominance is not transitive on a single latent scale, when A beats B, B beats C and C beats A reliably rather than by chance, no one-dimensional score describes it, and all four methods will return an order anyway. Counting circular triads before fitting is the check for that, and this post does not run it.

Unknown dyads are avoided here: every pair meets at least once by construction. In field data many pairs are never seen interacting, and that is where the linearity statistics and the ranking methods diverge most sharply. The chance corrected dyadic index handles rarely seen dyads but not unseen ones, and Bradley-Terry needs the digraph condition, which unseen dyads make much harder to satisfy.

The three hundred replicate summaries carry Monte Carlo error. The share of groups naming the wrong top animal has a standard error of at most 2.9 percentage points, so the ordering of the four methods on that measure is clear but the gap between the corrected David’s score and Bradley-Terry is not established by these runs.

The half win augmentation shown against separation is a demonstration, not a recommendation. It is a crude prior applied uniformly, it shrinks well estimated abilities as hard as unidentified ones, and the amount of shrinkage is a choice nobody has justified. Use it to see that the infinity is an artefact, then fit something with a stated prior.

References

David HA 1987 Biometrika 74(2):432-436 (10.1093/biomet/74.2.432)

Bradley RA, Terry ME 1952 Biometrika 39(3/4):324 (10.2307/2334029)

Sanchez-Tojar A, Schroeder J, Farine DR 2018 Journal of Animal Ecology 87(3):594-608 (10.1111/1365-2656.12776)

--- title: "Dominance hierarchies from wins" description: "Proportion of wins ranks who an animal fought, not how good it is. David's score, Bradley-Terry and Elo compared on a win matrix in base R, no packages." date: "2026-08-09 09:00" categories: [R, behaviour, social structure, statistics, ecology tutorial] image: thumbnail.png image-alt: "A scatter plot on a warm off-white background of estimated rank against true rank under proportion of wins, with one highlighted animal that is truly seventh placed first, far below the diagonal line." --- Watch a group for a season, write down who displaced whom, and you have a square matrix of wins. Turning that into an order looks like counting: whoever wins the highest share of their contests goes on top. That works when everybody fights everybody equally often, and no group does. Some pairs meet constantly and some almost never, and an animal's win rate is then a statement about who it happened to face as much as about what it can do. The corrections for this are old, they are short, and none of the standard ones needs a package. This post builds a group where the answer is known, shows the proportion of wins putting a middling animal on top, and works through what David's score, Bradley-Terry and Elo each do about it. All four are written out in base R, which is also the only way to see what each one is weighting. ## A group with a known order Ten animals with abilities on a logistic scale, and an encounter design of the sort field data actually has: two animals contest a lot with a few specific partners and once with everybody else. ```{r setup} #| message: false #| warning: false library(ggplot2) library(patchwork) te_paper <- "#f5f4ee" te_ink <- "#16241d" te_body <- "#2c3a31" te_forest <- "#275139" te_rust <- "#b5534e" te_gold <- "#c9b458" te_line <- "#dad9ca" theme_datasheet <- function() { theme_minimal(base_size = 12) + theme(plot.background = element_rect(fill = te_paper, colour = NA), panel.background = element_rect(fill = te_paper, colour = NA), panel.grid.major = element_line(colour = te_line, linewidth = 0.3), panel.grid.minor = element_blank(), text = element_text(colour = te_body), plot.title = element_text(colour = te_ink, face = "bold"), plot.subtitle = element_text(colour = te_body), axis.text = element_text(colour = te_body)) } n_ind <- 10 ids <- LETTERS[seq_len(n_ind)] ability <- setNames(seq(2.4, -2.4, length.out = n_ind), ids) n_many <- 20 n_few <- 1 encounters <- matrix(n_few, n_ind, n_ind, dimnames = list(ids, ids)) diag(encounters) <- 0 encounters["G", c("H", "I", "J")] <- n_many encounters[c("H", "I", "J"), "G"] <- n_many encounters["A", c("B", "C", "D")] <- n_many encounters[c("B", "C", "D"), "A"] <- n_many simulate_wins <- function(seed) { set.seed(seed) W <- matrix(0, n_ind, n_ind, dimnames = list(ids, ids)) for (i in seq_len(n_ind - 1)) for (j in (i + 1):n_ind) { k <- encounters[i, j] if (k == 0) next p <- 1 / (1 + exp(-(ability[i] - ability[j]))) w <- rbinom(1, k, p) W[i, j] <- w W[j, i] <- k - w } W } wins <- simulate_wins(31) total <- wins + t(wins) ``` A is the strongest animal and J the weakest. G sits seventh of ten, and G's contests are almost all against the three animals below it. A's are almost all against the three animals just below it, which are the hardest opponents in the group after A itself. ```{r opponents} prop_win <- rowSums(wins) / rowSums(total) mean_opp <- as.numeric(total %*% ability) / rowSums(total) names(mean_opp) <- ids n_events <- sum(wins) opp_gap <- mean_opp["A"] - mean_opp["G"] ``` `r sprintf("%d", n_events)` recorded contests. The average opponent A faced has an ability of `r sprintf("%+.2f", mean_opp["A"])` and the average opponent G faced has `r sprintf("%+.2f", mean_opp["G"])`, a gap of `r sprintf("%.2f", opp_gap)` on the same scale the abilities are measured on. ## The proportion of wins ranks the opponents ```{r naive} rank_of <- function(x) rank(-x, ties.method = "min") prop_top <- names(which.max(prop_win)) prop_rho <- cor(prop_win, ability, method = "spearman") g_rank_p <- rank_of(prop_win)["G"] ``` G wins `r sprintf("%.1f", 100 * prop_win["G"])` per cent of its contests and A wins `r sprintf("%.1f", 100 * prop_win["A"])` per cent, so the win rate ranks `r sprintf("%s", prop_top)` first, above A. G is seventh of ten. The Spearman correlation between the win rate and the true ability is `r sprintf("%.2f", prop_rho)`, which is not nothing, and is also not an ordering anyone should report. The mechanism is not observation effort. Effort inflates centrality scores in association networks, which is a separate problem covered elsewhere on the site, and it would show up here as a confounding between win rate and number of contests. What is happening instead is that a win rate has no way to express that one animal's wins were harder to get than another's. ## Two corrections, and what each one weights David's score replaces the raw count with a sum of dyadic proportions, then adds a second round weighted by how strong each defeated opponent was, and subtracts the same construction for losses. ```{r david} dyadic <- function(W, corrected) { tt <- W + t(W) P <- ifelse(tt > 0, W / pmax(tt, 1), 0) if (corrected) P <- ifelse(tt > 0, P - (P - 0.5) / (tt + 1), 0) P } david_score <- function(W, corrected = FALSE) { P <- dyadic(W, corrected) w1 <- rowSums(P) l1 <- colSums(P) setNames(w1 + as.numeric(P %*% w1) - l1 - as.numeric(t(P) %*% l1), ids) } ds_raw <- david_score(wins) ds_corr <- david_score(wins, corrected = TRUE) ``` Bradley-Terry goes further and writes a likelihood. Each animal gets a latent ability, the probability that `i` beats `j` is the logistic function of the difference, and the whole thing is an ordinary binomial regression on a design matrix of plus and minus ones. ```{r bradleyterry} bt_design <- function(W) { tt <- W + t(W) pr <- which(upper.tri(tt) & tt > 0, arr.ind = TRUE) X <- matrix(0, nrow(pr), n_ind, dimnames = list(NULL, ids)) for (r in seq_len(nrow(pr))) { X[r, pr[r, 1]] <- 1 X[r, pr[r, 2]] <- -1 } list(X = X[, -n_ind, drop = FALSE], y = cbind(W[cbind(pr[, 1], pr[, 2])], W[cbind(pr[, 2], pr[, 1])])) } bt_fit <- function(W) { d <- bt_design(W) suppressWarnings(glm(d$y ~ 0 + d$X, family = binomial)) } bt_ability <- function(W) { a <- c(coef(bt_fit(W)), 0) names(a) <- ids a - mean(a) } bt_hat <- bt_ability(wins) bt_rho <- cor(bt_hat, ability, method = "spearman") dsr_rho <- cor(ds_raw, ability, method = "spearman") dsc_rho <- cor(ds_corr, ability, method = "spearman") g_rank_dsr <- rank_of(ds_raw)["G"] g_rank_bt <- rank_of(bt_hat)["G"] bt_top <- names(which.max(bt_hat)) ``` Both corrections put G back where it belongs. David's score moves it from first to `r sprintf("%d", g_rank_dsr)`, Bradley-Terry to `r sprintf("%d", g_rank_bt)`, and the Bradley-Terry ordering here reproduces the true one exactly, with `r sprintf("%s", bt_top)` on top and a Spearman of `r sprintf("%.2f", bt_rho)`. ```{r fig-ranks} #| fig-cap: "Rank under each score against true rank, for one simulated group." #| fig-alt: "Four small scatter panels, one per score, each with true rank on the horizontal axis and estimated rank on the vertical and a faint diagonal line. In the proportion of wins panel the highlighted animal, truly seventh, sits at estimated rank one, far below the diagonal. In the three other panels it sits on or beside the diagonal." rank_tab <- data.frame( id = rep(ids, 4), truth = rep(rank_of(ability), 4), est = c(rank_of(prop_win), rank_of(ds_raw), rank_of(ds_corr), rank_of(bt_hat)), score = rep(c("proportion of wins", "David's score", "David's score, corrected", "Bradley-Terry"), each = n_ind)) rank_tab$score <- factor(rank_tab$score, levels = c("proportion of wins", "David's score", "David's score, corrected", "Bradley-Terry")) rank_tab$focal <- ifelse(rank_tab$id == "G", "G", "the rest") ggplot(rank_tab, aes(truth, est, colour = focal)) + geom_abline(slope = 1, intercept = 0, colour = te_line, linewidth = 0.8) + geom_point(size = 2.8) + scale_colour_manual(values = c(te_rust, te_forest), name = NULL) + scale_x_continuous(breaks = seq(1, n_ind, 3)) + scale_y_continuous(breaks = seq(1, n_ind, 3)) + facet_wrap(~score, nrow = 1) + labs(x = "true rank", y = "estimated rank", title = "One animal in the wrong place, and where each score puts it") + theme_datasheet() + theme(legend.position = "bottom", strip.text = element_text(colour = te_ink, face = "bold", size = 9)) ``` One group is an anecdote, so here is the same design simulated three hundred times. ```{r repeats} n_set <- 300 sweep_scores <- t(vapply(seq_len(n_set), function(s) { W <- simulate_wins(s) tt <- W + t(W) pw <- rowSums(W) / rowSums(tt) dr <- david_score(W) dc <- david_score(W, corrected = TRUE) bh <- bt_ability(W) c(cor(pw, ability, method = "spearman"), cor(dr, ability, method = "spearman"), cor(dc, ability, method = "spearman"), cor(bh, ability, method = "spearman"), rank_of(pw)["G"], rank_of(dr)["G"], rank_of(dc)["G"], rank_of(bh)["G"], names(which.max(pw)) != "A", names(which.max(dr)) != "A", names(which.max(dc)) != "A", names(which.max(bh)) != "A") }, numeric(12))) score_names <- c("proportion of wins", "David's score", "David's score, corrected", "Bradley-Terry") summ <- data.frame( score = factor(score_names, levels = score_names), rho = colMeans(sweep_scores[, 1:4]), g_rank = colMeans(sweep_scores[, 5:8]), top_miss = colMeans(sweep_scores[, 9:12]), g_top3 = colMeans(sweep_scores[, 5:8] <= 3)) top3_cut <- 3 g_top3_max <- max(summ$g_top3[2:4]) g_rank_mid <- mean(summ$g_rank[2:4]) mc_half <- 0.5 ``` The win rate places the seventh ranked animal in the top `r sprintf("%d", top3_cut)` in `r sprintf("%.0f", 100 * summ$g_top3[1])` per cent of the simulated groups, with a mean rank of `r sprintf("%.1f", summ$g_rank[1])`. Both versions of David's score and Bradley-Terry put it in the top `r sprintf("%d", top3_cut)` in at most `r sprintf("%.1f", 100 * g_top3_max)` per cent, with mean ranks around `r sprintf("%.1f", g_rank_mid)`. On the whole ordering the mean Spearman rises from `r sprintf("%.2f", summ$rho[1])` to `r sprintf("%.2f", summ$rho[2])` and `r sprintf("%.2f", summ$rho[4])`. The interesting column is the last one. Asked only who the top animal is, the raw David's score is wrong `r sprintf("%.0f", 100 * summ$top_miss[2])` per cent of the time, which is worse than the win rate at `r sprintf("%.0f", 100 * summ$top_miss[1])` per cent. The chance corrected version drops to `r sprintf("%.0f", 100 * summ$top_miss[3])` per cent and Bradley-Terry to `r sprintf("%.0f", 100 * summ$top_miss[4])` per cent. ```{r fig-methods} #| fig-cap: "Agreement with the true order, and error at the top, over 300 simulated groups." #| fig-alt: "Two panels. The left panel shows box plots of Spearman correlation with the truth for four scores, lowest and widest for proportion of wins and highest for Bradley-Terry. The right panel is a bar chart of how often each score names the wrong top animal, where the raw David's score bar is the tallest, taller than proportion of wins, and the corrected version and Bradley-Terry are the shortest." rho_long <- data.frame( rho = as.numeric(sweep_scores[, 1:4]), score = factor(rep(score_names, each = n_set), levels = score_names)) p_rho <- ggplot(rho_long, aes(score, rho, fill = score)) + geom_boxplot(outlier.size = 0.6, colour = te_ink, linewidth = 0.35) + scale_fill_manual(values = c(te_rust, te_gold, te_ink, te_forest)) + labs(x = NULL, y = "Spearman with true ability", title = "Whole ordering") + theme_datasheet() + theme(legend.position = "none", axis.text.x = element_text(angle = 25, hjust = 1, size = 9)) p_top <- ggplot(summ, aes(score, 100 * top_miss, fill = score)) + geom_col(width = 0.65) + scale_fill_manual(values = c(te_rust, te_gold, te_ink, te_forest)) + labs(x = NULL, y = "per cent of groups", title = "Wrong animal named first") + theme_datasheet() + theme(legend.position = "none", axis.text.x = element_text(angle = 25, hjust = 1, size = 9)) p_rho + p_top + plot_annotation(theme = theme_datasheet()) ``` The reason is visible in the code. David's score builds from a matrix of dyadic proportions, and a dyad seen once contributes a proportion of zero or one exactly as firmly as a dyad seen twenty times. In this design most dyads were seen once, so most of the score is carried by the least informative cells. Bradley-Terry never forms a proportion: it counts wins and losses in a likelihood, so a dyad contributes in proportion to how often it was contested. The chance corrected dyadic index, which pulls a proportion towards one half by an amount that depends on the number of encounters, closes most of the gap for the price of one extra term. ## When the likelihood has no answer Bradley-Terry has one structural requirement, and `glm()` will not tell you when it fails. The maximum likelihood estimate exists only if the digraph of wins is strongly connected: for every way of splitting the group into two parts, somebody in each part must have beaten somebody in the other. Otherwise one group of animals can be pushed arbitrarily far above the other and the likelihood keeps improving. ```{r connectivity} mat_pow <- function(A, k) { B <- diag(nrow(A)) for (i in seq_len(k)) B <- B %*% A B } strongly_connected <- function(W) all(mat_pow((W > 0) + diag(nrow(W)), nrow(W) - 1) > 0) conn <- t(vapply(seq_len(n_set), function(s) { W <- simulate_wins(s) m <- bt_fit(W) c(strongly_connected(W), any(rowSums(W) == 0) || any(colSums(W) == 0), max(abs(coef(m)))) }, numeric(3))) share_broken <- mean(conn[, 1] == 0) share_extreme <- mean(conn[, 2] == 1) big_coef <- median(conn[conn[, 1] == 0, 3]) ok_coef <- median(conn[conn[, 1] == 1, 3]) ``` Across the three hundred groups the digraph fails to be strongly connected in `r sprintf("%.0f", 100 * share_broken)` per cent of them, while an animal that never won or never lost turns up in only `r sprintf("%.1f", 100 * share_extreme)` per cent. Scanning for an undefeated alpha or a hopeless omega therefore catches almost none of the cases, which is the argument for running the digraph test instead. In the affected groups the largest fitted ability has a median absolute value of `r sprintf("%.0f", big_coef)` against `r sprintf("%.1f", ok_coef)` in the rest, and the standard errors go with them. ```{r separation} undefeated <- wins undefeated["A", ] <- undefeated["A", ] + undefeated[, "A"] undefeated[, "A"] <- 0 diag(undefeated) <- 0 m_sep <- bt_fit(undefeated) se_sep <- max(summary(m_sep)$coefficients[, 2]) ab_sep <- max(abs(coef(m_sep))) conn_sep <- strongly_connected(undefeated) half_add <- undefeated seen <- (undefeated + t(undefeated)) > 0 half_add[seen] <- half_add[seen] + 0.5 m_half <- bt_fit(half_add) se_half <- max(summary(m_half)$coefficients[, 2]) a_half <- c(coef(m_half), 0); names(a_half) <- ids a_half <- a_half - mean(a_half) ``` Handing A every one of its contests makes the point loudly. The digraph test returns `r sprintf("%s", conn_sep)`, `glm()` returns coefficients as large as `r sprintf("%.0f", ab_sep)` with standard errors up to `r sprintf("%.0f", se_sep)`, and it does so after converging without an error. Adding half a win to each side of every observed dyad brings the largest standard error down to `r sprintf("%.2f", se_half)` and A's ability to `r sprintf("%.2f", a_half["A"])`, at the cost of shrinking every other estimate as well. A penalty or a weakly informative prior is the same repair with a better justification, and the site covers what that looks like for logistic regression in general. ## Elo depends on the order you feed it Elo updates a rating after every contest, which makes it attractive for long observation records and gives it a property the other two do not have: the answer depends on the sequence. ```{r elo} elo_run <- function(seq_i, seq_j, k_factor, init_rating = 1000) { r <- setNames(rep(init_rating, n_ind), ids) for (t in seq_along(seq_i)) { i <- seq_i[t]; j <- seq_j[t] p <- 1 / (1 + 10^((r[j] - r[i]) / 400)) r[i] <- r[i] + k_factor * (1 - p) r[j] <- r[j] - k_factor * (1 - p) } r } event_i <- rep(row(wins)[wins > 0], wins[wins > 0]) event_j <- rep(col(wins)[wins > 0], wins[wins > 0]) n_shuffle <- 500 k_grid <- c(20, 50, 100, 200) set.seed(4) elo_sweep <- do.call(rbind, lapply(k_grid, function(kk) { out <- t(replicate(n_shuffle, { o <- sample(length(event_i)) r <- elo_run(event_i[o], event_j[o], kk) c(cor(r, ability, method = "spearman"), which.max(r), r["A"] - r["B"]) })) data.frame(k_factor = kk, rho = out[, 1], top_miss = out[, 2] != 1, gap = out[, 3]) })) elo_summ <- aggregate(cbind(rho, top_miss, gap) ~ k_factor, elo_sweep, mean) elo_sd <- aggregate(cbind(rho, gap) ~ k_factor, elo_sweep, sd) elo_neg <- aggregate(gap ~ k_factor, elo_sweep, function(v) mean(v < 0)) k_mid <- 100 row_mid <- which(elo_summ$k_factor == k_mid) ``` The contests are the same contests in every one of the `r sprintf("%d", n_shuffle)` runs; only the order changes. At an update constant of `r sprintf("%d", k_mid)` the mean Spearman with the truth is `r sprintf("%.2f", elo_summ$rho[row_mid])`, which is respectable, and the animal with the highest final rating is not A in `r sprintf("%.0f", 100 * elo_summ$top_miss[row_mid])` per cent of orderings. The rating gap between A and B averages `r sprintf("%.0f", elo_summ$gap[row_mid])` points with a standard deviation of `r sprintf("%.0f", elo_sd$gap[row_mid])`, and it comes out negative on `r sprintf("%.0f", 100 * elo_neg$gap[row_mid])` per cent of the orderings. ```{r fig-elo} #| fig-cap: "Elo across 500 random orderings of the same contests, at four update constants." #| fig-alt: "Two panels sharing a horizontal axis of update constant. The upper panel shows mean Spearman correlation with the truth rising from about 0.75 at the smallest constant and levelling off. The lower panel shows the percentage of orderings that name the wrong top animal rising steadily from under 10 per cent to nearly half." p_k1 <- ggplot(elo_summ, aes(k_factor, rho)) + geom_line(colour = te_forest, linewidth = 0.9) + geom_point(colour = te_forest, size = 2.6) + labs(x = NULL, y = "mean Spearman", title = "A larger update learns faster") + theme_datasheet() p_k2 <- ggplot(elo_summ, aes(k_factor, 100 * top_miss)) + geom_line(colour = te_rust, linewidth = 0.9) + geom_point(colour = te_rust, size = 2.6) + labs(x = "update constant", y = "per cent of orderings", title = "And forgets the past faster") + theme_datasheet() p_k1 / p_k2 + plot_annotation(theme = theme_datasheet()) ``` Raising the update constant helps the ordering, because the ratings reach a sensible spread within a short record, and hurts the top of it, because the last few contests then carry most of the weight. There is no setting that removes the dependence, only settings that trade one symptom for the other. If the ratings are meant to describe a stable hierarchy rather than track a changing one, that dependence is pure noise, and the only honest use of a single Elo ordering is one that reports how much it moves under reordering. ## What to report Give the encounter matrix, not only the ranking. The number of contests per dyad is what determines whether any of this is estimable, and it is one table. Do not rank by proportion of wins. It is the only method here that gets the order wrong for a structural reason rather than a statistical one, and the structure is present in every observational dataset. Fit Bradley-Terry with `glm()` and report the standard errors on the abilities. It is three lines to build the design matrix, the model is the same likelihood the specialist packages fit, and the errors tell you which adjacent ranks are not separated. Run the strong connectivity test before trusting the fit. It costs one matrix power and it is the only reliable warning that the estimates are on their way to infinity. If Elo is used, report the spread over random reorderings of the same contests, and the update constant. A single Elo ranking without both is one draw from a distribution the method does not show you. ## Honest limits Everything here assumes a single fixed ability per animal for the whole observation period. Real hierarchies change: animals mature, alliances form, an alpha is deposed. Bradley-Terry and David's score both average over that change and report something in the middle, while Elo tracks it, which is the one setting where Elo's order dependence is a feature rather than a defect. Nothing in this post tests for a changing order, and the methods that do are a different exercise. The simulation generates contests from exactly the Bradley-Terry model, so Bradley-Terry is being scored on its home ground. When real dominance is not transitive on a single latent scale, when A beats B, B beats C and C beats A reliably rather than by chance, no one-dimensional score describes it, and all four methods will return an order anyway. Counting circular triads before fitting is the check for that, and this post does not run it. Unknown dyads are avoided here: every pair meets at least once by construction. In field data many pairs are never seen interacting, and that is where the linearity statistics and the ranking methods diverge most sharply. The chance corrected dyadic index handles rarely seen dyads but not unseen ones, and Bradley-Terry needs the digraph condition, which unseen dyads make much harder to satisfy. The three hundred replicate summaries carry Monte Carlo error. The share of groups naming the wrong top animal has a standard error of at most `r sprintf("%.1f", 100 * mc_half / sqrt(n_set))` percentage points, so the ordering of the four methods on that measure is clear but the gap between the corrected David's score and Bradley-Terry is not established by these runs. The half win augmentation shown against separation is a demonstration, not a recommendation. It is a crude prior applied uniformly, it shrinks well estimated abilities as hard as unidentified ones, and the amount of shrinkage is a choice nobody has justified. Use it to see that the infinity is an artefact, then fit something with a stated prior. ## References David HA 1987 Biometrika 74(2):432-436 (10.1093/biomet/74.2.432) Bradley RA, Terry ME 1952 Biometrika 39(3/4):324 (10.2307/2334029) Sanchez-Tojar A, Schroeder J, Farine DR 2018 Journal of Animal Ecology 87(3):594-608 (10.1111/1365-2656.12776) ## Related tutorials - [Network centrality and sampling effort](../network-centrality-sampling-effort/) - [Permutation tests for social networks](../permutation-tests-social-networks/) - [Ranking sites when every estimate is noisy](../ranking-sites-with-uncertainty/) - [Bayesian logistic regression under separation](../bayesian-logistic-regression-separation/)