# Ridge Regression for Paired Comparisons: A Tractable New Approach, with
Application to Premier League Football

Cristiano Varin  Ca’ Foscari University of Venice, Italy  David Firth
University of Warwick, Coventry, UK

###### Abstract

Paired comparison models, such as Bradley-Terry and Thurstone-Mosteller, are
commonly used to estimate relative strengths of pairwise compared items in
tournament-style datasets. With predictive performance as primary criterion,
we discuss estimation of paired comparison models with a ridge penalty. A new
approach is derived which combines empirical Bayes and composite likelihoods
without any need to re-fit the model, as a convenient alternative to cross-
validation of the ridge tuning parameter. Simulation studies, together with
application to 28 seasons of English Premier League football, demonstrate much
better predictive accuracy of the new approach relative to ordinary maximum
likelihood. While the application of a standard bias-reducing penalty was
found to improve appreciably the performance of maximum likelihood, the ridge
penalty with tuning as developed here yields greater accuracy still.

Keywords: Bradley-Terry model; Composite likelihood; Pairwise empirical Bayes;
Probit regression; Rating; Shrinkage; Sport tournaments; Thurstone-Mosteller
model.

##  1 Introduction

The problem of rating a list of items on the basis of a set of paired
comparisons arises frequently, in a variety of fields including artificial
intelligence (e.g., Chen et al., 2013; Rafailov et al., 2023), bibliometrics
(Stigler, 1994; Varin, Cattelan and Firth, 2016), education (e.g., Wheadon et
al., 2020; Bartholomew and Jones, 2022), forensic science (Thompson et al.,
2018), genetics (e.g., Ma, Wong and Owen, 2012), politics (e.g., Loewen,
Rubenson and Spirling, 2012), psychometrics (e.g., Maydeu-Olivares and
Böckenholt, 2005) and sport (e.g., Glickman and Stern, 2016), to name just a
few. A classical reference for the statistical analysis of paired comparisons
is David (1988). A more recent review is Cattelan (2012), while Aldous (2017)
discusses various aspects of paired comparisons from the point of view of
applied probability.

Although the application areas amenable to paired-comparison modelling are
many, in what follows we adopt the language of sport; this is due to the
football application that motivated this work and which is described in
Section 5. Consider $p$ teams. In paired comparison models, the probability
that team $i$ beats team $j$ in a match between the two teams is

$\text{Pr}(i\;\text{beats}\;j)=F(\mu_{i}-\mu_{j}),$ (1)

where $F(\cdot)$ is a the cumulative distribution function of a zero-symmetric
continuous random variable and $\mu_{i}$ is the strength of team $i$. Popular
choices for $F$ are the logistic distribution which gives the Bradley-Terry
model (Bradley and Terry, 1952) and the normal distribution which corresponds
to the Thurstone-Mosteller model (Thurstone, 1927; Mosteller, 1951). These two
models are quite similar in practice, given the well-known correspondence
between the logit link assumed in the Bradley-Terry model and the probit link
of the Thurstone-Mosteller model (e.g., Agresti, 2002, pages 246–247).

We begin by assuming that the outcome of a match is independent of the order
in which the teams are presented, in other words that there is no home-field
effect; and also that there are no ties. These assumptions will be relaxed
later in Section 3.2.

Maximum likelihood estimation of the vector of strength parameters
$\bm{\mu}=(\mu_{1},\ldots,\mu_{p})^{\top}$ is well known to be problematic, in
part because the estimate of $\mu_{i}$ diverges if team $i$ wins or loses all
of its matches (e.g., Kosmidis and Firth, 2021). Bias-reduced maximum
likelihood (Firth, 1993) overcomes the limitations of maximum likelihood,
guaranteeing finiteness and also better frequentist properties of the strength
estimates (Kosmidis and Firth, 2021). Finite estimates can also be obtained
through other forms of penalized maximum likelihood (Mease, 2003) or Bayesian
methods (Caron and Doucet, 2012).

##  2 Ridge Regression for Paired Comparisons

Denote the outcome of the match between teams $i$ and $j$ as $Y_{ij}=1$ if $i$
wins and $Y_{ij}=-1$ if $j$ wins. In this paper we investigate a natural form
of regularization for estimation of the strength parameters, based on the
ridge-penalized log likelihood

$\ell_{\lambda}(\bm{\mu})=\ell(\bm{\mu})-\frac{\lambda}{2}\sum_{i=1}^{p}\mu_{i}^{2},$
(2)

where $\ell(\bm{\mu})$ is the log-likelihood of the paired comparison model,
_i.e._ ,

$\ell(\bm{\mu})=\sum_{(i,j)\in\mathcal{S}}\log
F\\{y_{ij}(\mu_{i}-\mu_{j})\\}.$ (3)

