1 Introduction
1.1 Modelling and predicting competitive sports
Competitive sports refers to any sport that involves two teams or individuals competing against each other to achieve higher scores. Competitive team sports includes some of the most popular and most watched games such as football, basketball and rugby. Such sports are played both in domestic professional leagues such as the National Basketball Association, and international competitions such as the FIFA World Cup. For football alone, there are over one hundred fully professional leagues in 71 countries globally. It is estimated that the Premier League, the top football league in the United Kingdom, attracted a (cumulative) television audience of 4.7 billion viewers in the last season
[47].The outcome of a match is determined by a large number of factors. Just to name a few, they might involve the competitive strength of each individual player in both teams, the smoothness of collaboration between players, and the team’s strategy of playing. Moreover, the composition of any team changes over the years, for example because players leave or join the team. The team composition may also change within the tournament season or even during a match because of injuries or penalties.
Understanding these factors is, by the predictionvalidation nature of the scientific method, closely linked to predicting the outcome of a pairing. By Occam’s razor, the factors which empirically help in prediction are exactly those that one may hypothesize to be relevant for the outcome.
Since keeping track of all relevant factors is unrealistic, of course one cannot expect a certain prediction of a competitive sports outcome. Moreover, it is also unreasonable to believe that all factors can be measured or controlled, hence it is reasonable to assume that unpredictable, or nondeterministic statistical “noise” is involved in the process of generating the outcome (or subsume the unknowns as such noise). A good prediction will, hence, not exactly predict the outcome, but will anticipate the “correct” odds more precisely. The extent to which the outcomes are predictable may hence be considered as a surrogate quantifier of how much the outcome of a match is influenced by “skill” (as surrogated by determinism/prediction), or by “chance”^{1}^{1}1We expressly avoid use of the word “luck” as in vernacular use it often means “chance”, jointly with the belief that it may be influenced by esoterical, magical or otherwise metaphysical means. While in the suggested surrogate use, it may well be that the “chance” component of a model subsumes possible points of influence which simply are not measured or observed in the data, an extremely strong corpus of scientific evidence implies that these will not be metaphysical, only unknown  two qualifiers which are obviously not the same, despite strong human tendencies to believe the contrary. (as surrogated by the noise/unknown factors).
Phenomena which can not be specified deterministically are in fact very common in nature. Statistics and probability theory provide ways to make inference under randomness. Therefore, modelling and predicting the results of competitive team sports naturally falls into the area of statistics and machine learning. Moreover, any interpretable predictive model yields a possible explanation of what constitutes factors influencing the outcome.
1.2 History of competitive sports modelling
Research of modeling competitive sports has a long history. In its early days, research was often closely related to sports betting or player/team ranking [22, 26]. The two most influential approaches are due to Bradley and Terry [3] and Élő [15]. The BradleyTerry and Élő models allow estimation of player rating; the Élő system additionally contains algorithmic heuristics to easily update a player’s rank, which have been in use for official chess rankings since the 1960s. The Élő system is also designed to predict the odds of a player winning or losing to the opponent. In contemporary practice, BradleyTerry and Élő type models are broadly used in modelling of sports outcomes and ranking of players, and it has been noted that they are very close mathematically.
In more recent days, relatively diverse modelling approaches originating from the Bayesian statistical framework [37, 13, 20], and also some inspired by machine learning principles [36, 23, 43] have been applied for modelling competitive sports. These models are more expressive and remove some of the BradleyTerry and Élő models’ limitations, though usually at the price of interpretability, computational efficiency, or both.
A more extensive literature overview on existing approaches will be given later in Section 3, as literature spans multiple communities and, in our opinion, a prior exposition of the technical setting and simultaneous straightening of thoughts benefits the understanding and allows us to give proper credit and context for the widely different ideas employed in competitive sports modelling.
1.3 Aim of competitive sports modelling
In literature, the study of competitive team sports may be seen to lie between two primary goals. The first goal is to design models that make good predictions for future match outcome. The second goal is to understand the key factors that influence the match outcome, mostly through retrospective analysis [45, 50]. As explained above, these two aspects are intrinsically connected, and in our view they are the two facets of a single problem: on one hand, proposed influential factors are only scientifically valid if confirmed by falsifiable experiments such as predictions on future matches. If the predictive performance does not increase when information about such factors enters the model, one should conclude by Occam’s razor that these factors are actually irrelevant^{2}^{2}2… to distinguish/characterize the observations, which in some cases may plausibly pertain to restrictions in set of observations, rather than to causative relevance. Hypothetical example: age of football players may be identified as unimportant for the outcome  which may plausibly be due to the fact that the data contained no players of ages 5 or 80, say, as opposed to player age being unimportant in general. Rephrased, it is only unimportant for cases that are plausible to be found in the data set in the first place.. On the other hand, it is plausible to assume that predictions are improved by making use of relevant factors (also known as “features”) become available, for example because they are capable of explaining unmodelled random effects (noise). In light of this, the main problem considered in this work is the (validable and falsifiable) prediction problem, which in machine learning terminology is also known as the supervised learning task.
1.4 Main questions and challenges in competitive sports outcomes prediction
Given the above discussion, the major challenges may be stated as follows:
On the methodological side, what are suitable models for competitive sports outcomes? Current models are not at the same time interpretable, easily computable, allow to use feature information on the teams/players, and allow to predict scores or ternary outcomes. It is an open question how to achieve this in the best way, and this manuscript attempts to highlight a possible path.
The main technical difficulty lies in the fact that offshelf methods do not apply
due to the structured nature of the data:
unlike in individual sports such as running and swimming where the outcome
depends only on the given team, and where the prediction task may
be dealt with classical statistics and machine learning technology
(see [2] for a discussion of this in the context of running),
in competitive team sports the outcome may
be determined by potentially complex interactions between two opposing teams.
In particular, the performance of any team is not measured directly using a simple metric,
but only in relation to the opposing team’s performance.
On the side of domain applications, which in this manuscript is premier league football, it is of great interest to determine the relevant factors determining the outcome, the best way to predict, and which ranking systems are fair and appropriate.
All these questions are related to predictive modelling, as well as the availability of suitable amounts of quality data. Unfortunately, the scarcity of features available in systematic presentation places a hurdle to academic research in competitive team sports, especially when it comes to assessing important factors such as team member characteristics, or strategic considerations during the match.
Moreover, closely linked is also the question to which extent the outcomes are determined by “chance” as opposed to “skill”. Since if on one hypothetical extreme, results would prove to be completely unpredictable, there would be no empirical evidence to distinguish the matches from a game of chance such as flipping a coin. On the other hand, importance of a measurement for predicting would strongly suggest its importance for winning (or losing), though without an experiment not necessarily a causative link.
We attempt to address these questions in the case of premier league football within the confines of readily available data.
1.5 Main contributions
Our main contributions in this manuscript are the following:

We give what we believe to be the first comprehensive literature review of stateofart competitive sports modelling that comprises the multiple communities (BradleyTerry models, Élő type models, Bayesian models, machine learning) in which research so far has been conducted mostly separately.

We present a unified BradleyTerryÉlő model which combines the statistical rigour of the BradleyTerry models with fitting and update strategies similar to that found in the Élő system. Mathematically only a small step, this joint view is essential in a predictive/supervised setting as it allows efficient training and application in an online learning situation. Practically, this step solves some problems of the Élő system (including ranking initialization and choice of Kfactor), and establishes close relations to logistic regression, lowrank matrix completion, and neural networks.

This unified view on BradleyTerryÉlő allows us to introduce classes of joint extensions, the structured logodds models, which unites desirable properties of the extensions found in the disjoint communities: probabilistic prediction of scores and wins/draws/losses, batch/epoch and online learning, as well as the possibility to incorporate features in the prediction, without having to sacrifice structural parsimony of the BradleyTerry models, or simplicity and computational efficiency of Élő’s original approach.

