# Obtain ranking from pairwise comparison with continuous outcome

I have some data of pairwise comparisons between players with a continuous value as outcome, for instance like this:

Player1 Player2 Score
---------------------
A       B       12.38
A       C        4.82
C       B        6.41
C       D        1.39
A       D       -1.05
...    ...       ...


Where a positive score means that Player1 won against Player2 by a meaningful amount of points quantified by "Score", and a negative score meaning that Player2 won. Note that there may be incongruences in the obtained scores, as in the table: A > C, C > D, but A < D.

I have seen there are a lot of techniques to score players when you get a binary outcome (just the information about who wins; see here or here) with the ELO scoring system or the Bradley-Terry Model. However, I believe it would be an unnecessary loss of information to dummy-code each game result as 0-1. Do you know of any way to somehow incorporate that information into a model? Ideally, what I am interested in would be a continuous measure of ability in the game (better) or a ranking of the players.

• Is the score symmetric? If A loses to D by -1.05, does this mean that D beats A by 1.05?
– clp
Dec 28, 2022 at 11:49
• Yes! Do you have something in mind? Dec 29, 2022 at 12:02

## 2 Answers

Note the Elo system is continuous. A rating difference is converted to an expected score percentage as described by Elo in The rating of chess players, chapter 8.94, Development of the Percentage Expectancy Table.

The normal probabilities may be taken directly from the standard tables of the areas under the normal curve when the difference in rating is expressed as a z score. Since the standard deviation σ of individual performances is defined as 200 points, the standard deviation σ' of the differences in performances becomes σ*sqrt(2) or 282.84. The z value of a difference then is D/282.84. This z will then divide the area under the curve into two parts, the larger giving P for the higher rated player and the smaller giving P for the lower rated player.

For example, let D = 160. Then z = 160/282.84 = .566. The table gives .7143 and .2857 as the areas of the two portions under the curve. These probabilities are rounded to two figures in table 2.11.

In chess, the size of these probabilities coincides with the actual scores, 1, ½, 0. An expected percentage score of 50% corresponds to a draw (½). In other scoring systems, the winning probabilities must be converted in score points. When the scoring system has a bounded score, as in matches (2, 1, 0), soccer (3,1,0), team scores, this conversion is straightforward. Additional information is needed to map the scoring system as described above.

A nice introduction to the Elo system and its predecessors can be found in Mark E. Glickman, Introductory Note to 1928.

The FIDE rating system is described in:

#### Update 1 - Current Rating Formula

A player with own rating = 2000, meets five players with ratings 1600, 1700, 1800, 1900, 2100 and scores 12,38, 4,82, 6,41, 1,39, -1,05 points respectively. The rating change is calculated as follows.

Suppose the above score is normally distributed with mean 0 and standard deviation 4. These scores can be translated into a 'normalized' score using the formula (in R), W = pnorm(score, 0, sd=4). Then the expected score We is calculated using the Elo formulas. Finally, when a rating system is conducted on a continuous basis, the rating change is determined using Elo's current rating formula Rn = R0 + K(W-We).

Calculation example Elo's continuous method of evaluation.

              Own rating = 2000
Opponent  Score     W (*)   D     We (**)
1600      12,38     1,00    400   0,92
1700       4,82     0,89    300   0,85
1800       6,41     0,95    200   0,76
1900       1,39     0,64    100   0,64
2100      -1,05     0,40   -100   0,36
Total               3,86     -    3,53  =    0,33
K-factor =   10
Rating change =    3,32
Rn = 2003
(*) pnorm(x, 0, sd=4)
(**) pnorm(D, 0, sd=2000/7)


#### Update 2 - Rating a group of unrated players

Suppose we want to determine the rating based on a single event as in the question, without assuming the rating of these players.

                                      Relative Logistic
Rating  Player  A     B     C     D  rating   dist.
1800     A      x  1.00  0.89  0.40   2017    2035
1800     B   0.00     x  0.05     .   1372    1317
1800     C   0.11  0.95     x  0.64   1878    1897
1800     D   0.60    .   0.36     .   1951    1951