Here $\mathcal{S}$ denotes the played tournament, _i.e._ , $\mathcal{S}$ is
the set of index pairs $(i,j)$ in which $y_{ij}$ is observed and $\lambda\geq
0$ a scalar tuning parameter. For a fixed value of $\lambda$, computation of
the ridge estimates of the strength parameters is a relatively simple
optimization that can be carried out, for example, with standard quasi-Newton
algorithms. The standard approach to determine a good value for $\lambda$ is
cross-validation. The choice of the splits used in cross-validation requires
care in the case of incomplete tournaments to guarantee that each assessment
set is comparable to the source data. This work is motivated by the estimation
of team strengths during ongoing round-robin tournaments, and in this case
appropriate balance in the assessment sets is obtained by using match-weeks as
cross-validation units.

An alternative to cross-validation is estimation of $\lambda$ using the
empirical Bayes method (Morris, 1983). Ridge regression coincides with
maximum-a-posteriori estimation under the assumption that the strength
parameters are uncorrelated zero-mean normal random variables with precision
$\lambda$,

$\bm{\mu}\sim\mathcal{N}(\bm{0},\lambda^{-1}{\text{I}}_{p}),$ (4)

where ${\text{I}}_{p}$ is the identity matrix of size $p$. The empirical Bayes
recipe estimates $\lambda$ by maximizing the marginal likelihood obtained by
integrating out the team strengths. This leads to a cumbersome integral of
dimension $p$,

$L(\lambda)=\lambda^{p/2}\int_{\mathbb{R}^{p}}\exp\\{\ell(\bm{\mu})\\}\prod_{i=1}^{p}\phi(\lambda^{1/2}\mu_{i})\text{d}{\mu},$
(5)

where $\phi(\cdot)$ is the density of the standard normal random variable.

##  3 The Pairwise Empirical Bayes Method

###  3.1 Estimation of the tuning parameter

We develop here a simple method that avoids the potentially high-dimensional
integration (5) by using composite likelihood methods (Lindsay, 1988; Varin,
Reid and Firth, 2011). The tuning parameter $\lambda$ is estimated by
maximizing the pairwise log-likelihood constructed from all pairs of
correlated matches, i.e., those pairs with one team in common, such as ($i$ vs
$j$) and ($i$ vs $k$). The pairwise log-likelihood for $\lambda$ is

$\ell_{\text{pair}}(\lambda)=\sum_{r,s\in\mathcal{O}}n_{rs}\log
p_{rs}(\lambda),$ (6)

where $\mathcal{O}=\\{-1,+1\\}$ is the set of possible match outcomes (which
will be later extended in Section 3.2 to include ties),
$n_{rs}=\sum_{ijk}1(Y_{ij}=r,Y_{ik}=s)$ is the number of correlated pairs of
matches with outcomes $r$ and $s$ and
$p_{rs}(\lambda)=\mathrm{Pr}\left(Y_{ij}=r,Y_{ik}=s\right)$ is the
corresponding marginal bivariate probability.

In the following we use the Thurstone-Mosteller model, because it is
mathematically more convenient than the Bradley-Terry model. Under the
Thurstone-Mosteller model, the bivariate probabilities $p_{rs}$ have closed-
form expressions by a well-known result sometimes referred to as the Theorem
of Median Dichotomy (Sheppard, 1898; Cox and Wermuth, 2002). The probability
of two successes for the common player $i$ is

$p_{11}(\lambda)=\int_{0}^{\infty}\int_{0}^{\infty}\phi_{2}\left(u,v;\frac{\lambda^{-1}}{1+2\lambda^{-1}}\right)\text{d}u\text{d}v=\frac{(1+\tau)}{4},$

where $\phi_{2}(u,v;\rho)$ is the density function of a bivariate standard
normal variable with correlation $\rho$, and the Kendall rank correlation

$\tau=\frac{2}{\pi}\text{arcsin}\left(\frac{\lambda^{-1}}{1+2\lambda^{-1}}\right)$
(7)

takes a value in the interval $(0,1/3)$. The general expression for $p_{rs}$
depends on whether there is concordance or not in the outcomes of the two
matches:

$p_{rs}(\lambda)=\begin{cases}\displaystyle{\frac{(1+\tau)}{4}},&\text{if}\;r=s,\\\
\displaystyle{\frac{(1-\tau)}{4}},&\text{if}\;r\neq s.\end{cases}$

A closed-form estimate of $\lambda$ is obtained via reparameterization of the
pairwise likelihood in $\tau,$

$\ell_{\text{pair}}(\lambda)=c\log(1+\tau)+d\log(1-\tau),$

where $c$ is the count of concordant correlated pairs (win-win, loss-loss) and
$d$ the count of discordant correlated pairs (win-loss, loss-win). The maximum
pairwise likelihood estimator of $\tau$ is then the familiar Kendall $\tau$
coefficient

$\hat{\tau}=\frac{c-d}{c+d},$

and inversion of equation (7) yields the maximum pairwise likelihood estimate
of $\lambda,$

$\hat{\lambda}=\frac{1-2\sin(\hat{\tau}\pi/2)}{\sin(\hat{\tau}\pi/2)}.$ (8)

In sparse tournaments, where the number of matches per team is small, it is
convenient to consider the classical adjustment to contingency tables that
consists here in adding one imaginary concordant pair and one imaginary
discordant pair to the counts for each of the $p$ teams. The corresponding
small-sample adjusted estimator of $\tau$ is then