We validate the practical usefulness of the structured logodds models in synthetic experiments and in answering domain questions on English Premier League data, most prominently on the importance of features, fairness of the ranking, as well as on the “chance”“skill” divide.
1.6 Manuscript structure
Section 2 gives an overview of the mathematical setting in competitive sports prediction. Building on the technical context, Section 3 presents a more extensive review of the literature related to the prediction problem of competitive sports, and introduces a joint view on BradleyTerry and Élő type models. Section 4 introduces the structured logodds models, which are validated in empirical experiments in Section 5. Our results and possible future directions for research are discussed in section 6.
Authors’ contributions
This manuscript is based on ZQ’s MSc thesis, submitted September 2016 at University College London, written under supervision of FK. FK provided the ideas of reinterpretation and possible extensions of the Élő model. Literature overview is jointly due to ZQ an FQ, and in parts follows some very helpful pointers by I. Kosmidis (see below). Novel technical ideas in Sections 4.2 to 4.4, and experiments (setup and implementation) are mostly due to ZQ.
The present manuscript is a substantial reworking of the thesis manuscript, jointly done by FK and ZQ.
Acknowledgements
We are thankful to Ioannis Kosmidis for comments on an earlier form of the manuscript, for pointing out some earlier occurrences of ideas presented in it but not given proper credit, as well as relevant literature in the “BradleyTerry” branch.
2 The MathematicalStatistical Setting
This section formulates the prediction task in competitive sports and fixes notation, considering as an instance of supervised learning with several nonstandard structural aspects being of relevance.
2.1 Supervised prediction of competitive outcomes
We introduce the mathematical setting for outcome prediction in competitive team sports. As outlined in the introductory Section 1.1, three crucial features need to be taken into account in this setting:

The outcome of a pairing cannot be exactly predicted prior to the game, even with perfect knowledge of all determinates. Hence it is preferable to predict a probabilistic estimate for all possible match outcomes (win/draw/loss) rather than deterministically choosing one of them.

In a pairing, two teams play against each other, one as a home team and the other as the away or guest team. Not all pairs may play against each other, while others may play multiple times. As a mathematically prototypical (though inaccurate) subcase one may consider all pairs playing exactly once, which gives the observations an implicit matrix structure (row = home team, column = away team). Outcome labels and features crucially depend on the teams constituting the pairing.

Pairings take place over time, and the expected outcomes are plausibly expected to change with (possibly hidden) characteristics of the teams. Hence we will model the temporal dependence explicitly to be able to take it into account when building and checking predictive strategies.
2.1.1 The Generative Model.
Following the above discussion, we will fix a generative model as follows: as in the standard supervised learning setting, we will consider a generative joint random variable
taking values in , where is the set of features (or covariates, independent variables) for each pairing, while is the set of labels (or outcome variables, dependent variables).In our setting, we will consider only the cases and , in which case an observation from is a socalled match outcome, as well as the case , in which case an observation is a socalled final score (in which case, by convention, the first component of is of the home team), or the case of score differences where (in which case, by convention, a positive number is in favour of the home team). From the official rule set of a game (such as football), the match outcome is uniquely determined by a score or score difference. As all the above sets are discrete, predicting will amount to supervised classification (the score difference problem may be phrased as a regression problem, but we will abstain from doing so for technical reasons that become apparent later).
The random variable and its domain shall include information on the teams playing, as well as on the time of the match.
We will suppose there is a set of teams, and for we will denote by the random variable conditioned on the knowledge that is the home team, and is the away team. Note that information in can include any knowledge on either single team or , but also information corresponding uniquely to the pairing .
We will assume that there are teams, which means that the and may be arranged in matrices each.
Further there will be a set of time points at which matches are observed. For we will denote by or an additional conditioning that the outcome is observed at time point .
Note that the indexing and formally amounts to a double conditioning and could be written as and , where are random variables denoting the home team, the away team, and the time of the pairing. Though we do believe that the index/bracket notation is easier to carry through and to follow (including an explicit mirroring of the the “matrix structure”) than the conditional or “graphical models” type notation, which is our main reason for adopting the former and not the latter.
2.1.2 The Observation Model.
By construction, the generative random variable contains all information on having any pairing playing at any time, However, observations in practice will concern two teams playing at a certain time, hence observations in practice will only include independent samples of for some , and never full observations of which can be interpreted as a latent variable.
Note that the observations can be, inprinciple, correlated (or unconditionally dependent) if the pairing or the time is not made explicit (by conditioning which is implicit in the indices ).
An important aspect of our observation model will be that whenever a value of or is observed, it will always come together with the information of the playing teams and the time at which it was observed. This fact will be implicitly made use of in description of algorithms and validation methodology. (formally this could be achieved by explicitly exhibiting/adding as a Cartesian factor of the sampling domains or which we will not do for reasons of clarity and readability)
Two independent batches of data will be observed in the exposition. We will consider:
where and are i.i.d. samples from .
Note that unfortunately (from a notational perspective), one cannot omit the superscripts as in when defining the samples, since the figurative “dies” should be cast anew for each pairing taking place. In particular, if all games would consist of a single pair of teams playing where the results are independent of time, they would all be the same (and not only identically distributed) without the superindex, i.e., without distinguishing different games as different samples from .
2.1.3 The Learning Task.
As set out in the beginning, the main task we will be concerned with is predicting future outcomes given past outcomes and features, observed from the process above. In this work, the features will be assumed to change over time slowly. It is not our primary goal to identify the hidden features in , as they are never observed and hence not accessible as ground truth which can validate our models. However, these will be of secondary interest and considered empirically validated by a wellpredicting model.
More precisely, we will describe methodology for learning and validating predictive models of the type
where
is the set of (discrete probability) distributions on
. That is, given a pairing and a time point at which the teams and play, and information of type , make a probabilistic prediction of the outcome.Most algorithms we discuss will not use added information in , hence will be of type . Some will disregard the time in . Indeed, the latter algorithms are to be considered scientific baselines above which any algorithm using information in and/or has to improve.
The models above will be learnt on a training set , and validated on an independent test set as defined above. In this scenario, will be a random variable which may implicitly depend on but will be independent of . The learning strategy  which is depending on  may take any form and is considered in a full blackbox sense. In the exposition, it will in fact take the form of various parametric and nonparametric prediction algorithms.
The goodness of such an will be evaluated by a loss which compares a probabilistic prediction to the true observation. The best will have a small expected generalization loss
at any future time point and for any pairing . Under mild assumptions, we will argue below that this quantity is estimable from and only mildly dependent on .
Though a good form for is not apriori clear. Also, it is unclear under which assumptions is estimable, due do the conditioning on in the training set. These special aspects of the competitive sports prediction settings will be addressed in the subsequent sections.
2.2 Losses for probablistic classification
In order to evaluate different models, we need a criterion to measure the goodness of probabilistic predictions. The most common error metric used in supervised classification problems is the prediction accuracy. However, the accuracy is often insensitive to probabilistic predictions.
For example, on a certain test case model A predicts a win probability of 60%, while model B predicts a win probability of 95%. If the actual outcome is not win, both models are wrong. In terms of prediction accuracy (or any other nonprobabilistic metric), they are equally wrong because both of them made one mistake. However, model B should be considered better than model A since it predicted the “true” outcome with higher accuracy.
Similarly, if a large number of outcomes of a fair coin toss have been observed as training data, a model that predicts 50% percent for both outcomes on any test data point should be considered more accurate than a model that predicts 100% percent for either outcome 50% of the time.
There exists two commonly used criteria that take into account the probabilistic nature of predictions which we adopt. The first one is the Brier score (Equation 1 below) and the second is the logloss or loglikelihood loss (Equation 2 below). Both losses compare a distribution to an observation, hence mathematically have the signature of a function . By (very slight) abuse of notation, we will identify distributions on (discrete) with its probability mass function; for a distribution , for we write for mass on the observation (= the probability to observe in a random experiment following ).
With this convention, logloss and Brier loss are defined as follows:
(1)  
(2) 
The logloss and the Brier loss functions have the following properties:

the Brier Score is only defined on with an addition/subtraction and a norm defined. This is not necessarily the case in our setting where it may be that . In literature, this is often identified with , though this identification is arbitrary, and the Brier score may change depending on which numbers are used.
On the other hand, the logloss is defined for any and remains unchanged under any renaming or renumbering of a discrete .