For example, let x be an arbitrary rating vector x = (1800, 1800, 1800, 1800)T. For rating difference D, the scoring probability is set to pnorm(D, 0, 2000 / 7, 1). On this basis, calculate the expected score We(x). Now define the function f(x) = We(x) - W. The roots of this function are the required ratings. The above function forms a system of non-linear equations that can be solved by numerical methods, e.g. in Calc, Excel, Python or R.

For small systems an iterative procedure is practical:

• R0 = (c, ..., c ).
• Ri+1 = Ri - delta.

Where delta is either

• (We - W) / 800. (Newton Raphson) or
• ln(We/W) * 400 / ln(10), calibrated to sum = 0. (Zermelo)

This tournament's solution, the relative ratings, are unique to a constant.

#### Update 3 - Indivisible domains

Suppose player C in the above example beats player B by 1 to 0. Then no non-losing path exists between B and the other players in the tournament. When using a sigmoid curve between zero and 100% as a scoring probability as in the Elo system, no bounded solution exists. Player B is "infinitely" weaker than the other players. A bounded solution exists if and only if the underlying tournament graph is strongly connected.

                                     Relative Logistic Linear
Rating  Player  A     B     C     D  rating   dist.   dist.
1800     A      x  1.00  0.89  0.40   1871    1870    1958
1800     B   0.00     x  0.00     .    -∞      -∞     1504
1800     C   0.11  1.00     x  0.64   1739    1740    1850
1800     D   0.60    .   0.36     .   1790    1790    1818


If the scoring probability P is postulated to be linearly dependent on the rating difference D, for example P = D / 4C + 50%, then all ratings become finite, even if the domain is divisible, as in this example. C = 200 is the Elo class interval. The constant 1/4 is the slope of the logistic expectation at x = 0. Note that these ratings are equivalent to least squares ratings.

One obvious solution (esp. if you only need to calculate this once) is to analyze your data with a mixed effects model, where there is a random player effect. Each match has two players, so you need a multi-membership (always exactly 2 memberships) version of such a model. See e.g. this blog post on how to fit such a model in R, but its even easier to do in something like the brms R package:

library(tidyverse)
library(brms)
library(posterior)
library(patchwork)

set.seed(1234)
data <- tibble(player1 = sample(x=1L:20L, size=500, replace = T),
player2 = sample(x=1L:20L, size=500, replace = T),
score = rnorm(n=500,mean=player1-player2,sd=1),
w1=1L, w2=-1L) %>%
filter(player1 != player2)

p1 <- data %>%
ggplot(aes(x=player1, y=score, col=factor(player2))) +
theme_bw(base_size=18) +
geom_jitter(height=0, width=0.1, alpha=0.5) +
ylab("Observed score difference")

# multi-membership model with two members per group and equal weights
fit1 <- brm(score ~ 0 + (1|mm(player1, player2, weights=cbind(w1, w2), scale=FALSE)),
family = gaussian(),
data = data)

summary(fit1)

p2 <- fit1 %>%
as_draws_df() %>%
pivot_longer(cols=everything(), names_to="parameter", values_to="value") %>%
filter(str_detect(parameter, "r_mmplayer")) %>%
mutate(playerno = as.integer(str_remove_all(str_extract(parameter, "\\[[0-9]+\\,"), "\\[|\\,"))) %>%
ggplot(aes(x=playerno, y=value)) +
theme_bw(base_size=18) +
geom_jitter(width=0.1, height=0, alpha=0.3) +
ylab("Estimated player skill")

p1 + p2


This version fits this as a random effects model, which induces shrinkage on the skill of players with very little data, which may often make sense. If you don't like that, just set a prior that fixes to between player standard deviation to a very large number (prior = prior(class="sd", constant(1000)) or the like).

Of course, one concern could be that the residuals would not be normal, e.g. that there would be more large scores than one would expect, then a Student-t error could be chosen. Or perhaps you worry that in very one-sided matches the better player tries less hard (for which they should perhaps have their estimated skill reduced, but you might also feel that you'd not want to do that). In such a setting, you could introduce a term based on the estimated skill differential that says that the expected score is reduced based on the difference in skill (that might require a more bespoke model, I guess) or something like that. You could also add an intercept to the model, if being player 1 is an advantage vs. being player 2 (e.g. like in chess where it's better to be white rather than black).