$\hat{\tau}=\frac{c-d}{c+d+p},$

which coincides with the maximizer of the penalized pairwise log-likelihood

$\tilde{\ell}_{\text{pair}}(\lambda)={\ell}_{\text{pair}}(\lambda)+p\log(1-\tau^{2}).$
(9)

###  3.2 Home-field effect and ties

The order of presentation of the teams or other items is important in many
paired comparison contexts, notable examples being the home-field effect in
sport competitions or the white player advantage in chess. In paired-
comparison models, such an order effect is typically described through the
inclusion of an intercept term (_e.g._ , Agresti, 2002),

$\Pr(Y_{ij}=1|\bm{\mu})=\Phi(\delta+\mu_{i}-\mu_{j}),$ (10)

where $\Phi(\cdot)$ is the distribution function of a standard normal
variable.

Another common complication in paired comparisons is the presence of ties,
represented by the null outcome $Y_{ij}=0$ whenever $i$ and $j$ compete but
neither of them wins. As discussed in Agresti (1992), ties can be incorporated
in paired comparison models through a cumulative link model. The Thurstone-
Mosteller model with home-field effect corresponds to the latent continuous
process $Z_{ij}=\delta+\mu_{i}-\mu_{j}+\epsilon_{ij},$ where $\epsilon_{ij}$
are uncorrelated standard normal variables. The outcome of the comparison
depends on a threshold parameter $\gamma\geq 0$, so that $j$ beats $i$ if
$Z_{ij}<-\gamma$, $i$ ties with $j$ if $-\gamma\leq Z_{ij}<\gamma$ and $i$
beats $j$ if $Z_{ij}\geq\gamma.$ The corresponding cumulative probabilities
for the three outcomes given the strengths are

$\Pr(Y_{ij}\leq x|\bm{\mu})=\Phi(c_{x}-\mu_{i}+\mu_{j}),\quad
x\in\\{-1,0,1\\},$

for cutpoints $-\infty=c_{-2}<c_{-1}\leq c_{0}<c_{1}=+\infty,$ with
$c_{-1}=-\gamma-\delta$ and $c_{0}=\gamma-\delta$. Ridge estimates of the
strength parameters are still computed by maximizing a penalized likelihood as
in (2), but now with the log-likelihood

$\ell(\bm{\mu})=\sum_{(i,j)\in\mathcal{S}}\log\\{\Phi(c_{y_{ij}}-\mu_{i}+\mu_{j})-\Phi(c_{y_{ij}-1}-\mu_{i}+\mu_{j})\\}.$

As before, we proceed by empirical Bayes estimation of the model parameters
$\gamma$, $\delta$ and $\lambda$ from the marginal distribution of the
matches, and then compute the ridge estimates of the strength parameters.

The marginal probabilities of losing, drawing and winning for the home team
are simply calculated by setting all strength parameters to zero, giving
$\Pr(Y_{ij}=x)=\Phi(c_{x})-\Phi(c_{x-1})$, for $x\in\\{-1,0,1\\}$. From these
marginal probabilities, consistent estimates of $\gamma$ and $\delta$ are
readily obtained as

$\hat{\gamma}=\frac{1}{2}\left\\{\Phi^{-1}(1-\hat{p}_{1})-\Phi^{-1}(\hat{p}_{-1})\right\\},\quad\hat{\delta}=\frac{1}{2}\left\\{\Phi^{-1}(\hat{p}_{1})-\Phi^{-1}(\hat{p}_{-1})\right\\},$
(11)

where $\hat{p}_{-1}=\sum_{ij}1(Y_{ij}=-1)/(n+1)$ and
$\hat{p}_{1}=\sum_{ij}1(Y_{ij}=1)/(n+1)$ are the sample estimates of the
probability of a win away and at home, respectively. In tournaments without
ties, $\hat{\delta}=0$ as expected.

Having estimated $\delta$ and $\gamma$, we proceed with estimation of the
tuning parameter $\lambda$ by maximizing the pairwise log-likelihood (6) where
$p_{rs}(\lambda)$ now denotes the joint probability of the outcomes of two
matches with either the same home team or the same away team. For example, the
probability of a victory and a tie in two matches played at home by $i$ is

$\displaystyle p_{10}(\lambda)$ $\displaystyle=\Pr(Y_{ij}=1,Y_{ik}=0)$
$\displaystyle=\int_{c_{0}}^{c_{1}}\int_{c_{-1}}^{c_{0}}\phi_{2}\left(u,v;\frac{\lambda^{-1}}{1+2\lambda^{-1}}\right)\text{d}u\text{d}v$
$\displaystyle=\int_{\gamma-\delta}^{\infty}\int_{-\gamma-\delta}^{\gamma-\delta}\phi_{2}\left(u,v;\frac{\lambda^{-1}}{1+2\lambda^{-1}}\right)\text{d}u\text{d}v.$

The other bivariate probabilities involved in the pairwise likelihood are
similarly computed. In these probabilities $\gamma$ and $\delta$ are replaced
with their estimates (11) so that the pairwise likelihood depends on $\lambda$
only.