For a joint random variable taking values in , it can be shown that the expected losses are minimized by the “correct” prediction .
The two loss functions usually are introduced as empirical losses on a test set , i.e.,
The empirical logloss is the (negative log)likelihood of the test predictions.
The empirical Brier loss, usually called the “Brier score”, is a straightforward translation of the mean squared error used in regression problems to the classification setting, as the expected mean squared error of predicted confidence scores. However, in certain cases, the Brier score is hard to interpret and may behave in unintuitive ways [27], which may partly be seen as a phenomenon caused by abovementioned lack of invariance under class relabelling.
Given this and the interpretability of the empirical logloss as a likelihood, we will use the logloss as principal evaluation metric in the competitive outcome prediction setting.
2.3 Learning with structured and sequential data
The dependency of the observed data on pairing and time makes the prediction task at hand nonstandard. We outline the major consequences for learning and model validation, as well as the implicit assumptions which allow us to tackle these. We will do this separately for the pairing and the temporal structure, as these behave slightly differently.
2.3.1 Conditioning on the pairing
Match outcomes are observed for given pairings , that is, each featurelabelpair will be of form , where as above the subscripts denote conditioning on the pairing. Multiple pairings may be observed in the training set, but not all; some pairings may never be observed.
This has consequences for both learning and validating models.
For model learning, it needs to be made sure that the pairings to be predicted can be predicted from the pairings observed. With other words, the label in the test set that we want to predict is (in a practically substantial way) dependent on the training set . Note that smart models will be able to predict the outcome of a pairing even if it has not been observed before, and even if it has, it will use information from other pairings to improve its predictions
For various parametric models, “predictability” can be related to completability of a data matrix with
as entries. In section 4, we will relate Élő type models to lowrank matrix completion algorithms; completion can be understood as lowrank completion, hence predictability corresponds to completability. Though, exactly working completability out is not the main is not the primary aim of this manuscript, and for our data of interest, the English Premier League, all pairings are observed in any given year, so completability is not an issue. Hence we refer to [33] for a study of lowrank matrix completability. General parametric models may be treated along similar lines.For modelagnostic model validation, it should hold that the expected generalization loss
can be wellestimated by empirical estimation on the test data. For league level team sports data sets, this can be achieved by having multiple years of data available. Since even if not all pairings are observed, usually the set of pairings which is observed is (almost) the same in each year, hence the pairings will be similar in the training and test set if whole years (or halfseasons) are included. Further we will consider an average over all observed pairings, i.e., we will compute the empirical loss on the training set as
By the above argument, the set of all observed pairings in any given year is plausibly modelled as similar, hence it is plausible to conclude that this empirical loss estimates some expected generalization loss
where (possibly dependent) are random variables that select teams which are paired.
Note that this type of aggregate evaluation does not exclude the possibility that predictions for single teams (e.g., newcomers or after restructuring) may be inaccurate, but only that the “average” prediction is good. Further, the assumption itself may be violated if the whole league changes between training and test set.
2.3.2 Conditioning on time
As a second complication, match outcome data is gathered through time. The data
set might display temporal structure and correlation with time. Again, this has consequences for learning and validating the models.
For model learning, models should be able to intrinsically take into account the temporal structure  though as a baseline, timeagnostic models should be tried. A common approach for statistical models is to assume a temporal structure in the latent variables that determine a team’s strength. A different and somewhat adhoc approach proposed by Dixon and Coles [13] is to assign lower weights to earlier observations and estimate parameter by maximizing the weighted loglikelihood function. For machine learning models, the temporal structure is often encoded with handcrafted features.
Similarly, one may opt to choose a model that can be updated as time progresses.
A common adhoc solution is to retrain the model after a certain amount
of time (a week, a month, etc), possibly with temporal discounting, though there is no general consensus about
how frequently the retraining should be performed.
Further there are genuinely updating models, socalled online learning models, which update model
parameters after each new match outcome is revealed.
For model evaluation, the sequential nature of the data poses a severe restriction: Any two data points were measured at certain time points, and one can not assume that they are not correlated through time information. That such correlation exists is quite plausible in the domain application, as a team would be expected to perform more similarly at close time points than at distant time points. Also, we would like to make sure that we fairly test the models for their prediction accuracy  hence the validation experiment needs to mimic the “real world” prediction process, in which the predicted outcomes will be in the temporal future of the training data. Hence the test set, in a validation experiment that should quantify goodness of such prediction, also needs to be in the temporal future of the training set.
In particular, the common independence assumption that allows application of resampling strategies such as the Kfold crossvalidation method [61]
, which guarantees the expected loss to be estimated by the empirical loss, is violated. In the presence of temporal correlation, the variance of the error metric may be underestimated, and the error metric itself will, in general, be misestimated. Moreover, the validation method will need to accommodate the fact that the model may be updated online during testing. In literature, modelindependent validation strategies for data with temporal structure is largely an unexplored (since technically difficult) area. Nevertheless, developing a reasonable validation method is crucial for scientific model assessment. A plausible validation method is introduced in section
5.2.2 in detail. It follows similar lines as the oftenseen “temporal crossvalidation” where training/test splits are always temporal, i.e., the training data points are in the temporal past of the test data points, for multiple splits. An earlier occurrence of such a validation strategy may be found in [25].This strategy comes without strong estimation guarantees and is part heuristic; the empirical loss will estimate the generalization loss as long as statistical properties do not change as time shifts forward, for example under stationarity assumptions. While this implicit assumption may be plausible for the English Premier League, this condition is routinely violated in financial time series, for example.
3 Approaches to competitive sports prediction
In this section, we give a brief overview over the major approaches to prediction in competitive sports found in literature. Briefly, these are:

The BradleyTerry models and extensions.

The Élő model and extensions.

Bayesian models, especially latent variable models and/or graphical models for the outcome and score distribution.