Maximization of the pairwise likelihood for $\lambda$ is conveniently
performed by using the robust algorithm of Brent (1973, ch. 4) on the finite
interval domain of $\tau$ from equation (7), with subsequent back-
transformation as in equation (8). The same small-sample adjustment as before
remains available through the pairwise likelihood penalty of equation (9).

##  4 Simulations

###  4.1 True strengths normally distributed

The performance of the pairwise empirical Bayes method is illustrated by
simulations from the Thurstone-Mosteller model with normally distributed
strength parameters. We considered 120 scenarios given by the combination of:
(1) three values for the number of teams $p\in\\{20,40,60\\}$, (2) eight
values for the tuning parameter $\lambda\in\\{2^{-1},2^{0},\ldots,\allowbreak
2^{6}\\}$ and (3) five sample sizes corresponding to five increasing fractions
of the matches in a double round-robin tournament. More specifically, we
simulated the full double-round-robin tournament and used the first $m$ match-
weeks for training and the remaining $2(p-1)-m$ match-weeks for testing. The
number of training weeks was chosen so as to approximate a specified
proportion of the tournament matches. For example, when we consider 20% of
matches for training and the number of teams is $p=30$, then the training set
is taken to be the first 12 match-weeks because $0.2\times(2\times 29)=11.6$.
The five sample fractions considered for training in our simulations are
$20\%,30\%,40\%,50\%$ and $80\%$ of the matches in the double round-robin
tournament. For each scenario, $1000$ replications were generated and
analysed.

The performance of our pairwise empirical Bayes method is compared here with:
(1) maximum likelihood estimation, (2) bias-reduced maximum likelihood and (3)
ridge estimation with $\lambda$ estimated by cross-validation. Bias-reduced
maximum likelihood estimation is computed through the R (R Core Team, 2024)
package brglm2 (Kosmidis, 2021). Since paired-comparison models are specified
through pairwise differences of strengths, an arbitrary constraint is needed
for identification of the maximum likelihood and bias-reduced maximum
likelihood estimates. In the simulations we used the constraint $\mu_{1}=0$.
Cross-validated ridge estimates were calculated via the glmnet package
(Friedman, Tibshirani and Hastie, 2010; Tay, Narasimhan & Hastie, 2023), using
match-weeks as the assessment unit. Paraphrasing the typical language used in
model validation, we could call such a design a ‘leave one match-week out’
cross-validation.

We simulated matches without ties because we decided to focus on the
estimation of the strength parameters. Estimation of the tie threshold
parameter $\gamma$ would not be problematic, but bias-reduced maximum
likelihood and cross-validated logistic regression are easier to compute in
paired comparison models without ties. In all simulations the home-field
effect parameter was set to $\delta=0.2$, so that approximately 58% of matches
are won by the home team since $\Phi(0.2)\approx 0.58$.

The predictive performance of each method is measured by the classical
logarithmic score LS (Good, 1952), defined as the negative average log-
likelihood for the testing matches computed at the strengths estimated with
the training matches. As reference we consider the naive forecast that
predicts the result of every match using the assumed probability of a home
win, that is $\text{LS}_{\text{naive}}=-0.58\log(0.58)-0.42\log(0.42)=0.68$.
The naive forecast corresponds to predictions for the paired comparison model
with strength parameters all zero, $\mu_{i}=0$ for all $i=1,\ldots,p$. With
this reference, we define the logarithmic skill score as

$\text{LSS}=1-\frac{\text{LS}}{\text{LS}_{\text{naive}}}=1-\frac{\text{LS}}{0.68}.$

The higher the logarithmic skill score, the better the predictive performance.
Better forecasts on average than those obtained by the naive forecast have a
logarithmic skill score above zero.

The simulation results are summarized in Figure 1 which displays the average
of the simulated logarithmic skill scores. We do not report results based on
cross-validation in Figure 1, because they were almost everywhere visually
indistinguishable from those based on the adjusted pairwise empirical Bayes
method. This finding confirms that pairwise empirical Bayes for the Thurstone-
Mosteller model is essentially equivalent to standard cross-validation, but
with the substantial computational advantage that no repeated refitting of the
model is needed.

Figure 1: Logarithmic skill scores computed for incomplete double round-robin
tournaments of different sizes. The columns correspond to various levels of
incompleteness (20%, 30%, 50% or 80% of all possible matches used for training
and the rest of matches for prediction). The rows correspond to different
numbers of teams (20, 40, 60). The team strengths are drawn from a zero-mean
normal distribution with precision $\lambda$. The strengths are estimated by
bias reduced maximum likelihood (BRMLE), maximum likelihood (MLE) and pairwise
empirical Bayes (PEB). The reference method is a naive forecast that uses only
the expected probability of win and corresponds to a zero logarithmic skill
score (grey solid line).