Supervised machine learning type models that use domain features for prediction.
(a) The BradleyTerry model is the most influential statistical approach to ranking based on competitive
observations [3].
With its original applications in psychometrics, the goal of the class of BradleyTerry models is to
estimate a hypothesized rank or skill level from observations of pairwise competition outcomes (win/loss).
Literature in this branch of research is, usually, primarily concerned not with prediction, but estimation of
a “true” rank or skill, existence of which is hypothesized, though prediction
of (binary) outcome probabilities or odds is well possible within the paradigm.
A notable exception is the work of [60] where the problem is in essence formulated
as supervised prediction, similar to our work.
Mathematically, BradleyTerry models may be seen as loglinear twofactor models that, at the stateofart are usually
estimated by (analytic or semianalytic) likelihood maximization [24].
Recent work has seen many extensions of the BradleyTerry models, most notably for modelling of ties [48],
making use of features [18] or for explicit modelling the time dependency of skill [7].
(b) The Élő system is one of the earliest attempts to model competitive sports
and, due to its mathematical simplicity, wellknown and widelyused by practitioners [15].
Historically, the Élő system is used for chess rankings, to assign a rank score to chess players.
Mathematically, the Élő system only uses information about the historical match outcomes. The Élő
system assigns to each team a parameter, the socalled Élő rating.
The rating reflects a team’s competitive skills: the team with higher
rating is stronger.
As such, the Élő system is, originally, not a predictive model or a statistical model in the usual sense.
However, the Élő system also gives a probabilistic prediction for the binary match outcome based
on the ratings of two teams.
After what appears to have been a period of parallel development that is still partly ongoing,
it has been recently noted by members of the BradleyTerry community that the Élő prediction heuristic
is mathematically equivalent to the prediction via the simple BradleyTerry
model [see 10, , section 2.1].
The Élő ratings are learnt via an update rule that is applied whenever a new outcome is observed.
This suggested update strategy is inherently algorithmic and later shown to be closely related to
online learning strategies in neural network; to our knowledge it appears first in Élő’s work
and is not found in the BradleyTerry strain.
(c) The Bayesian paradigm
offers a natural framework to model match outcomes probabilistically, and to obtain probabilistic predictions as the posterior predictive distribution. Bayesian parametric models also allow researchers to inject expert knowledge through the prior distribution. The prediction function is naturally given by the posterior distribution of the scores, which can be updated as more observations become available.
Often, such models explicitly model not only the outcome but also the score distribution,
such as Maher’s model [37] which models outcome scores
based on independent Poisson random variables with teamspecific means.
Dixon and Coles [13]
extend Maher’s model by introducing a correlation effect between
the two final scores.
More recent models also include dynamic components to model
temporal dependence [20, 50, 11].
Most models of this type only use historical match outcomes as features,
see Constantinou et al. [9] for an exception.
(d) More recently, the methodagnostic supervised machine learning paradigm has been
applied to prediction of match outcomes [36, 23, 43].
The main rationale in this branch of research is that the best model is not known, hence
a number of offshelf predictors are tried and compared in a benchmarking experiment.
Further, these models are able to make use of features other than previous outcomes easily.
However, usually, the machine learning models are trained inbatch, i.e., not following a dynamic update or online learning strategy,
and they need to be retrained periodically to incorporate new observations.
In this manuscript, we will reinterpret the Élő model and its update rule as the simplest case of a structured extension of predictive logistic (or generalized linear) regression models, and the canonical gradient ascent update of its likelihood  hence, in fact, giving it a parametric form not entirely unlike the models mentioned in (b), In the subsequent sections, this will allow us to complement it with the beneficial properties of the machine learning approach (c), most notably the addition of possibly complex features, paired with the Élő update rule which can be shown generalize to an online update strategy.
More detailed literature and technical overview is given given in the subsequent sections. The Élő model and its extensions, as well as its novel parametric interpretation, are reviewed in Section 3.1. Section 3.2 reviews other parametric models for predicting final scores. Section 3.3 reviews the use of machine learning predictors and feature engineering for sports prediction.
3.1 The BradleyTerryÉlő models
This section reviews the BradleyTerry models, the Élő system, and closely related variants.
We give the abovementioned joint formulation, following the modern rationale of considering as a “model” not only a generative specification, but also algorithms for training, predicting and updating its parameters. As the first seems to originate with the work of [3], and the second in the online update heuristic of [15], we argue that for giving proper credit, it is probably more appropriate to talk about BradleyTerryÉlő models (except in the specific hypothesis testing scenario covered in the original work of Bradley and Terry).
Later, we will attempt to understand the Élő system as an online update of a structured logistic odds model.
3.1.1 The original formulation of the Élő model
We will first introduce the original version of the Élő model, following [15]. As stated above, its original form which is still applied for determining the official chess ratings (with minor domainspecific modifications), is neither a statistical model nor a predictive model in the usual sense.
Instead, the original version is centered around the ratings for each team . These ratings are updated via the Élő model rule, which we explain (for sake of clarity) for the case of no draws: After observing a match between (home) team and (away) team , the ratings of teams and are updated as
(3)  
where , often called “the K factor”, is an arbitrarily chosen constant, that is, a model parameter usually set per hand. is if team/player has been observed to win, and otherwise.
Further, is the probability of winning against which is predicted from the ratings prior to the update by
(4) 
where is the logistic function (which has a sigmoid shape, hence is also often called “the sigmoid”). Sometimes a home team parameter is added to account for home advantage, and the predictive equation becomes
(5) 
Élő’s update rule (Equation 3) makes sense intuitively because the term can be thought of as the discrepancy between what is expected, , and what is observed, . The update will be larger if the current parameter setting produces a large discrepancy. However, a concise theoretical justification has not been articulated in literature. If fact, Élő himself commented that “the logic of the equation is evident without algebraic demonstration” [15]  which may be true in his case, but not satisfactory in an applied scientific nor a theoretical/mathematical sense.
As an initial issue, it has been noted that the whole model is invariant under joint rescaling of the , and the parameters , as well as under arbitrary choice of zero for the (i.e., adding of a fixed constant to all ). Hence, fixed domain models will usually choose zero and scale arbitrarily. In chess rankings, for example, the formula includes additional scaling constants of the form ; scale and zero are set through fixing some historical chess players’ rating, which happens to set the “interesting” range in the positive thousands^{3}^{3}3A common misunderstanding here is that no Élő ratings below zero may occur. This is, inprinciple, wrong, though it may be extremely unlikely in practice if the arbitrarily chosen zero is chosen low enough.. One can show that there are no more parameter redundancies, hence scaling/zeroing turns out not to be a problem if kept in mind.
However, three issues are left open in this formulation:

How the ratings for players/teams are determined who have never played a game before.

The choice of the constant/parameter , the “Kfactor”.

If a home parameter is present, its size.
These issues are usually addressed in everyday practice by (more or less welljustified) heuristics.
The parametric and probabilistic supervised setting in the following sections yields more principled ways to address this. step (i) will become unnecessary by pointing out a batch learning method; the constant in (ii) will turn out to be the learning rate in a gradient update, hence it can be crossvalidated or entirely replaced by a different strategy for learning the model. Parameters such as in (iii) will be interpretable as a logistic regression coefficient.
3.1.2 BradleyTerryÉlő models
As outlined in the initial discussion, the class of BradleyTerry models introduced by [3] may be interpreted as a proper statistical model formulation of the Élő prediction heuristic.
Despite their close mathematical vicinity, it should be noted that classically BradleyTerry and Élő models are usually applied and interpreted differently, and consequently fitted/learnt differently: while both models estimate a rank or score, the primary (historical) purpose of the BradleyTerry is to estimate the rank, while the Élő system is additionally intended to supply easytocompute updates as new outcomes are observed, a feature for which it has historically paid for by lack of mathematical rigour. The Élő system is often invoked to predict future outcome probabilities, while the BradleyTerry models usually do not see predictive use (despite their capability to do so, and the mathematical equivalence of both predictive rules).
However, as mentioned above and as noted for example by [10], a joint mathematical formulation can be found, and as we will show, the different methods of training the model may be interpreted as variants of likelihoodbased batch or online strategies.
The parametric formulation is quite similar to logistic regression models, or generalized linear models, in that we will use a link function and define a model for the outcome odds. Recall, the odds for a probability are
, and the logit function is
(sometimes also called the “logodds function” for obvious reasons). A straightforward calculation shows that , or equivalently, for any , i.e., the logistic function is the inverse of the logit (and vice versa by the symmetry theorem for the inverse function).Hence we can posit the following two equivalent equations in latent parameters as definition of a predictive model:
(6)  
That is, in the first equation is interpreted as a predictive probability; i.e., . The second equation interprets this prediction in terms of a generalized linear model with a response function that is linear in the . We will write
for the vector of
; hence the second equation could also be written, in vector notation, as . Hence, in particular, the matrix with entries has rank (at most) two.Fitting the above model means estimating its latent variables . This may be done by considering the likelihood of the latent parameters given the training data. For a single observed match outcome , the loglikelihood of and is
where the on the right hand side need to be interpreted as functions of (namely, as in equation 6). We call the oneoutcome loglikelihood as it is based on a single data point. Similarly, if multiple training outcomes are observed, the loglikelihood of the vector is
We will call the batch loglikelihood as the training set contains more than one data point.
The derivative of the oneoutcome loglikelihood is
hence the in the Élő update rule (see equation 3) may be updated as a gradient ascent rate or learning coefficient in an online likelihood update. We also obtain a batch gradient from the batch loglikelihood:
where, is team ’s number of wins minus number of losses observed in , and is the (multi)set of (unordered) pairings team has participated in . The batch gradient directly gives rise to a batch gradient update
Note that the above model highlights several novel, interconnected, and possibly so far unknown (or at least not jointly observed) aspects of BradleyTerry and Élő type models:

The Élő system can be seen as a learning algorithm for a logistic odds model with latent variables, the BradleyTerry model (and hence, by extension, as a full fit/predict specification of a certain onelayer neural network).

The BradleyTerry and Élő model may simultaneously be interpreted as Bernoulli observation models of a rank two matrix.