We now comment on the results reported in Figure 1. As expected, as the number
of teams or the fraction of matches used for training increases, all methods
converge to the same predictive performance. Furthermore, as $\lambda$
increases, estimates of team strengths converge to zero and therefore all
methods produce predictions that are on average equivalent to those of the
naive forecast. The predictive performance of maximum likelihood estimation is
completely unsatisfactory when the training set is small and the number of
teams is 20, with its average logarithmic skill scores being worse than the
naive forecast. This result is relevant for us because twenty is the number of
teams in various major football leagues, such as the English Premier League
which motivated this research, which is considered later in Section 5. The
bias-reduced maximum likelihood estimator improves substantially on maximum
likelihood in most settings, but with small training sets and large values of
$\lambda$ it still leads to worse logarithmic skill scores on average than the
naive forecast.

As shown in Figure 1, pairwise empirical Bayes has uniformly the best
predictive performance of all methods considered. Pairwise empirical Bayes was
even found to perform relatively well in a still more extreme setting, with
$p=20$ teams and only 4 matches per team (approximately 10% of the tournament)
used for training; that extreme setting is not displayed in Figure 1 because
maximum likelihood performs so poorly that its logarithmic skill score values
fall too far below zero.

###  4.2 True strengths distributed as Student t

The results just described were obtained under the assumption of normally
distributed strengths, as was used to construct the pairwise empirical Bayes
method. Less benign settings for pairwise empirical Bayes are considered next,
by drawing the strengths from Student’s $t$ distribution with $\nu$ degrees of
freedom. The simulated strengths were scaled so that their precision continues
to be $\lambda$,

$\mu_{i}\sim\sqrt{\frac{\nu-2}{\lambda\nu}}t_{\nu},\quad\nu>2,\;i=1,\ldots,p,$

with $t_{\nu}$ denoting a Student’s $t$-distributed variate with $\nu$ degrees
of freedom. We have considered degrees of freedom $\nu>2$ because otherwise
the variance would be undefined. Drawing the strengths from a $t$ distribution
is of interest because it allows a small number of teams to be distinctly
stronger or weaker than the rest. Figures S1 and S2 of the Supplementary
Material display the average logarithmic scores for the same 120 scenarios
considered previously, but with strengths simulated from Student’s $t$
distributions with $\nu=8$ and $\nu=3$ degrees of freedom. There is no
perceptible difference between the $t_{8}$ and normal distributions, as far as
predictive performance is concerned; and even in the more extreme scenario
with $t_{3}$-distributed strengths, the differences are rather small. While
such simulation evidence is inevitably limited, it does provide some
reassurance that the pairwise empirical Bayes method is robust to excess
kurtosis in the strengths. Departures from normality due instead to skewness
need not be studied, because the Thurstone-Mosteller model is identified
through pairwise differences of the strengths which have a symmetrical
distribution by construction.

##  5 Premier League

Since season 1995–1996, the English Premier League has been a double round
robin between twenty teams. The current format of the tournament therefore
consists of 38 match-weeks. We consider here all match results from the 28
seasons between 1995 and 2023. For each season, we trained the paired
comparison models with the first 10, 15, 20, 25 and 30 match-weeks and
predicted the results of the rest of the tournament. The paired comparison
models were estimated by maximum likelihood estimation using the polr function
from the MASS package (Venables and Ripley, 2002), by bias-reduced maximum
likelihood for the cumulative logistic regression model developed in Kosmidis
(2014) and implemented in the R function bpolr available through the
supplementary materials of that paper, and by our pairwise empirical Bayes
method.

As a reference forecast, we considered the naive predictions calculated from
long-term frequencies of home wins, draws and away wins pooled across all 28
seasons of the Premier League. Those seasons saw 46% home wins, 25% draws and
29% away wins out of a total number of $10,640$ matches. The logarithmic score
for the naive forecast is thus

$\text{LS}_{\text{naive}}=-0.46\log(0.46)-0.25\log(0.25)-0.29\log(0.29)=1.06.$

Figure 2 displays boxplots of the collected logarithmic skill scores: each
boxplot is calculated from the 28 logarithmic skill scores for the various
seasons, in a given match-week.

Figure 2: Boxplot summary of the distributions, across 28 Premier League
seasons, of the logarithmic skill scores computed after 10, 15, 20, 25 and 30
matchweeks. The predictions are computed with strengths estimated using
maximum likelihood (MLE), bias-reduced maximum likelihood (BRMLE) and pairwise
empirical Bayes (PEB). The reference method is a naive forecast that uses only
long-term frequencies of home wins, draws and away wins and corresponds to a
zero logarithmic skill score (grey solid line).

The boxplots show that the predictive performance of pairwise empirical Bayes
is substantially better than that of maximum likelihood estimation, in
particular when the training data are match results from only the first 10 or
15 weeks of the season. At those early stages of the tournament, just as we
saw in the simulation studies of Section 4, maximum likelihood very often
produces strength-parameter estimates that predict future matches worse than
does the naive forecast.

The results here also show that the bias-reduction method of Firth (1993)
improves upon maximum likelihood estimation in terms of prediction: in all 28
seasons and for all sizes of the training data, the logarithmic skill scores
based on bias-reduced maximum likelihood are better than those of standard
maximum likelihood. However, pairwise empirical Bayes typically performs even
better than bias-reduced maximum likelihood. For prediction in the early weeks
of the season, especially, pairwise empirical Bayes out-performs the other
methods by a substantial margin.