The gradient of the BradleyTerry model’s loglikelihood gives rise to a (novel) batch gradient and a singleoutcome gradient ascent update. A single iteration persample of the latter (with a fixed update constant) is Élő’s original update rule.
These observations give rise to a new family of models: the structured logodds models that will be discussed in Section 4 and 4.1, together with concomitant gradient update strategies of batch and online type. This joint view also makes extensions straightforward, for example, the “home team parameter” in the common extension of the Élő system may be interpreted as BradleyTerry model with an intercept term, with logodds , that is updated by the oneoutcome Élő update rule.
Since more generally, the structured logodds models arise by combining the parametric form of the BradleyTerry model with Élő’s update strategy, we also argue for synonymous use of the term “BradleyTerryÉlő models” whenever BradleyTerry models are updated batch, or epochwise, or whenever they are, more generally, used in a predictive, supervised, or online setting.
3.1.3 Glickman’s BradleyTerryÉlő model
For sake of completeness and comparison, we discuss the probabilistic formulation of Glickman [19]. In this fully Bayesian take on the BradleyTerryÉlő model, it is assumed that there is a latent random variable associating with team . The latent variables are statistically independent and they follow a specific generalized extreme value (GEV) distribution:
where the mean parameter varies across teams, and the other two parameters are fixed at one and zero. The density function of , is
The model further assumes that team wins over team in a match if and only if a random sample (, ) from the associated latent variables satisfies . It can be shown that the difference variables then happen to follow a logistic distribution with mean and scale parameter 1, see [49].
Hence, the (predictive) winning probability for team is eventually given by Élő’s original equation 4 which is equivalent to the BradleyTerryodds. In fact, the arguably strange parametric form for the distribution of the makes the impression of being chosen for this particular, singular reason.
We argue, that Glickman’s model makes unnecessary assumptions through the latent random variables which furthermore carry an unnatural distribution . This is certainly true in the frequentist interpretation, as the parametric model in Section 3.1.2 is not only more parsimonious as it does not assume a process that generates the , but also it avoids to assume random variables that are never directly observed (such as the ). This is also true in the Bayesian interpretation, where a prior is assumed on the which then indirectly gives rise to the outcome via the .
Hence, one may argue by Occam’s razor, that modelling the is unnecessary, and, as we believe, may put obstacles on the path to the existing and novel extensions in Section 4 that would otherwise appear natural.
3.1.4 Limitations of the BradleyTerryÉlő model and existing remedies
We point out some limitations of the original BradleyTerry and Élő models which we attempt to address in Section 4.
Modelling draws
The original BradleyTerry and Élő models do not model the possibility of a draw. This might be reasonable in official chess tournaments where players play on until draws are resolved. However, in many competitive sports a significant number of matches end up as a draw  for example, in the English Premier League about twenty percent of the matches. Modelling the possibility of draw outcome is therefore very relevant. One of the first extensions of the BradleyTerry model, the ternary outcome model by Rao and Kupper [48], was suggested to address exactly this shortcoming. The strategy for modelling draws in the joint framework, closely following this work, is outlined in Section 4.2.2.
Using final scores in the model
The BradleyTerryÉlő model only takes into account the binary outcome of the match. In sports such as football, the final scores for both teams may contain more information. Generalizations exist to tackle this problem. One approach is adopted by the official FIFA Women’s football ranking [17], where the actual outcome of the match is replaced by the "Actual Match Percentage", a quantity that depends on the final scores. FiveThirtyEight, an online media, proposed another approach [52]. It introduces the “Margin of Victory Multiplier” in the rating system to adjust the Kfactor for different final scores.
In a survey paper, Lasek et al. [35] showed empirical evidence that rating methods that take into account the final scores often outperform those that do not. However, it is worth noticing that the existing methods often rely on heuristics and their mathematical justifications are often unpublished or unknown. We describe a principled way to incorporate final scores in Section 4.2.3 into the framework, following ideas of Dixon and Coles [13].
Using additional features
The BradleyTerryÉlő model only takes into account very limited information. Apart from previous match outcomes, the only feature it uses is the identity of home and away teams. There are many other potentially useful features. For example, whether the team is recently promoted from a lowerdivision league, or whether a key player is absent from the match. These features may help make better prediction if they are properly modeled. In Section 4.2.1, we extend the BradleyTerryÉlő model to a logistic odds model that can also make use of features, along lines similar to the featuredependent models of Firth and Turner [18].
3.2 Domainspecific parametric models
We review a number of parametric and Bayesian models that have been considered in literature to model competitive sports outcomes. A predominant property of this branch of modelling is that the final scores are explicitly modelled.
3.2.1 Bivariate Poisson regression and extensions
Maher [37] proposed to model the final scores as independent Poisson random variables. If team is playing at home field against team , then the final scores and follows
where and measure the ’attack’ rates, and and measure the ’defense’ rates of the teams. The parameter
is an adjustment term for home advantage. The model further assumes that all historical match outcomes are independent. The parameters are estimated from maximizing the loglikelihood function of all historical data. Empirical evidence suggests that the Poisson distribution fits the data well. Moreover, the Poisson distribution can be derived as the expected number of events during a fixed time period at a constant risk. This interpretation fits into the framework of competitive team sports.
Dixon and Coles [13] proposed two modifications to Maher’s model. First, the final scores and are allowed to be correlated when they are both less than two. The model employs a free parameter to capture this effect. The joint probability function of is given by the bivariate Poisson distribution 7:
(7) 
where
and
The function adjusts the probability function so that drawing becomes less likely when both scores are low. The second modification is that the DixonColes model no longer assumes match outcomes are independent through time. The modified loglikelihood function of all historical data is represented as a weighted sum of loglikelihood of individual matches illustrated in equation 8, where represents the time index. The weights are heuristically chosen to decay exponentially through time in order to emphasize more recent matches.
(8) 
The parameter estimation procedure is the same as Maher’s model. Estimates are obtained from batch optimization of modified loglikelihood.
Karlis and Ntzoufras [30] explored several other possible parametrization of the bivariate Poisson distribution including those proposed by Kocherlakota and Kocherlakota [34], and Johnson et al. [28]. The authors performed a model comparison between Maher’s independent Poisson model and various bivariate Poisson models based on AIC and BIC. However, the comparison did not include the DixonColes model. Goddard [21] performed a more comprehensive model comparison based on their forecasting performance.
3.2.2 Bayesian latent variable models
Rue and Salvesen [50] proposed a Bayesian parametric model based on the bivariate Poisson model. In addition to the paradigm change, there are three major modifications on the parameterization. First of all, the distribution for scores are truncated: scores greater than four are treated as the same category. The authors argued that the truncation reduces the extreme case where one team scores many goals. Secondly, the final scores and are assumed to be drawn from a mixture model:
The component is the truncated version of the DixonColes model, and the component is a truncated bivariate Poisson distribution (7) with and equal to the average value across all teams. Thus, the mixture model encourages a reversion to the mean. Lastly, the attack parameters and defense parameters
for each team changes over time following a Brownian motion. The temporal dependence between match outcomes are reflected by the change in parameters. This model does not have an analytical posterior for parameters. The Bayesian inference procedure is carried out via Markov Chain Monte Carlo method.
Crowder et al. [11] proposed another Bayesian formulation of the bivariate Poisson model based on the DixonColes model. The parametric form remains unchanged, but the attack parameters ’s and defense parameter changes over time following an AR(1) process. Again, the model does not have an analytical posterior. The authors proposed a fast variational inference procedure to conduct the inference.
Baio and Blangiardo [1] proposed a further extension to the bivariate Poisson model proposed by Karlis and Ntzoufras [30]. The authors noted that the correlation between final scores are parametrized explicitly in previous models, which seems unnecessary in the Bayesian setting. In their proposed model, both scores are conditionally independent given an unobserved latent variable. This hierarchical structure naturally encodes the marginal dependence between the scores.
3.3 Featurebased machine learning predictors
In recent publications, researchers reported that machine learning models achieved good prediction results for the outcomes of competitive team sports. The strengths of the machine learning approach lie in the modelagnostic and datacentric modelling using available offshelf methodology, as well as the ability to incorporate features in model building.
In this branch of research, the prediction problems are usually studied as a supervised classification problem, either binary (home team win/lose or win/other), or ternary, i.e., where the outcome of a match falls into three distinct classes: home team win, draw, and home team lose.
Liu and Lai [36]
applied logistic regression, support vector machines with different kernels, and AdaBoost to predict NCAA football outcomes. For this prediction problem, the researchers hand crafted 210 features.
Hucaljuk and Rakipović [23]
explored more machine learning predictors in the context of sports prediction. The predictors include naïve Bayes classifiers, Bayes networks, LogitBoost, knearest neighbors, Random forest, and artificial neural networks. The models are trained on 20 features derived from previous match outcomes and 10 features designed subjectively by experts (such as team’s morale).
Odachowski and Grekow [43]
conducted a similar study. The predictors are commercial implementations of various Decision Tree and ensembled trees algorithms as well as a handcrafted Bayes Network. The models are trained on a subset of 320 features derived form the time series of betting odds. In fact, this is the only study so far where the predictors have no access to previous match outcomes.
Kampakis and Adamides [29] explored the possibility of predicting match outcome from Tweets. The authors applied naïve Bayes classifiers, Random forests, logistic regression, and support vector machines to a feature set composed of 12 match outcome features and a number of Tweets features. The Tweets features are derived from unigrams and bigrams of the Tweets.
3.4 Evaluation methods used in previous studies
In all studies mentioned in this section, the authors validated their new model on a real data set and showed that the new model performs better than an existing model. However, complication arises when we would like to aggregate and compare the findings made in different papers. Different studies may employ different validation settings, different evaluation metrics, and different data sets. We report on this with a focus on the following, methodologically crucial aspects:

Studies may or may not include a wellchosen benchmark for comparison. If this is not done, then it may not be concluded that the new method outperforms the stateofart, or a random guess.

Variable selection or hyperparameter tuning procedures may or may not be described explicitly. This may raise doubts about the validity of conclusions, as “handtuning” parameters is implicit overfitting, and may lead to underestimate the generalization error in validation.

Last but equally importantly, some studies do not report the error measure on evaluation metrics (standard deviation or confidence interval). In these studies, we cannot rule out the possibility that the new model is outperforming the baselines just by chance.
In table 1, we summarize the benchmark evaluation methodology used in previous studies. One may remark that the size of testing data sets vary considerably across different studies, and most studies do not provide a quantitative assessment on the evaluation metric. We also note that some studies perform the evaluation on the training data (i.e., insample). Without further argument, these evaluation results only show the goodnessoffit of the model on the training data, as they do not provide a reliable estimate of the expected predictive performance (on unseen data).
Study  Validation  Tuning  Task  Metrics  Baseline  Error  Train  Test 

Lasek et al. [35]  Online  Yes  Binary  Brier score, Binomial divergence  Yes  Yes  NA  979 
Maher [37]  Insample  No  Scores  statistic  No  No  5544  NA 
Dixon and Coles [13]  No  No  Scores  Nonstandard  No  No  NA  NA 
Karlis and Ntzoufras [30]  Insample  Bayes  Scores  AIC, BIC  No  No  615  NA 
Goddard [21]  Custom  Bayes  Scores  logloss  No  No  6930  4200 
Rue and Salvesen [50]  Custom  Bayes  Scores  logloss  Yes  No  280  280 
Crowder et al. [11]  Online  Bayes  Tenary  Accuracy  No  No  1680  1680 
Baio and Blangiardo [1]  Holdout  Bayes  Scores  Not reported  No  No  4590  306 
Liu and Lai [36]  Holdout  No  Binary  Accuracy  Yes  No  480  240 
Hucaljuk and Rakipović [23]  Custom  Yes  Binary  Accuracy, F1  Yes  No  96  96 
Odachowski and Grekow [43]  10fold CV  No  Tenary  Accuracy  Yes  No  1116  1116 
Kampakis and Adamides [29]  LOOCV  No  Binary  Accuracy, Cohen’s kappa  No  Yes  NR  NR 
4 Extending the BradleyTerryÉlő model
In this section, we propose a new family of models for the outcome of competitive team sports, the structured logodds models. We will show that both BradleyTerry and Élő models belong to this family (section 4.1), as well as logistic regression. We then propose several new models with added flexibility (section 4.2) and introduce various training algorithms (section 4.3 and 4.4).
4.1 The structured logodds model
Recall our principal observations obtained from the joint discussion of BradleyTerry and Élő models in Section 3.1.2:

The Élő system can be seen as a learning algorithm for a logistic odds model with latent variables, the BradleyTerry model (and hence, by extension, as a full fit/predict specification of a certain onelayer neural network).

The BradleyTerry and Élő model may simultaneously be interpreted as Bernoulli observation models of a rank two matrix.