We now illustrate the above findings in more detail for the most recent
Premier League season included in our study. Figure 3 displays the logarithmic
skill scores for the 2022–23 season. When only the first ten Premier League
matchweeks are used for training, the predictions calculated from maximum
likelihood estimates are much worse even than the naive forecast, whereas the
pairwise empirical Bayes method substantially out-performs the naive
forecaster. The predictions calculated with bias reduced maximum likelihood
clearly improve upon maximum likelihood, but still not enough to produce
predictions with acceptable precision when calculated from only ten
matchweeks. The performance of maximum likelihood estimation does improve as
the number of matchweeks used for training increases, but in any case pairwise
empirical Bayes remains superior to maximum likelihood with or without bias
correction.

Figure 3: Logarithmic skill scores for Premier League 2022–2023. The
predictions are computed with strengths estimated using maximum likelihood
(MLE), bias-reduced maximum likelihood (BRMLE) and pairwise empirical Bayes
(PEB). The reference method is a naive forecast that uses only long-term
frequencies of home wins, draws and away wins and corresponds to a zero
logarithmic skill score (grey solid line).

Further insight into the relative predictive performance of the three methods
can be obtained from Figure 4, which shows estimated strength parameters for
the twenty Premier League teams after ten weeks of matches and at the end of
the season. The plots show clearly the strong shrinkage effect of the ridge
estimator tuned by the pairwise empirical Bayes method, relative to maximum
likelihood estimation which over-emphasises match results seen in the early
stages of the tournament.

After ten weeks of the 2022–23 season, Arsenal were at the top of the table
with 27 points (nine wins and one defeat), followed by Manchester City with 26
points (eight wins and two draws). With these results, the maximum likelihood
estimated strengths were $3.22$ (Arsenal) and $2.50$ (Manchester City). At the
end of the season the league champions were Manchester City with a five-point
lead over Arsenal, with maximum likelihood estimated strengths $1.35$
(Arsenal) and $1.66$ (Manchester City). While bias-reduced maximum likelihood
estimates after ten matchweeks are sensibly attenuated to $2.55$ (Arsenal) and
$2.02$ (Manchester City), still they are substantially larger then the
corresponding strength estimates at the end of the season.

The ridge-estimated strengths are instead quite stable from the beginning to
the end of the season: the strength estimates obtained by pairwise empirical
Bayes moved only from $0.64$ to $0.62$ for Arsenal and from $0.60$ to $0.75$
for Manchester City, for example, between the 10th and 38th matchweek. Thus
over-confident prediction of future match outcomes, based on only a few match
results early in the season, is avoided. This example illustrates well the
general importance of shrinkage estimation, for prediction with good
precision.

At the other end of the Premier League table the picture is similar. After ten
matchweeks of the 2022–23 season, Leicester City had only 4 league points
(from one win, one draw and eight defeats). Leicester’s corresponding strength
estimate via maximum likelihood at that time was very low, at $-2.99$. By the
end of the season, the maximum likelihood estimate of Leicester’s strength had
improved substantially to $-0.86$. The ridge estimates of strength are again
much more stable, with Leicester moving only from $-0.59$ after ten weeks to
$-0.40$ at the end of the season. Leicester were relegated, still, but not by
the large margin that would have been predicted after ten weeks based on
maximum likelihood (with or without bias reduction). Ridge estimation via
pairwise empirical Bayes, on the other hand, gave a much more realistic
picture of Leicester City’s early-season standing.

Figure 4: Estimated strength parameters of the twenty teams in Premier League
season 2022–2023 after ten matchweeks and at the end of the season (38
matchweeks). Strengths are estimated by the maximum likelihood (MLE), bias-
reduced maximum likelihood (BRMLE) and pairwise empirical Bayes (PEB) methods,
in all cases with their sum constrained to zero. Highlighted are the estimated
strengths of Arsenal (ARS), Manchester City (MNC) and Leicester City (LEI).

##  6 Final Remarks

The main purposes of this paper were to discuss ridge estimation of paired-
comparison models motivated by applications in sport tournaments, and to
develop the simple and effective pairwise empirical Bayes method for tuning
the ridge penalty in the Thurstone-Mosteller model. As expected, the shrinkage
provided by a ridge penalty yields predictions with substantially higher
precision than standard maximum likelihood estimation. Indeed when a paired-
comparison model is trained with relatively few matches, as is the case for
very sparse tournaments or in the early stages of a round robin or similar
tournament, prediction based on maximum likelihood is actually worse than the
naive predictor that uses only the the long-term proportions of wins, draws
and losses.