The gradient of the BradleyTerry model’s loglikelihood gives rise to a (novel) batch gradient and a singleoutcome gradient ascent update. A single iteration persample of the latter (with a fixed update constant) is Élő’s original update rule.
We collate these observations in a mathematical model, and highlight relations to wellknown model classes, including the BradleyTerryÉlő model, logistic regression, and neural networks.
4.1.1 Statistical definition of structured logodds models
In the definition below, we separate added assumptions and notations for the general setup, given in the paragraph “Setup and notation”, from modelspecific assumptions, given in the paragraph “model definition”. Modelspecific assumptions, as usual, need not hold for the “true” generative process, and the mismatch of the assumed model structure to the true generative process may be (and should be) quantified in a benchmark experiment.
Setup and notation.
We keep the notation of Section 2; for the time being, we assume that there is no dependence on time, i.e., the observations follow a generative joint random variable . The variable models the outcomes of a pairing where home team plays against away team . We will further assume that the outcomes are binary home team win/lose = 1/0, i.e., . The variable models features relevant to the pairing. From it, we may single out features that pertain to a single team , as a variable . Without loss of generality (for example, through introduction of indicator variables), we will assume that takes values in , and takes values in . We will write and for the components.
The two restrictive assumptions (independence of time, binary outcome) are temporary and are made for expository reasons. We will discuss in subsequent sections how these assumptions may be removed.
We have noted that the double subindex notation easily allows to consider in matrix form. We will denote by to the (real) matrix with entry in the th row and th column. Similarly, we will denote by the matrix with entries . We do not fix a particular ordering of the entries in as the numbering of teams does not matter, however the indexing needs to be consistent across and any matrix of this format that we may define later.
A crucial observation is that the entries of the matrix can be plausibly expected to not be arbitrary. For example, if team is a strong team, we should expect to be larger for all ’s. We can make a similar argument if we know team is a weak team. This means the entries in matrix are not completely independent from each other (in an algebraic sense); in other words, the matrix can be plausibly assumed to have an inherent structure.
Hence, prediction of should be more accurate if the correct structural assumption is made on , which will be one of the cornerstones of the structured logodds models.
For mathematical convenience (and for reasons of scientific parsimony which we will discuss), we will not directly endow the matrix with structure, but the matrix where as usual and as in the following, univariate functions are applied entrywise (e.g., is also a valid statement and equivalent to the above).
Model definition.
We are now ready to introduce the structured logodds models for competitive team sports. As the name says, the main assumption of the model is that the logodds matrix is a structured matrix, alongside with the other assumptions of the BradleyTerryÉlő model in Section 3.1.2.
More explicitly, all assumptions of the structured logodds model may be written as
(9)  
satisfies certain structural assumptions 
where we have not made the structural assumptions on explicit yet. The matrix may depend on , though a sensible model may be already obtained from a constant matrix with restricted structure. We will show that the BradleyTerry and Élő models are of this subtype.
Structural assumptions for the logodds.
We list a few structural assumptions that may or may not be present in some form,
and will be key in understanding important cases of the structured logodds models.
These may be applied to as a constant matrix to obtain the simplest class of logodds models,
such as the BradleyTerryÉlő model as we will explain in the subsequent section.
Lowrankness. A common structural restriction for a matrix (and arguably the most scientifically or mathematically parsimonious one) is the assumption of low rank: namely, that the rank of the matrix of relevance is less than or equal to a specified value . Typically,
is far less than either size of the matrix, which heavily restricts the number of (model/algebraic) degrees of freedom in an
matrix from to . The lowrank assumption essentially reflects a belief that the unknown matrix is determined by only a small number of factors, corresponding to a small number of prototypical rows/columns, with the small number being equal to. By the singular value decomposition theorem, any rank
matrix may be written asfor some , pairwise orthogonal , pairwise orthogonal ;
equivalently, in matrix notation, where is diagonal, and (and where , and are the rows of ).
Antisymmetry. A further structural assumption is symmetry or antisymmetry of a matrix. Antisymmetry arises in competitive outcome prediction naturally as follows: if all matches were played on neutral fields (or if home advantage is modelled separately), one should expect that , which means the probability for team to beat team is the same regardless of where the match is played (i.e., which one is the home team). Hence,
that is, is an antisymmetric matrix, i.e., .
Antisymmetry and lowrankness. It is known that any real antisymmetric matrix always has even rank [16]. That is, if a matrix is assumed to be lowrank and antisymmetric simultaneously, it will have rank or or etc. In particular, the simplest (nontrivial) antisymmetric lowrank matrices have rank . One can also show that any real antisymmetric matrix with rank can be decomposed as
(10) 
for some , pairwise orthogonal , pairwise orthogonal ;
equivalently, in matrix notation,
where is diagonal, and (and where , and are the rows of ).
Separation. In the above, in general, the factors give rise to interaction constants (namely: ) that are specific to the pairing. To obtain interaction constants that only depend on one of the teams, one may additionally assume that one of the factors is constant, or a vector of ones (without loss of generality from the constant vector). Similarly, a matrix with constant entries corresponds to an effect independent of the pairing.
Learning/fitting of structured logodds models
will be discussed in Section 4.3. after we have established a number of important subcases and the full formulation of the model.
In a brief preview summary, it will be shown that the loglikelihood function has in essence the same form for all structured logodds models. Namely, for any parameter on which or may depend, it holds for the (oneoutcome loglikelihood) that
Similarly, for its derivative one obtains
where the partial derivatives on the right hand side will have a different form for different structural assumptions, while the general form of the formula above is the same for any such assumption.
Section 4.3 will expand on this for the full model class.
4.1.2 Important special cases
We highlight a few important special types of structured logodds models that we have already seen, or that are prototypical
for our subsequent discussion:
The BradleyTerrymodel and via identification the Élő system are obtained under the structural assumption that is antisymmetric and of rank 2 with one factor vector of ones.
Namely, recalling equation 6, we recognize that the logodds matrix in the BradleyTerry model is given by , where and are the Élő ratings. Using the rule of matrix multiplication, one can verify that this is equivalent to
where is a vector of ones and is the vector
of Élő ratings. For general , the
logodds matrix will have rank two (general = except if for all ).
By the exposition above, making the three assumptions is equivalent to positing the BradleyTerry or Élő model.
Two interesting observations may be made:
First, the onesvector being a factor entails that the winning chance
depends only on the difference between the teamspecific ratings , without any further interaction term.
Second, the entrywise exponential of is a matrix of rank (at most) one.
The popular Élő model with home advantage is obtained from the BradleyTerryÉlő model under the structural assumption that is a sum of lowrank matrix and a constant; equivalently, from an assumption of rank 3 which is further restricted by fixing some factors to each other or to vectors of ones.
More precisely, from equation 5, one can recognize that for the Élő model with home advantage, the logodds matrix decomposes as
Note that the logodds matrix is no longer antisymmetric due to the constant term
with home advantage parameter that is (algebraically) independent of the playing teams.
Also note that the antisymmetric part, i.e., ,
is equivalent to the constantfree Élő model’s logodds, while the symmetric
part, i.e., is exactly the new constant home advantage term.
More factors: full twofactor BradleyTerryÉlő models may be obtained by dropping the separation assumption from either BradleyTerryÉlő model, i.e., keeping the assumption of antisymmetric rank two, but allowing an arbitrary second factor not necessarily being the vector of ones. The team’s competitive strength is then determined by two interacting factors , , as
(11) 
Intuitively, this may cover, for example, a situation where the benefit from being much better may be
smaller (or larger) than being a little better, akin to a discounting of extremes.
If the full twofactor model predicts better than the BradleyTerryÉlő model, it may certify for
different interaction in different ranges of the Élő scores.
A home advantage factor (a constant) may or may not be added, yielding a model of total rank 3.
Raising the rank: higherrank BradleyTerryÉlő models may be obtained by model by relaxing assumption of rank 2 (or 3) to higher rank. We will consider the next more expressive model, of rank four. The rank four BradleyTerryÉlő model which we will consider will add a full antisymmetric rank two summand to the logodds matrix, which hence is assumed to have the following structure:
(12) 
The team’s competitive strength is captured by three factors , and ; note that we have kept the vector of ones as a factor. Also note that setting either of to would not result in a model extension as the resulting matrix would still have rank two. The rankfour model may intuitively make sense if there are (at least) two distinguishable qualities determining the outcome  for example physical fitness of the team and strategic competence. Whether there is evidence for the existence of more than one factor, as opposed to assuming just a single one (as a single summary quantifier for good vs bad) may be checked by comparing predictive capabilities of the respective models. Again, a home advantage factor may be added, yielding a logodds matrix of total rank 5.
We would like to note that a mathematically equivalent model, as well as models with more factors, have already been considered by Stanescu [60],
though without making explicit the connection to matrices which are of low rank, antisymmetric or structured in any other way.
Logistic regression may also be obtained as a special case of structured logodds models. In the simplest form of logistic regression, the logodds matrix is a linear functional in the features. Recall that in the case of competitive outcome prediction, we consider pairing features taking values in , and team features taking values in . We may model the logodds matrix as a linear functional in these, i.e., model that
where . If , we obtain a simple twofactor logistic regression model. In the case that there is only two teams playing only with each other, or (the mathematical correlate of) a single team playing only with itself, the standard logistic regression model is recovered.
Conversely, a way to obtain the BradleyTerry model as a special case of classical logistic regression is as follows: consider the indicator feature . With a coefficient vector , the logistic odds will be . In this case, the coefficient vector corresponds to a vector of Élő ratings.
Note that in the above formulation, the coefficient vectors are explicitly allowed to depend on the teams. If we further allow to depend on both teams, the model includes the BradleyTerryÉlő models above as well; we could also make the depend on both teams. However, allowing the coefficients to vary in full generality is not very sensible, and as for the constant term which may yield the Élő model under specific structural assumptions, we need to endow all model parameters with structural assumptions to prevent combinatorial explosion of parameters and overfitting.
These subtleties in incorporating features, and more generally how to combine features with hidden factors will be discussed in the separate, subsequent Section 4.2.1.
4.1.3 Connection to existing model classes
Close connections to three important classes of models become apparent through the discussion in the previous sections:
Generalized Linear Models generalize both linear and loglinear models (such as the BradleyTerry model) through socalled link functions, or more generally (and less classically) link distributions, combined with flexible structural assumptions on the target variable. The generalization aims at extending prediction with linear functionals through the choice of link which is most suitable for the target [for an overview, see 39].
Particularly relevant for us are generalized linear models for ordinal outcomes which includes the
ternary (win/draw/lose) case, as well as link distributions for scores. Some existing extensions of this type,
such as the ternay outcome model of [48] and the score model of [37],
may be interpreted as specific choices of suitable linking distributions.
How these ideas may be used as a component of structured logodds models will be discussed in Section 4.2.
Neural Networks
(vulgo “deep learning”) may be seen as a generalization of logistic regression which is mathematically equivalent to a singlelayer network with softmax activation function. The generalization is achieved through functional nesting which allows for nonlinear prediction functionals, and greatly expands the capability of regression models to handle nonlinear featurestargetrelations
[for an overview, see 51].A family of ideas which immediately transfers to our setting are strategies for training and model fitting. In particular, online update strategies as well as training in batches and epochs yields a natural and principled way to learn BradleyTerryÉlő and logodds models in an online setting or to potentially improve its predictive power in a static supervised learning setting. A selection of such training strategies for structured logodds models will be explored in Section 4.3
. This will not include variants of stochastic gradient descent which we leave to future investigations.
It is also beyond the scope of this manuscript to explore the implications of using multiple layers
in a competitive outcome setting, though it seems to be a natural idea given the closeness of the model classes
which certainly might be worth exploring in further research.
Lowrank Matrix Completion is the supervised task of filling in some missing entries of a lowrank matrix, given others and the information that the rank is small. Many machine learning applications can be viewed as estimation or completion of a lowrank matrix, and different solution strategies exist [4, 6, 42, 32, 40, 55, 63, 33].
The featurefree variant of structured logodds models (see Section 4.1.1) may be regarded as a lowrank matrix completion problem: from observations of for where the set of observed pairings may be considered as the set of observed positions, estimate the underlying lowrank matrix , or predict for some which is possibly not contained in .
One popular lowrank matrix completion strategy in estimating model parameters or completing missing entries uses the idea of replacing the discrete rank constraint by a continuous spectral surrogate constraint, penalizing not rank but the nuclear norm ( = trace norm = 1Schattennorm) of the matrix modelled to have low rank [an early occurrence of this idea may be found in 57]. The advantage of this strategy is that no particular rank needs to be apriori assumed, instead the objective implicitly selects a low rank through a tradeoff with model fit. This strategy will be explored in Section 4.4 for the structured logodds models.
Further, identifiability of the structured logodds models is closely linked to the question whether a given entry of a lowrank matrix may be reconstructed from those which have been observed. Somewhat straightforwardly, one may see that reconstructability in the algebraic sense, see [33], is a necessary condition for identifiability under respective structure assumptions. However, even though many results of [33] directly generalize, completability of antisymmetric lowrank matrices with or without vectors of ones being factors has not been studied explicitly in literature to our knowledge, hence we only point this out as an interesting avenue for future research.
4.2 Predicting nonbinary labels with structured logodds models
In Section 4.1, we have not introduced all aspects of structured logodds models in favour of a clearer exposition. In this section, we discuss these aspects that are useful for the domain application more precisely, namely:

How to use features in the prediction.

How to model ternary match outcomes (win/draw/lose) or score outcomes.

How to train the model in an online setting with a batch/epoch strategy.
For point (i) “using features”, we will draw from the structured logodds models’ closeness to logistic regression; the approach to (ii) “general outcomes” may be treated by choosing an appropriate link function as with generalized linear models; for (iii), parallels may be drawn to training strategies for neural networks.
4.2.1 The structured logodds model with features
As highlighted in Section 4.1.2, pairing features taking values in , and team features taking values in may be incorporated by modelling the logodds matrix as
(13) 
where . Note that differently from the simpler exposition in Section 4.1.2, we allow all coefficients, including , to vary with and .
Though, allowing and to vary completely freely may lead to overparameterisation or overfitting, similarly to an unrestricted (full rank) logodds matrix of in the lowrank Élő model, especially if the number of distinct observed pairings is of similar magnitude as the number of total observed outcomes.
Hence, structural restriction of the degrees of freedom may be as important for the feature coefficients as for the constant term. The simplest such assumption is that all are equal, all are equal, and all are equal, i.e., assuming that
for some and where may follow the assumptions of the featurefree logodds models. This will be the main variant which will refer to as the structured logodds model with features.
However, the assumption that constants are independent of the pairing may be too restrictive, as it may be plausible that, for example, teams of different strength profit differently from or are impaired differently by the same circumstance, e.g., injury of a key player.
To address such a situation, it is helpful to rewrite Equation 13 in matrix form:
where is the matrix whose rows are the , where and are matrices whose rows are the , and where is the matrix with entries . The symbols and
denote tensors of degree 3 (= 3Darrays) whose
th elements are and . The symbol stands for the indexwise product of degree3tensors which eliminates the third index and yields a matrix, i.e.,A natural parsimony assumption for , and is, again, that of lowrank. For the matrices, , one can explore the same structural assumptions as in Section 4.1.1: lowrankness and factors of one are reasonable to assume for all three, while antisymmetry seems natural for but not for .
A low tensor rank (Tucker or Waring) appears to be a reasonable assumption for . As an adhoc definition of tensor (decomposition) rank of , one may take the minimal such that there is a decomposition into real vectors such that
Further reasonable assumptions are antisymmetry in the first two indices, i.e., , as well as some factors being vectors of ones.
Exploring these possible structural assumptions on the coefficients of features in experiments is possibly interesting both from a theoretical and practical perspective, but beyond the scope of this manuscript. Instead, we will restrict ourselves to the case of , of and having the same entry each, and following one of the lowrank assumptions in structural assumptions as in Section 4.1.1 as in the featurefree model.
We would like to note that variants of the BradleyTerry model with features have already been proposed and implemented in the BradleyTerry2 package for R [18], though isolated from other aspects of the BradleyTerryÉlő model class such as modelling draws, or structural restrictions on hidden variables or the coefficient matrices and tensors, and the Élő online update.
4.2.2 Predicting ternary outcomes
This section addresses the issue of modeling draws raised in 3.1.4. When it is necessary to model draws, we assume that the outcome of a match is an ordinal random variable of three socalled levels: win draw lose. The draw is treated as a middle outcome. The extension of structured logodds model is inspired by an extension of logistic regression: the Proportional Odds model.
The Proportional Odds model is a wellknown family of models for ordinal random variables [38]. It extends the logistic regression to model ordinary target variables. The model parameterizes the logit transformation of the cumulative probability as a linear function of features. The coefficients associated with feature variables are shared across all levels, but there is an intercept term which is specific to a certain level. For a generic featurelabel distribution , where takes values in and takes values in a discrete set of ordered levels, the proportional odds model may be written as
where , and . The model is called Proportional Odds model because the odds for any two different levels , , given an observed feature set, are proportional with a constant that does not depend on features; mathematically,
Using a similar formulation in which we closely follow Rao and Kupper [48], the structured logodds model can be extended to model draws, namely by setting
where is the entry in structured logodds matrix and is a free parameter that affects the estimated probability of a draw. Under this formulation, the probabilities for different outcomes are given by
Note that this may be seen as a choice of ordinal link distribution in a “generalized” structured odds model, and may be readily combined with feature terms as in Section 4.2.1.
4.2.3 Predicting score outcomes
Several models have been considered in Section 3.1.4 that use score differences to update the Élő ratings. In this section, we derive a principled way to predict scores, score differences and/or learn from scores or score differences.
Following the analogy to generalized linear models, we will be able to tackle this by using a suitable linking distribution, the model can utilize additional information in final scores.
The simplest natural assumption one may make on scores is obtained from assuming a dependent scoring process, i.e., both home and away team’s scores are Poissondistributed with a teamdependent parameter and possible correlation. This assumption is frequently made in literature [37, 13, 11] and eventually leads to a (double) Poisson regression when combined with structured logodds models.
The natural linking distributions for differences of scores are Skellam distributions which are obtained as difference distributions of two (possibly correlated) Poisson distributions [54], as it has been suggested by Karlis and Ntzoufras [31].
In the following, we discuss only the case of score differences in detail, predicting both team’s score distributions can be obtained similarly as predicting the correlated Poisson variables with the respective parameters instead of the Skellam difference distribution.
We first introduce some notation. As a difference of Poisson distributions whose support is , the support of a Skellam distribution is the set of integers . The probability mass function of Skellam distributions takes two positive parameters and , and is given by
where is the modified Bessel function of first kind with parameter , given by
If random variables and follow Poisson distributions with mean parameters and respectively, and their correlation is , then their difference follows a Skellam distribution with mean parameters and .
Now we are ready to extend the structured logodds model to incorporate historical final scores. We will use a Skellam distribution as the linking distribution: we assume that the score difference of a match between team and team , that is, (taking values in ), follows a Skellam distribution with (unknown) parameter and .
Note that hence there are now two structured , each of which may be subject to constraints such as in Section 4.1.1, or constraints connecting them to each other, and each of which may depend on features as outlined in Section 4.2.1.
A simple (and arguably the simplest sensible) structural assumption is that , is rank two, with factors of ones, as follows:
equivalently, that has rank one and only nonnegative entries.
As mentioned above, features such as home advantage may be added to the structured parameter matrix or using the way introduced in Section 4.2.1.
Also note that the above yields a strategy to make ternary predictions while training on the scores. Namely, a prediction for ternary match outcomes may simply be derived from predicted score differences , through defining
Comments
There are no comments yet.