The focus here was on ‘static’ models, with team strength assumed constant
through each season. A potentially fruitful direction for future work is to
explore penalized estimation for dynamic paired-comparison models, in which
the strength parameters vary during the season. In earlier work Fahrmeir and
Tutz (1994) and Knorr-Held (2000) discuss dynamic models where the strength
parameters evolve according to random walks, whereas Cattelan, Varin and Firth
(2013) describe serial dependence in the strengths through an exponential
smoothing scheme. Other dynamic models for football have used the number of
goals scored by the teams, or their difference, as the response variable,
rather than the match outcome as was used in this work; an example of dynamic
modelling of the number of goals scored is Koopman and Lit (2015).

The simulations of Section 4 and the application to English Premier League
show that the amount of shrinkage needed for good predictive performance when
the training set and the number of teams are both small is appreciably greater
than the shrinkage inherent in the bias-reduced maximum likelihood method of
Firth (1993). This is unsurprising: it is well known that unbiasedness
typically is not the main consideration for predictive performance. Kosmidis
and Firth (2021) suggest that using a suitably scaled-up version of the bias-
reducing penalty might be used instead, to improve predictive performance at
the cost of some additional bias. In sport tournaments such as those
considered in this paper, an appropriate scale factor to use could perhaps be
computed through ‘leave-one-matchweek-out’ cross-validation.

## Funding

The work of both authors was supported by the IRIDE grant 2023 from DAIS, Ca’
Foscari University.

## Supplementary Material

Supplementary material provides the simulation results described in Section
4.2 where strengths are drawn from a Student $t$ distribution. The Premier
League data and R code to reproduce the analysis and graphs of Section 5 are
available in the public repository https://github.com/crisvarin/peb.

## References

  * Agresti (1992) Agresti, A. (1992).  Analysis of ordinal paired comparison data.  Applied Statistics 41, 287–297.
  * Agresti (2002) Agresti, A. (2002).  Categorical Data Analysis. Second Edition.  Wiley, New York.
  * Aldous (2017) Aldous, D. (2017).  Elo ratings and the sports model: A neglected topic in applied probability?  Statistical Science 32, 616–629.
  * Baker and McHale (2017) Baker, R.D. & McHale, I.G. (2017)  An empirical Bayes model for time-varying paired comparisons ratings: Who is the greatest women’s tennis player?  European Journal of Operational Research 258, 328–333.
  * Bartholomew and Jones (2022) Bartholomew, S.R. and Jones, M.D.(2022).  A systematized review of research with adaptive comparative judgment (ACJ) in higher education.  International Journal of Technology and Design Education 32(2), 1159–1190.
  * Bradley and Terry (1952) Bradley, R.A. & Terry, M.E. (1952).  Rank analysis of incomplete block designs I: The method of paired comparisons.  Biometrika 39, 324–345.
  * Brent (1973) Brent, R. (1973).  Algorithms for Minimization without Derivatives.  Prentice-Hall, Englewood Cliffs N.J.
  * Caron and Doucet (2012) Caron, F., & Doucet, A. (2012).  Efficient Bayesian inference for generalized Bradley-Terry models.  Journal of Computational and Graphical Statistics 21, 174–196.
  * Cattelan (2012) Cattelan, M. (2012).  Models for paired comparison data: A review with emphasis on dependent data.  Statistical Science 27, 412–433.
  * Cattelan, Varin and Firth (2013) Cattelan, M., Varin, C. & Firth, D. (2013)  Dynamic Bradley–Terry modelling of sports tournaments.  Journal of the Royal Statistical Society Series C: Applied Statistics 62, 135–150.
  * Chen et al. (2013) Chen, X., Bennett, P.N., Collins-Thompson, K. & Horvitz, E. (2013).  Pairwise ranking aggregation in a crowdsourced setting.  In Proceedings of the sixth ACM international conference on Web search and data mining, 193–202.
  * Cox and Wermuth (2002) Cox, D.R. & Wermuth, N. (2002).  On some models for multivariate binary variables parallel in complexity with the multivariate Gaussian distribution.  Biometrika 89, 462–469.
  * David (1988) David, H. A. (1988).  The Method of Paired Comparisons. 2nd edition.  Griffin’s Statistical Monographs and Courses 41. Griffin, London.
  * Fahrmeir and Tutz (1994) Fahrmeir, L. & Tutz, G. (1994).  Dynamic stochastic models for time-dependent ordered paired comparison systems.  Journal of the American Statistical Association 89, 1438–1449.
  * Firth (1993) Firth, D. (1993).  Bias reduction of maximum likelihood estimates.  Biometrika 80, 27–38.
  * Friedman, Tibshirani and Hastie (2010) Friedman J., Tibshirani R. & Hastie T. (2010).  Regularization Paths for Generalized Linear Models via Coordinate Descent.  Journal of Statistical Software 33(1), 1–22
  * Glickman and Stern (2016) Glickman, M.E. & Stern, H.S. (2017).  Estimating team strength in the NFL.  In Handbook of Statistical Methods and Analyses in Sports, J. Albert, M E. Glickman, T. B. Swartz and R. U. Koning, eds. Boca Raton, USA: Chapman and Hall/CRC, 113–136.
  * Good (1952) Good, I.J. (1952).  Rational decisions.  Journal of the Royal Statistical Society, Series B 14, 107–114.
  * Knorr-Held (2000) Knorr-Held, L. (2000).  Dynamic rating of sports teams.  Journal of the Royal Statistical Society: Series D (The Statistician) 49, 261–276.
  * Koopman and Lit (2015) Koopman, S. J. & R. Lit (2015).  A dynamic bivariate Poisson model for analysing and forecasting match results in the English Premier League.  Journal of the Royal Statistical Society Series A: Statistics in Society 178, 167–186.
  * Kosmidis (2014) Kosmidis, I. (2014).  Improved estimation in cumulative link models.  Journal of the Royal Statistical Society Series B: Statistical Methodology 76(1), 169–196.
  * Kosmidis (2021) Kosmidis, I. (2021).  brglm: Bias reduction in binary-response generalized linear models.  R package version 0.8.2.
  * Kosmidis and Firth (2021) Kosmidis, I. & Firth, D. (2021).  Jeffreys-prior penalty, finiteness and shrinkage in binomial-response generalized linear models.  Biometrika 108, 71–82.
  * Lindsay (1988) Lindsay, B.G. (1988).  Composite likelihood methods.  Contemporary Mathematics 80, 221–239.
  * Loewen, Rubenson and Spirling (2012) Loewen, P.J., Rubenson, D. & Spirling, A.(2012).  Testing the power of arguments in referendums: A Bradley–Terry approach.  Electoral Studies 31(1), 212–221.
  * Ma, Wong and Owen (2012) Ma, L. and Wong, W.H. & Owen, A.B.(2012).  A sparse transmission disequilibrium test for haplotypes based on Bradley–Terry graphs.  Human Heredity 73(1), 52–61.
  * Maydeu-Olivares and Böckenholt (2005) Maydeu-Olivares, A. & Böckenholt, U. (2005).  Structural equation modeling of paired-comparison and ranking data.  Psychological Methods 10, 285–304.
  * Mease (2003) Mease, D. (2003).  A penalized maximum likelihood approach for the ranking of college football teams independent of victory margins.  American Statistician 57, 241–248.
  * Morris (1983) Morris, C.N. (1983).  Parametric empirical Bayes inference: Theory and applications.  Journal of the American Statistical Association 78, 47–65.
  * Mosteller (1951) Mosteller, F. (1951).  Remarks on the method of paired comparisons: I. The least squares solution assuming equal standard deviations and equal correlations.  Psychometrika 16, 3–9.
  * Rafailov et al. (2023) Rafailov, R., Sharma, A., Mitchell, E., Manning, C.D., Ermon, S. & Finn, C.(2023).  Direct preference optimization: Your language model is secretly a reward model.  Advances in Neural Information Processing Systems 36, 53728–53741.
  * R Core Team (2024) R Core Team (2024).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.  https://www.R-project.org/
  * Sheppard (1898) Sheppard, W.F. (1898).  On the geometric treatment of the ‘normal curve’ of statistics with particularly reference to correlation and the theory of errors.  Proceedings of the Royal Society of London 62, 170–173.
  * Stigler (1994) Stigler, S.M. (1994).  Citation patterns in the journals of statistics and probability.  Statistical Science 9, 94–108.
  * Tay, Narasimhan & Hastie (2023) Tay, J.K., Narasimhan, B. & Hastie, T. (2023).  Elastic net regularization paths for all generalized linear models.  Journal of Statistical Software 106(1), 1–31.
  * Thompson et al. (2018) Thompson, W.C., Grady, R.H., Lai, E. & Stern, H.S. (2018).  Perceived strength of forensic scientists’ reporting statements about source conclusions.  Law, Probability and Risk 17(2), 133–155.
  * Thurstone (1927) Thurstone, L.L. (1927).  A law of comparative judgment.  Psychological Review 34, 273–286.
  * Turner and Firth (2012) Turner, H. & Firth, D. (2012).  Bradley-Terry models in R: The BradleyTerry2 package.  Journal of Statistical Software 48 (9), 1–12.
  * Varin, Cattelan and Firth (2016) Varin, C., Cattelan, M. & Firth, D. (2016).  Statistical modelling of citation exchange between Statistics journals.  Journal of the Royal Statistical Society Series A: Statistics in Society 179, 1–63 (with discussion and rejoinder).
  * Varin, Reid and Firth (2011) Varin, C., Reid, N. & Firth, D. (2011).  An overview of composite likelihood methods.  Statistica Sinica 21, 5–42.
  * Venables and Ripley (2002) Venables, W. N. & Ripley, B. D. (2002).  Modern Applied Statistics with S.  Fourth Edition. Springer, New York.
  * Wheadon et al. (2020) Wheadon, C., Barmby, P., Christodoulou, D. & Henderson, B. (2020).  A comparative judgement approach to the large-scale assessment of primary writing in England.  Assessment in Education: Principles, Policy & Practice 1, 46–64.
  * Whiting et al. (2006) Whiting, M.J., Stuart-Fox, D.M., O’Connor, D., Firth, D., Bennett, N.C. & Blomberg, S.P. (2006).  Ultraviolet signals ultra-aggression in a lizard.  Animal Behaviour 72, 353–363.