Modeling Intransitivity in Pairwise Comparisons with Application to Baseball Data

Abstract

The seminal Bradley-Terry model exhibits transitivity, that is, the property that the probabilities of player A beating B and B beating C give the probability of A beating C, with these probabilities determined by a skill parameter for each player. Such transitive models do not account for different strategies of play between each pair of players, which gives rise to intransitivity. Various intransitive parametric models have been proposed but they lack the flexibility to cover the different strategies across n players, with the $O(n^2)$ values of intransitivity modeled using $O(n)$ parameters, while they are not parsimonious when the intransitivity is simple. We overcome their lack of adaptability by allocating each pair of players to one of a random number of K intransitivity levels, each level representing a different strategy. Our novel approach for the skill parameters involves having the n players allocated to a random number of $A<n$ distinct skill levels, to improve efficiency and avoid false rankings. Although we may have to estimate up to $O(n^2)$ unknown parameters for $(A,K)$ we anticipate that in many practical contexts $A+K<n$. Our semiparametric model, which gives the Bradley-Terry model when $(A=n-1,K=0)$, is shown to have an improved fit relative to the Bradley-Terry, and the existing intransitivity models, in out-of-sample testing when applied to simulated and American League baseball data. Supplementary materials for the article areavailable online.

Keywords:

• Baseball
• Bayesian hierarchical modeling
• Bradley-Terry
• Clustering
• Intransitivity
• Pairwise Comparisons
• Ranking
• Reversible jump Markov chain Monte Carlo
• Tournament structure

1 Introduction

The seminal Bradley-Terry model (Bradley and Terry 1952) is commonly used to rank objects from paired comparison data. Given a set $\mathcal{I}$ of n objects with each object $i\in \mathcal{I}$ having skill $r_i\in \mathbb{R}$, then the Bradley-Terry model gives, for $i\ne j\in \mathcal{I}$,(1) $$p_{ij}^{(\text{BT})}=\mathrm{Pr}\{i\succ j\}:={\{1+\hspace{0.17em}\exp [-(r_i-r_j)\left]\right\}}^{-1},$$(1) where $a\succ b$ denotes preference for object a over b, and $r_1=0$ to avoid identifiability issues. A ranking of the objects is given by sorting estimates of $r:=\{r_i\in \mathbb{R}:i\in \mathcal{I}\}$. This model is transitive, that is, $p_{jk}^{(\text{BT})}$ is given by $p_{ij}^{(\text{BT})}$ and $p_{ik}^{(\text{BT})}$, for all $i\ne j\ne k\in \mathcal{I}$, see Section 3.

Now consider the game of Rock-Paper-Scissors, a zero-sum game in which Rock beats Scissors, Scissors beats Paper, and Paper beats Rock, and specifically consider the deterministic scenario where players (r,p,s) always pick (Rock, Paper, Scissors), respectively. In this scenario, all win probabilities in a game are either 0 or 1 depending on the opponent, and each player wins their next game with probability $1/2$ if their next opponent is to be selected at random. Whatever way the skill of a player is defined, the symmetry of this game setup unquestionably leads to the conclusion that the three players have equalskill levels.

Conclusions drawn from a Bradley-Terry model fitted to data from this simple game are surprisingly poor. Given a round-robin tournament, where each player plays all other players an equal number of times, the model will correctly estimate that all players are equally ranked in terms of skills; however, it would also estimate all pairwise win probabilities to be 1/2, which couldn’t be more wrong. Even worse, is that any illusory ranking can result when the tournament is not round-robin, for example, if the most common pairing of players is (r,s) and the other two pairings occur equally often then the Bradley-Terry model will rank player r as top. The key reason for the failure of the Bradley-Terry model is its transitive nature, a trait shared by almost all commonly used ranking systems.

Here we develop a novel pairwise comparison model, and an associated ranking system, which accounts for intransitivity. Thus, it describes how specific pairwise probabilities differ from probabilities given by overall skill levels alone, that is, how probabilities differ from those given by the Bradley-Terry model. The Rock-Paper-Scissors game also illustrates that ranking can involve ties, where subsets of players can have equal skill levels, and that tournament structure can effect the subsequent inference. We also address some aspects associated withthese issues.

The concept and associated modeling of intransitivity is not new. Makowski and Piotrowski (2006) present many examples of competitions exhibiting intransitivity and argue that it can occur whenever the best strategy in a given comparison depends on the strategy of the opponent, and Smead (2019) provides a philosophical argument as to why intransitivity is particularly likely to occur in sports. Given this, it is not surprising to find cases of intransitivity in e-sports (Makhijani and Ugander 2019; Chen and Joachims 2016; Duan et al. 2017). Other applications include social choice, real sensory analysis, and electiondatasets.

With n competitors there are $n(n-1)/2$ interactions, or intransitivities, so even in round-robin competitions, with m rounds, there are too many terms to estimate efficiently using empirical methods, unless $m/n$ is large. Causeur and Husson (2005) proposed an $O(n^2)$ parameter extension of the Bradley-Terry model to address intransitivity. Subsequently $O(nd)$ parametric models have been studied for some fixed $d\in \mathbb{N}$ $(d\ll n)$, see all the models in Section 2, but they lack the flexibility to cover the potentially $O(n^2)$ different intransitivities across n players, leading to bias; whilst they are not parsimonious when the intransitivity is simple, leading to inefficiency.

Although the concept of intransitivity is quite clear, there is no established measure of the amount of intransitivity in a dataset. In this work, we propose a definition of intransitivity through a distance metric between the assumed probability of paired comparisons under a Bradley-Terry model, and the empirical or model-based probability estimate, such that for any given dataset the magnitude of the intransitivity present is unambiguous. A flexible model then, is one which is capable of exploring the space of all possible combinations of intransitivity, as defined by this measure. Any parametric model is restricted to a subset of this space by definition, with this restriction being most obviously revealed when assessing predictive performance.

We then develop a novel semiparametric extension of the Bradley-Terry model, allocating the $n(n-1)/2$ pairs of objects to a random number K, with $0\le K\le n(n-1)/2$, of distinct intransitivity levels, each level representing a different strategy. We term this model the Intransitive Clustered Bradley-Terry (ICBT) model. Relative to the aforementioned parametric models, this ICBT model provides greater flexibility to enable the incorporation of varying structures, and degrees of, intransitivity. As many of these strategies will have similar effects, we anticipate that K should be small, yet the random property of K provides the potential for it to be large when required. This flexibility ensures that our model is parsimonious, whatever the complexity of the data. For our Rock-Paper-Scissors illustration $K=1$.

Moreover, our novel approach for the objects’ skills is to allocate the n objects into a random number of $A+1\le n$ distinct skill levels, to improve efficiency and avoid false rankings. This constraint recognizes that from paired comparison data there will be objects that are indistinguishable as having statistically significantly different skill levels, for example, for our Rock-Paper-Scissors illustration $A=0$. So clustering skills avoids overinterpretation of misinformed rankings, a feature Masarotto and Varin (2012) address by clustering skills via a lasso procedure.

The basis of our model is the belief that in practice there are likely to small subsets of skill and intransitivity levels, namely $A\le n-1$ and $K\ll n(n-1)/2$, respectively. As we have little prior knowledge about the skills of the objects or the intransitivities of the pairs of objects, we allow the clustering of objects into different skill levels, and of the pairs of objects into separate intransitivity levels, to be determined entirely through a Bayesian hierarchical model. We take each of $(A,K)$, the allocations of objects to skill levels, and the allocations of the pairs of objects to intransitivity levels as unknown, with inference being conducted via a reversible jump Markov chain Monte Carlo (RJMCMC) algorithm. This formulation does offer computational challenges; however, we anticipate that typically the posterior will give a high probability that $A+K<n$ and that many of the cluster allocations also will be strongly identified. Our inference framework offers the opportunity to select a highly simplified model, with the values of A, K and allocations fixed at values given by posterior means/modes if these are found to align with known structure about the paired comparison. In the absence of such knowledge our results allow for the full uncertainty of these features to be accounted for.

In certain circumstances our model has the potential to identify and correct for imbalanced tournament structure on overall rankings since teams are not penalized if they (unfairly) compete most frequently against those whom they perform systematically worse to relative to what is expected based on respective skills alone.

We use American League Baseball data to illustrate the performance of our methods in comparison to existing models for a range of reasons. First, each game results in a win or a loss for a team. Second, it is known to be a highly strategic game, see Section 5, so we anticipate that the level of intransitivity will be high. Finally, although the tournament structure is not round robin, each team plays each other team often, and so the existence of intransitivity should become apparent in inference. Indeed this is found in Section 5, where our model is shown to have an improved fit over the Bradley-Terry model and existing parametric intransitivity models in out of sample testing for each of the nine seasons we study.

The layout is as follows. Section 2 introduces other approaches to modeling intransitivity. Section 3 then introduces our novel measure of intransitivity, the ICBT model, and the ranking formulation. Section 4 contains details of the inference, including prior specification, our full Bayesian hierarchical modeling strategy, an overview of the RJMCMC algorithm and its novel features, and an overview of a simulation study. Section 5 compares this model with the Bradley-Terry model and other competitor models, using American League baseball data. Section 6 is a discussion. Full details of the RJMCMC algorithm, simulation study, and extended analysis of the baseball application are in the supplementary material.

2 Literature on Intransitive Modeling

The blade-chest model of Chen and Joachims (2016), extends the Bradley-Terry model into d-dimensions by incorporating so-called blade and chest vectors $b_i,c_i\in {\mathbb{R}}^d$ for each object $i\in \mathcal{I}$. There are two versions: the -dist and -inner variants, given, respectively, by$$\begin{array}{cc}\text{logit}\left(p_{ij}^{(\text{BCD})}\right)\hfill & :=\left|\right|b_j-c_i|{|}_2^2-||b_i-c_j|{|}_2^2+r_i-r_j,\hfill \\ \mathit{and}\hspace{1em}\text{logit}\left(p_{ij}^{(\text{BCI})}\right):=b_i^T\cdot c_j-b_j^T\cdot c_i+r_i-r_j.\hfill & \hfill \end{array}$$

The blade and chest parameters of all the objects can be viewed as features on a d-dimensional latent space. Then, if an object i’s blade is close to object j’s chest, and simultaneously object i’s chest is far from object j’s blade, then object i has an additional advantage over object j. If $d=2$ this model can represent a deterministic Rock-Paper-Scissors game, by placing the blade of Rock at the chest of Scissors, the blade of Scissors at the chest of Paper, and the blade of Paper at the chest of Rock. By increasing d, ever more complex relationships can be captured between the pairs of objects. Given n objects and $b_i,c_i\in {\mathbb{R}}^d$ for each object i, the model contains $2d(n-1)$ identifiable parameters. The r parameters can be absorbed into the blade and chest parameters; however, the above parametrisation makes it clear that the Bradley-Terry model is a special case of the blade-chest model, when $b_i=b_j=c_i=c_j,\forall i,j\in \mathcal{I}$.

Duan et al. (2017) introduce a generalized model for intransitivity, with$$\text{logit}\left(p_{ij}^{(G)}\right)={\mu }_i^T\Sigma {\mu }_j+{\mu }_i^T\Gamma {\mu }_i-{\mu }_j^T\Gamma {\mu }_j,$$with ${\mu }_i\in {\mathbb{R}}^d$, where d is even, being a d-dimensional strength vector, for an object $i\in \mathcal{I}$, and $\Sigma ,\Gamma \in {\mathbb{R}}^{d\times d}$ are so-called transitive matrices. The first matrix $\Sigma$ represents the interactions between objects, and $\Gamma$ controls how an individual object’s strength components form the object’s overall strength. The number of identifiable parameters is $d(3d/2+n-1)$, since there are two $d\times d$ matrices $(\Sigma ,\Gamma )$, of which $\Sigma$ is skew-symmetric, and n d-dimensional vectors. They show that their model is a generalization of the blade-chest and Bradley-Terry models. Specifically, the blade-chest-inner model arises when ${\mu }_i=(b_i,c_i),$ ${\mu }_j=(b_j,c_j)$, and $\left|\right|{\mu }_i|{|}_2^2=|\left|{\mu }_j\right|{|}_2^2$, that is, the objects all have equal skill, and $\Sigma$ is a block diagonal matrix with two $(d/2)\times (d/2)$ matrices of zeros on the diagonal and matrices $I_{d/2}$ and $-I_{d/2}$ on the off-diagonals where $I_m$ is the $m\times m$ identity matrix, then ${\mu }_i^T\Sigma {\mu }_j=b_i^T\cdot c_j-b_j^T\cdot c_i.$ The degrees of freedom are restricted by regularization, using an $L_2$ norm for the object strength vectors and Frobenius norm for both transitivity matrices. The tuning parameters are selected via cross-validation.

Makhijani and Ugander (2019) introduced a majority vote model with object i having a vector of d skill attributes, $({\mu }_{i,1},\dots ,{\mu }_{i,d})$, where d is odd. Then, given a suitable choice of mapping function f, for example, logistic or Gaussian, define $q_{ij}^l=f({\mu }_{i,l}-{\mu }_{j,l}),$ $\forall l\in \{1,\dots ,d\}$ to be the probability of i beating j based only on their lth attribute. Then, majority vote model says that the probability of i being preferred to j overall, is the probability that it wins across the majority of attributes. For $d=1$ the model is linearly transitive, but not when $d=3$, as$$\mathrm{Pr}\{i\succ j\}=q_{ij}^1{q}_{ij}^2{q}_{ij}^3+(1-q_{ij}^1)q_{ij}^2{q}_{ij}^3+q_{ij}^1(1-q_{ij}^2)q_{ij}^3+q_{ij}^1{q}_{ij}^2(1-q_{ij}^3).$$

3 Modeling

3.1 Measure of Intransitivity

From the model definition (1), the Bradley-Terry model assumes linear transitivity. This assumption constrains the pairwise probabilities of the model such that, given $p_{ij}^{(\text{BT})}$ and $p_{jk}^{(\text{BT})}$ from (1) for any i ≠ j ≠ k ∈ ℐ, the probability $p_{ik}^{(\text{BT})}$ is completely determined. It is straightforward to show that the form of $p_{ik}^{(\text{BT})}$ is given as$$p_{ik}^{(\text{BT})}=\frac{p_{ij}^{(\text{BT})}{p}_{jk}^{(\text{BT})}}{1+2p_{ij}^{(\text{BT})}{p}_{jk}^{(\text{BT})}-\left(p_{ij}^{(\text{BT})}+p_{jk}^{(\text{BT})}\right)},\forall j\ne i,k,$$noting it is independent of the choice of bridge object j. Therefore, there can be no interaction that is specific to the pair $\{i,k\}$, that is not already captured between all other pairs.

Including intransitivity; however, allows for some pairwise probabilities to depart from those assumed by the Bradley-Terry model. This can be formalized by supposing that for all $i\ne k\in \mathcal{I}$ the true probability of preference $i\succ k$ is given as some function $f:\{[0,1],\mathbb{R}\}\to [0,1]$ of the Bradley-Terry probability and the intransitivity, ${\theta }_{ik}$, of the pair $\{i,k\}$, then we can write(2) $$p_{ik}:=f\left(p_{ik}^{(\text{BT})},{\theta }_{ik}\right),\forall i\ne k\in \mathcal{I},$$(2) where we identify the form of f in Section 3.2. We define the intransitivity to be the displacement of the true probabilities from the Bradley-Terry probabilities on the log-odds scale, so that(3) $${\theta }_{ik}:=\hspace{0.17em}\log \hspace{0.17em}\left(\frac{p_{ik}/(1-p_{ik})}{p_{ik}^{(\text{BT})}/(1-p_{ik}^{(\text{BT})})}\right),\forall i\ne k\in \mathcal{I},$$(3) is the amount of intransitivity between the pair of objects $\{i,k\}$. A value of ${\theta }_{ik}=0$ indicates that the comparison is transitive, that is, the pairwise probabilities could be modeled by the Bradley-Terry model. As a consequence we require $f(x,0)=x,x\in [0,1]$. The choice of log-odds ratio in Equation (3)(3) $${\theta }_{ik}:=\hspace{0.17em}\log \hspace{0.17em}\left(\frac{p_{ik}/(1-p_{ik})}{p_{ik}^{(\text{BT})}/(1-p_{ik}^{(\text{BT})})}\right),\forall i\ne k\in \mathcal{I},$$(3) reflects the nonlinearity of probabilities. For example, if $ϵ=0.099$, then a small linear shift in probability from 0.5 to $0.5+ϵ$ has little impact on the odds, which remain at approximately 1:2. However, a linear shift in probability from 0.9 to $0.9+ϵ$ has a huge impact on the odds, which move from 1:10 to 1:1000. Moreover, our definition (3) for the intransitivity ${\theta }_{ik}$ also imposes rotational symmetry for pairs of objects, that is(4) $$\begin{array}{cc}{\theta }_{ki}\hfill & =\hspace{0.17em}\log \hspace{0.17em}\left(\frac{p_{ki}/(1-p_{ki})}{p_{ki}^{(\text{BT})}/(1-p_{ki}^{(\text{BT})})}\right)\hfill \end{array}$$(4) $$=\hspace{0.17em}\log \hspace{0.17em}\left(\frac{(1-p_{ik})/p_{ik}}{(1-p_{ik}^{(\text{BT})})/p_{ik}^{(\text{BT})}}\right)=-{\theta }_{ik},\forall i\ne k\in \mathcal{I},$$ so we need to find $\{{\theta }_{ik},\forall i>k\in \mathcal{I}\}$ only, in order to completely define $\{{\theta }_{ik},\forall i\ne k\in \mathcal{I}\}$.

3.2 Model Formulation

To find the function f in Equation (2)(2) $$p_{ik}:=f\left(p_{ik}^{(\text{BT})},{\theta }_{ik}\right),\forall i\ne k\in \mathcal{I},$$(2) , Equation (3)(3) $${\theta }_{ik}:=\hspace{0.17em}\log \hspace{0.17em}\left(\frac{p_{ik}/(1-p_{ik})}{p_{ik}^{(\text{BT})}/(1-p_{ik}^{(\text{BT})})}\right),\forall i\ne k\in \mathcal{I},$$(3) can be simply rearranged which gives(5) $$p_{ik}=\frac{p_{ik}^{(\text{BT})}\hspace{0.17em}\exp \hspace{0.17em}({\theta }_{ik})}{p_{ik}^{(\text{BT})}\hspace{0.17em}\exp \hspace{0.17em}({\theta }_{ik})+1-p_{ik}^{(\text{BT})}},\forall i\ne k\in \mathcal{I},$$(5) and so for any pair $\{i,k\}$, Equation (5)(5) $$p_{ik}=\frac{p_{ik}^{(\text{BT})}\hspace{0.17em}\exp \hspace{0.17em}({\theta }_{ik})}{p_{ik}^{(\text{BT})}\hspace{0.17em}\exp \hspace{0.17em}({\theta }_{ik})+1-p_{ik}^{(\text{BT})}},\forall i\ne k\in \mathcal{I},$$(5) can be rewritten as(6) $$p_{ik}=\frac{1}{1+\hspace{0.17em}\exp [-\left({\theta }_{ik}+r_i-r_k\right)]},\forall i\ne k\in \mathcal{I}.$$(6)

Here the effect of ${\theta }_{ik}$ is clear, positive (negative) ${\theta }_{ik}$, increases (decreases) the probability of team i beating team k relative to their skills alone, that is, relative to the Bradley-Terry model.

Thus far, the model contains the flexibility to describe $P:=\{p_{ik}\in [0,1],\forall i\ne k\in \mathcal{I}\}$ completely. Here P contains $n(n-1)/2$ degrees of freedom, because $p_{ki}=1-p_{ik}$; however, the model contains $n(n-1)/2+n$ parameters: $n(n-1)/2$ values of intransitivity between pairs, and n skill parameters from r, and thus, the model parameters are not identifiable. One way of ensuring identifiability in the standard Bradley-Terry model is to fix one object’s skill level, and here it is chosen that $r_1=0$. As well as this constraint on the objects’ skill parameters, the intransitivity parameters require constraints for parameter identifiability. The minimal set of required constraints is identified in Proposition 1, see the Appendix for the proof.

Proposition 1.

Consider a round-robin tournament with pairs of objects $(i,j)$ being compared, with $i,j\in \mathcal{I}$, with $i\ne j$ where $\left|\mathcal{I}\right|=n$. Suppose that the probabilities of i beating j are $p_{ij}$ where these probabilities are given by expression (6), with $r_1=0$ and intransitivity values ${\theta }_{ij}$. If a set of $n-1$ pairs of objects, indexed by ${\mathcal{J}}_{n-1}$, have their intransitivity values set to arbitrary specified values, then all the rest of the $\left\{r_i\right\}$ and $\left\{{\theta }_{ij}\right\}$ parameters in expression (6) are identifiable if ${\mathcal{J}}_{n-1}$ forms a connected graph over $\mathcal{I}$. Furthermore, if less than $n-1$ pairs’ intransitivity values are specified or if ${\mathcal{J}}_{n-1}$ is not a connected graph over $\mathcal{I}$ then identifiability is not achievable.

We choose the $n-1$ constraints to be ${\theta }_{ij}=0,\forall (i,j):i=1,j\in \mathcal{I}\setminus \left\{1\right\}$, that is, all pairs involving object 1 have intransitivity set to 0. Proposition 1 gives that if any further constraints are imposed on the intransitivity values the flexibility of model (6) will be compromised.

With the above constraints, the minimal conditions for parameter identifiability are satisfied, but the model is still likely to overfit with so many parameters. To rectify this we restrict the total number of degrees of freedom, by restricting both the number of intransitivity values to only $K\le n(n-1)/2$ unique values and restricting the number of unknown skill values to be $A<n$, where $A+K\le n(n-1)/2$ and ideally $A+K\ll n(n-1)/2.$ In this fashion our ICBT model embraces intransitivity in a parsimonious way.

First, consider the $A+1$ unique skill values, which ensures parsimony in the model by clustering the objects’ skills r into distinct values which are sufficiently statistically significantly different. Since $r_1=0$ is fixed, there are only A unknown skill levels, $\varphi \in {\mathbb{R}}^A$. By defining the labels of the set of skill levels to be $\mathcal{A}:=\{-A_{-},\dots ,0,\dots ,A_{+}\}$ with $A_{+}$ being the number of skill levels which are greater than 0 and $A_{-}$ the number of skill levels less than 0 such that $A_{-}+A_{+}+1=A+1=\left|\mathcal{A}\right|\le n$, we impose the equivalent condition in our model by fixing the skill level with label $\left\{0\right\}$ to be ${\varphi }_0=0$, and fixing object 1 to always be allocated to this cluster. The possible skill values an object can take are, therefore, defined as$$\begin{array}{c}\{{\varphi }_0=0,\varphi :=\{{\varphi }_a\in \mathbb{R},\forall a\in \mathcal{A}\setminus \left\{0\right\}\}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}:{\varphi }_{-A_{-}}<\dots <{\varphi }_0<\dots <{\varphi }_{A_{+}}\},\end{array}$$where $\varphi$ are the unknown skill levels, and the ordering helps with label switching problems in the inference. The skill cluster allocation of object i, denoted $y_{\left\{i\right\}}\in {\{0,1\}}^{A+1}$, is a binary vector which takes the value 1 at position $s\in \mathcal{A}$ and 0 everywhere else, if object $i\in \mathcal{I}$ belongs to cluster s. The set $Y:=\{y_{\left\{i\right\}}:i\in \mathcal{I}\setminus \{1\left\}\right\}$ then contains all the objects’ skill cluster allocations except object $\left\{1\right\}$ which has fixed cluster allocation. Therefore, by defining ${\mathcal{S}}_{\left\{i\right\}}(Y):={\text{argmax}}_s{y}_{\left\{i\right\},s},\forall i\in \mathcal{I}\setminus \left\{1\right\}$, then the objects’ skills can be written as(7) $$r_i=\{\begin{array}{cc}{\varphi }_{{\mathcal{S}}_{\left\{i\right\}}\left(Y\right)},& i\in \mathcal{I}\setminus \left\{1\right\}\\ 0,& i=1\end{array}\hspace{1em}:=f_r\left(\varphi ,Y,i\right),\forall i\in \mathcal{I}.$$(7)

Now consider the K unique values of intransitivity to describe the different inter-object strategies. Of the $n(n-1)/2$ pairs of objects, many will adopt similar strategies depending on their opponents. These similar strategies are translated by the model as having similar departures from transitivity, and are, thus, clustered together. For example, suppose some group of objects $\mathcal{V}:j\notin \mathcal{V}$ competed against object j in the same way. Then it would be reasonable to assume that ${\theta }_{ij}$ is the same for all $i\in \mathcal{V}$. This creates clusters of pairs of objects, such that the pairs are clustered according to them having identical intransitivity.

In order to measure the departure from a Bradley-Terry model, a linearly transitive cluster is imposed, which contains the set of pairs ${\mathcal{J}}_T\subseteq \left\{\right\{i,k\}:i\ne k\in \mathcal{I}\},$ which have an intransitivity level ${\theta }_0=0$. Thus, the Bradley-Terry modeling assumption (1) holds for these pairs, such that $p_{ik}=p_{ik}^{(\text{BT})}$, for all $\{i,k\}\in {\mathcal{J}}_T$. Given the existence of this cluster, there must be strong evidence from the data to produce an additional cluster with an intransitivity level close to 0. This choice does not impose transitivity of pairs as the linearly transitive cluster ${\mathcal{J}}_T$ may be empty, except for all pairs with object 1 which are classified as transitive due to our imposition of constraints for identifiability from Proposition 1. Let the distinct set of intransivity levels be$${\theta }_{\mathcal{K}}:=\{{\theta }_0=0,\theta =\{{\theta }_k\in {\mathbb{R}}_{+},\forall k\in \mathcal{K}\}:0<{\theta }_1<\dots <{\theta }_K\},$$where $\mathcal{K}=\{1,\dots ,K\}$ and the levels of intransitivity are ordered from smallest to largest. The levels of intransitivity, $\theta$, contain the set of positive values of intransitivity which, due to symmetry and the completely transitive cluster with intransitivity value ${\theta }_0$, then define the full $2K+1$ possible values of intransitivity between any pair of objects.

We define the intransitivity cluster allocation of a given pair $\{i,k\}$ to be another binary matrix $z_{\{i,k\}}$, which takes the value 1 at position $s\in \{-K,\dots ,K\}$ and 0 everywhere else, if the pair $\{i,k\}$ belongs to cluster s. The clusters are, therefore, labeled from $-K$ to K, where a cluster labeled $k\in \{1,\dots K\}$ has cluster level ${\theta }_k$, a cluster labeled $k\in \{-K,\dots ,-1\}$ has cluster level $-{\theta }_{-k}$, and a cluster with label 0 has cluster level ${\theta }_0=0$. The set $Z:=\{z_{\{i,k\}},\forall i>k\in \mathcal{I}\setminus \{1\left\}\right\}$ then defines all the cluster allocations for all the free pairs $i\ne k\in \mathcal{I}\setminus \left\{1\right\}$, because of the rotational symmetry. For example, if the Kth index of $z_{\{i,k\}}$ has value $z_{\{i,k\},K}=1$, then this indicates that the pair $\{i,k\}$ belongs to the cluster with label K, whose cluster level is the largest level of intransitivity ${\theta }_K$, and this enforces that the pair $\{k,i\}$ belongs to cluster $-K$ and has the smallest level of intransitivity $-{\theta }_K$. If the cluster allocation of the pair $\{i,k\}\in \mathcal{I}\setminus \left\{1\right\}$ is(8) $${\mathcal{S}}_{\{i,k\}}(Z):=\{\begin{array}{cc}{\text{argmax}}_s{z}_{\{i,k\},s}& \text{if}\mathrm{ }i>k,\\ -{\text{argmax}}_s{z}_{\{k,i\},s}& \text{if}\mathrm{ }i<k,\end{array}$$(8) then the level of intransitivity for a pair $\{i,k\}$, ${\theta }_{ik}$ can be redefined as(9) $$\begin{array}{cc}{\theta }_{ik}\hfill & =\{\begin{array}{cc}{\theta }_{{\mathcal{S}}_{\{i,k\}}(Z)}\mathrm{ }\mathbb{1}\{{\mathcal{S}}_{\{i,k\}}(Z)\ge 0\}& \\ -{\theta }_{-{\mathcal{S}}_{\{i,k\}}(Z)}\mathbb{1}\{{\mathcal{S}}_{\{i,k\}}(Z)<0\},& \{i,k\}\in \mathcal{I}\setminus \left\{1\right\}\\ 0,& \text{otherwise}\end{array}\hfill \\ :=f_{\theta }\left(\theta ,Z,\{i,k\}\right)\hfill & \hfill \end{array}$$(9) where $\mathbb{1}$ is the indicator function, and remembering that ${\theta }_0=0$.

The full model can be written either in terms of Equation (6)(6) $$p_{ik}=\frac{1}{1+\hspace{0.17em}\exp [-\left({\theta }_{ik}+r_i-r_k\right)]},\forall i\ne k\in \mathcal{I}.$$(6) , noting that the parameters will be clustered, or can be written in terms of the levels and the cluster allocations,(10) $$p_{ik}={(1+\hspace{0.17em}\exp \hspace{0.17em}\left\{-\left[f_{\theta }\left(\theta ,Z,\{i,k\}\right)+f_r\left(\varphi ,Y,i\right)-f_r\left(\varphi ,Y,k\right)\right]\right\})}^{-1}.$$(10)

So the ICBT model is defined by $\psi =\{\varphi =\{{\varphi }_a:a\in \mathcal{A}\setminus \left\{0\right\}\},$ $\theta =\{{\theta }_k:k\in \mathcal{K}\}\}$.

Due to the intransitivity levels being fixed to 0 for all pairs of objects involving object 1, an adjustment is required to get a more interpretable value of intransitivity between the pairs. We define the adjusted intransitivity to be(11) $${\theta }_{ij}^{*}:=\text{logit}\left(p_{ij}\right)-\text{logit}\left(p_{ij}^{(\text{BT})}\right)={\theta }_{ij}+r_i-r_j-\left(r_i^{(\text{BT})}-r_j^{(\text{BT})}\right),$$(11) that is, the difference between the logits of the pairwise probability between our ICBT model and the Bradley-Terry model. Note that the rotational symmetry of $\left\{{\theta }_{ij}\right\}$ (4) also imposes rotational symmetry on $\left\{{\theta }_{ij}^{*}\right\}$, that is, ${\theta }_{ij}^{*}={\theta }_{ji}^{*},\forall i\ne j\in \mathcal{I}$.

To help see the value of this reparametrisation, consider then the earlier example of a deterministic game of Rock-Paper-Scissors. Take Rock as the constrained object, then Rock has fixed skill level $r_r=0$, and that pairs involving Rock have intransitivity 0, that is ${\theta }_{rp}={\theta }_{rs}=0$, where the p and s subscripts denote Paper and Scissors. To maintain that Rock always beats Scissors $p_{rs}=1$, then from the constraints, we get an excellent approximation from the ICBT model when $r_s=-M$ for some large M, with the approximation improving as $M\to \infty$. Likewise $r_p=M$, and ${\theta }_{ps}=-3M$. With this model there is only one skill level M, and one nonzero intransitivity level $-3M$. This parametrisation somewhat hides the symmetry of the intransitivity over pairs. However, with definition (11), then ${\theta }_{rs}^{*}={\theta }_{sp}^{*}={\theta }_{pr}^{*}=M$, resulting in an intuitive and easy interpretation of the intransitivity, reflecting the symmetry of the game, no-matter the choice of the fixed parameters.

3.3 Model Ranking

In the Bradley-Terry model, the skill parameters can simply be ordered to give a rank since a greater skill always results in higher win probabilities against all other objects. In our ICBT model this is not the case, because both the intransitivity parameters of each pair and the skill parameters of the objects impact the win probability between any pair. However, below we present two intuitive methods for determining overall ability, and therefore, ranking.

First, if $p_{ij}=\mathrm{Pr}\{i\succ j\}$ is the probability of an object i beating object j according to our model, then we can rank the objects by ordering(12) $$p_{.}:=\left\{p_{i.}:=\frac{1}{n-1}\sum _{j\in \mathcal{I}:\mathrm{ }j\ne i}{p}_{ij}:i\in \mathcal{I}\right\},$$(12) that is, $p_{i.}$ is the average probability of object i beating any other object $j\ne i\in \mathcal{I}$.

Second, if we consider the intransitivity between an object i and an opposing object $j\ne i$ as some “boost” which contributes to the overall ability (which could be negative), then the overall ability $a_i$ of object i could be defined by(13) $$a_i:=r_i+\frac{1}{n}\sum _{j\in \mathcal{I}}{\theta }_{ij},\text{\hspace{0.17em}where\hspace{0.17em}}{\theta }_{ii}=0,\forall i\in \mathcal{I},$$(13) that is, the object skill plus its average intransitivity level. Definition (13) is equivalent to the Bradley-Terry definition of “ability.” Defining $\text{logit}(p_{ii}^{(\text{BT})})=0,\forall i\in \mathcal{I}$, a Bradley-Terry gives(14) $$\frac{1}{n}\sum _{j\in \mathcal{I}}\text{logit}\left(p_{ij}^{(\text{BT})}\right)=r_i-\frac{1}{n}\sum _{j\in \mathcal{I}}{r}_j,$$(14) where the sum on the right hand side does not depend on i, so the skill of object i is entirely determined by $r_i$. Similarly, in our model(15) $$\frac{1}{n}\sum _{j\in \mathcal{I}}\text{logit}\left(p_{ij}\right)=r_i+\frac{1}{n}\sum _{j\in \mathcal{I}}{\theta }_{ij}-\frac{1}{n}\sum _{j\in \mathcal{I}}{r}_j=a_i-\frac{1}{n}\sum _{j\in \mathcal{I}}{r}_j,$$(15) then given definition (13), both (14) and (15) have the same form but with $a_i$ replacing $r_i$. Then a ranking can be formed by ordering the set of abilities $a:=\{a_i:i\in \mathcal{I}\}$. We argue that the first method, using the probabilities $p_{.}$ to rank the objects, is more meaningful since it is directly associated with the pairwise probabilities, the modeling of which is our ultimate aim. The application to baseball data of both methods is discussed in the supplementary material.

4 Inference

4.1 Likelihood

The data, $x:=\{x_c:c\in \mathcal{C}\}$, are binary, and $i\succ j$ denotes that i is preferred to j. Then $x_c=1$ if $i_c\succ j_c$, and $x_c=0$ otherwise, where $i_c,j_c\in \mathcal{I}$ are the objects being compared in comparison c. Then, the log-likelihood for the ICBT model is(16) $$\mathcal{l}(x|\varphi ,Y,A,\theta ,Z,K)=\sum _{c\in \mathcal{C}}\left[x_c\hspace{0.17em}\log \hspace{0.17em}\left(p_{i_c{j}_c}\right)+(1-x_c)\hspace{0.17em}\log \hspace{0.17em}\left(1-p_{i_c{j}_c}\right)\right],$$(16) where $p_{i_c{j}_c}$ is given by the ICBT model for all $c\in \mathcal{C}$ and is calculated from the set of parameters $(\varphi ,Y,A,\theta ,Z,K)$. All pairs’ intransitivities $\{{\theta }_{ik}:i\ne j\in \mathcal{I}\}$ can be found from the intransitivity levels $\theta$ and the cluster allocations $Z$, using equation (3.2), so it is only necessary to do inference on these parameters, rather than the full $2K+1$ separate clusters. Therefore, from here onwards the term intransitivity levels refers only to those K values which have positive intransitivity. Similarly, any individual object’s skill $r_i\forall i\in \mathcal{I}$ can be found from knowing the ability levels $\varphi$ and the cluster allocations $Y$, using Equation (7)(7) $$r_i=\{\begin{array}{cc}{\varphi }_{{\mathcal{S}}_{\left\{i\right\}}\left(Y\right)},& i\in \mathcal{I}\setminus \left\{1\right\}\\ 0,& i=1\end{array}\hspace{1em}:=f_r\left(\varphi ,Y,i\right),\forall i\in \mathcal{I}.$$(7) . We formulate a Bayesian hierarchical model, which treats both K and A as unknown parameters, thus, accounting for uncertainty in the number of clusters. The posterior is, therefore, written as$$\pi \left(\varphi ,Y,A,\theta ,Z,K|x\right)\propto L\left(x|\varphi ,Y,A,\theta ,Z,K\right)\pi \left(\varphi ,Y,A,\theta ,Z,K\right)$$ where $L(\cdot )=\hspace{0.17em}\exp [\mathcal{l}(\cdot )]$ is the likelihood and $\pi (\varphi ,Y,A,\theta ,Z,K)$ is the prior.

4.2 Prior Specification

Formulating the prior, we make the assumption that $Z\perp \perp$ $\theta |K$ that is, the intransitivity level allocations and intransitivity levels are independent from one another given the number of intransitivity levels K. Likewise, it is assumed that $Y\perp \perp$\$\varphi |A$. Furthermore, we assume that the clustering of the objects’ skills and the clustering of the pairs’ intransitivities are independent systems, that is, $A\perp \perp$ $K$, $\varphi \perp \perp$ $\theta$, and $Y\perp \perp$ $Z$. This means that the prior specification for the two features we are clustering, skills and intransitivities, can be approached separately.

Consider first the prior specification for the clustering of the intransitivity values of the pairs. Remember that labels $z_{\{i,j\}},\forall i>j\in \mathcal{I}\setminus \left\{1\right\}$ have domain $\{-K,\dots ,K\}$, that is, $z_{\{i,j\}}:\{-K,\dots ,K\}\to \{0,1\},\forall i>j\in \mathcal{I}\setminus \left\{1\right\}$, and also that $z_{\{i,1\}},\forall i\in \mathcal{I}\setminus \left\{1\right\}$ (and by symmetry $z_{\{1,j\}},\forall j\in \mathcal{I}\setminus \left\{1\right\}$ too) are fixed in the transitive level $\left\{0\right\}$ for identifiability purposes, see Section 3.2. Let the prior on the cluster allocation be$$z_{\{i,j\}}|{\omega }_K\sim \text{multinomial}\left(1,{\omega }_K\right),\forall i>j\in \mathcal{I}\setminus \{1\},$$where ${\omega }_K$ is on $\{-K,\dots ,K\}$ such that$${\omega }_K=\{{\omega }_{K,s}\in [0,1]:s\in \{-K,\dots ,K\},\sum _{s=-K}^K{\omega }_{K,s}=1\}.$$

The distribution of $Z|({\omega }_K,K)$ is assumed independent over all pairs $i>j\in \mathcal{I}\setminus \left\{1\right\}$, that is,$$f\left(Z|{\omega }_K,K\right)=\frac{\left(\sum _{k=-K}^K\left|b_k\right|\right)!}{\prod _{k=-K}^K\left|b_k\right|!}\prod _{s=-K}^K{\omega }_{K,s}^{\left|b_s\right|},$$where $b_k=\{(i,j):i>j\in \mathcal{I}\setminus \{1\}:z_{\{i,j\},k}=1\}\forall k\in \{-K,\dots ,K\}$ is the set of allocated pairs of objects belonging to cluster k. We set ${\omega }_K|K\sim \text{Dirichlet}({\overline{\gamma }}_K)$ to come from $2K+1$ dimensional Dirichlet prior distribution, and ${\overline{\gamma }}_K\in {\mathbb{R}}_{+}^{2K+1}$ is the hyperparameters vector. We use an uninformative prior, setting ${\overline{\gamma }}_K={\gamma }_K{1}_{2K+1}$ where $1_{2K+1}$ is a vector of ones of length $2K+1$ and ${\gamma }_K\in {\mathbb{R}}_{+}$. In this case, the ${\omega }_K$ parameter can be marginalized out, by(17) $$\begin{array}{cc}f(Z|{\gamma }_K,K)=\hfill & {\int }_{{\omega }_K}f(Z|{\omega }_K,K)f({\omega }_K|{\gamma }_K,K)\mathrm{}\mathrm{d}{\omega }_K\hfill \\ =\frac{\left(\sum _{k=-K}^K\left|b_k\right|\right)!}{\prod _{k=-K}^K\left|b_k\right|!}\frac{\Gamma ((2K+1){\gamma }_K)}{\Gamma {({\gamma }_K)}^{2K+1}}\hfill & \hfill \\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\times \frac{\prod _{k=-K}^K\Gamma \left({\gamma }_K+\left|b_k\right|\right)}{\Gamma \left((2K+1){\gamma }_K+\sum _{k=-K}^K\left|b_k\right|\right)}\hfill & \hfill \end{array}$$(17) where integration on ${\omega }_K$ is taken over the $2K+1$ simplex. This is referred to as a Dirichlet-multinomial allocation prior. The prior for K is a Poisson $({\lambda }_K)$ distribution with probability mass function denoted $g_0(k|{\lambda }_K)$ so that $\mathbb{E}\left[K|{\lambda }_K\right]={\lambda }_K$, with ${\lambda }_K>0$. Note that $K=0$ is feasible, as this corresponds to the Bradley-Terry model since $\theta |(K=0)=\varnothing$ and so only the transitive cluster exists, that is ${\theta }_{\mathcal{K}}|(K=0)={\theta }_0$, and $\{i,j\}\in {\mathcal{J}}_T,\forall i\ne j\in \mathcal{I}$, so all pairs belong to the transitive cluster ${\mathcal{J}}_T$. Formally $K<n(n-1)/2-n$ but as this is large relative to our prior beliefs on K, for simplicity we ignore this constraint in the inference.

As the $\theta$ elements are ordered in increasing order and are positive, the prior on the $\theta$ parameters is taken to be the joint distribution of K order statistics drawn from independent gamma random variables, such that(18) $$\begin{array}{c}{h}_0\left(\theta |K\right)=K!\prod _{i=1}^K{h}_0\left({\theta }_i|\alpha ,\beta \right),\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{with}\mathrm{ }\hspace{1em}0<{\theta }_1<\dots <{\theta }_K,K\ge 1,\end{array}$$(18) and where $h_0(x|\alpha ,\beta )$ is the Gamma $(\alpha ,\beta )$ density with shape and scale $\alpha ,\beta >0$, respectively.

Consider the prior for the skill levels clustering. The set of skill cluster allocations has distribution $Y=\left\{y_{\left\{i\right\}}\right|{\omega }_A,A\sim \text{multinomial}(1,{\omega }_A),\forall i\in \{2,\dots ,n\}\},$ where ${\omega }_A$ has domain on $\{-A_{-},\dots ,A_{+}\}$. The distribution of $y_{\left\{i\right\}}|({\omega }_A,A)$ is assumed independent over all objects $i\in \mathcal{I}\setminus \left\{1\right\}$ such that$$f(Y|{\omega }_A,A)=\prod _{i\in \mathcal{I}\setminus \left\{1\right\}}f(y_{\left\{i\right\}}|{\omega }_A,A),$$where$${\omega }_A=\{{\omega }_{A,s}\in [0,1]:s\in \{-A_{-},\dots ,A_{+}\},\sum _{s=-A_{-}}^{A_{+}}{\omega }_{A,s}=1\}.$$

Again, ${\omega }_A|(A,{\gamma }_A)\sim \text{Dirichlet}({\overline{\gamma }}_A)$ is modeled to come from an $A+1$ dimensional Dirichlet prior distribution, with ${\overline{\gamma }}_A={\gamma }_A{1}_{A+1}$ where ${\gamma }_A\in {\mathbb{R}}_{+}$. Marginalizing out as in derivation (4.2), another Dirichlet-multinomial allocation prior is obtained by integrating ${\omega }_A$ over the $A+1$ dimensional simplex. The prior density for the skill allocations is, therefore, given as(19) $$f(Y|{\gamma }_A,A)=\frac{n!}{\prod _{a=-A_{-}}^{A_{+}}\left|c_a\right|!}\frac{\Gamma (A{\gamma }_A)}{\Gamma {({\gamma }_A)}^A}\frac{\prod _{a=-A_{-}}^{A_{+}}\Gamma ({\gamma }_A+|c_a|)}{\Gamma (A{\gamma }_A+n)},$$(19) because ${\sum }_{a=-A_{-}}^{A_{+}}\left|c_a\right|=n$, where $c_a:=\{i,\forall i\in \mathcal{I}\setminus \{1\}:y_{\left\{i\right\},a}=1\}$ is the set of objects belonging to skill cluster $a\in \mathcal{A}$.

The prior for the number of unknown skill levels A is taken to be a truncated Poisson distribution with parameter $({\lambda }_A)$, ${\lambda }_A>0$ with probability mass function$$g_A\left(a|{\lambda }_A\right)=\frac{{\lambda }_A^a}{a!}{\left(\sum _{i=0}^{n-1}\frac{{\lambda }_A^i}{i!}\right)}^{-1}a=0,1,\dots ,n-1.$$

Similarly to $\theta$, the prior choice for $\varphi$ is taken to be the joint distribution of order statistics of independent and identically distributed $A+1$ Gaussian random variables such that$$\pi (\varphi |A)=(A+1)!\mathrm{}\prod _{a\in \mathcal{A}\setminus \left\{0\right\}}\pi ({\varphi }_a)\mathrm{ }\text{for}\mathrm{ }{\varphi }_{A_{-}}<\dots <{\varphi }_0<\dots {\varphi }_{A_{+}},$$where ${\varphi }_a\sim \mathcal{N}(0,{\nu }_A^2)\forall a\in \mathcal{A}\setminus \left\{0\right\}$, and with ${\nu }_A\in {\mathbb{R}}_{+}$. The $(A+1)!$ term arises as ${\varphi }_0$ can occur anywhere in the sequence of $\varphi$.

In summary, the prior $\pi (\varphi ,Y,A,\theta ,Z,K)$ is equal to(20) $$\begin{array}{c}(A+1)!\left[\prod _{a\in {\mathcal{A}}_{\{-0\}}}\pi \left({\varphi }_a\right)\right]f(Y|{\gamma }_A,A)g_A\left(A|{\lambda }_A\right)K!\\ \hspace{1em}\times \left[\prod _{i=1}^K{h}_0\left({\theta }_i|\alpha ,\beta \right)\right]f(Z|{\gamma }_K,K)g_0\left(K|{\lambda }_K\right),\end{array}$$(20) where, ${\lambda }_K,{\lambda }_A,{\gamma }_K,{\gamma }_A,{\nu }_A$, and $\alpha ,\beta$ are the hyperparameters.

4.3 Reversible Jump Markov Chain Monte Carlo Sampler

Inference is made via a reversible jump Markov chain Monte Carlo sampler (Green 1995), which provides samples from the posterior distribution $\pi (\varphi ,Y,A,\theta ,Z,K|x)$, that is, the intransitivity and skill levels, the allocations to the these levels, and the number of levels. Since the number of skill and intransitivity levels (A, K) are assumed to be unknown, the uncertainty in these parameters must be accounted for, thus, motivating the use of a reversible jump sampler. In a sense, the reversible jump sampler mixes over models as well as parameters, and thus, fully accounts for this uncertainty in the final inference.

The ICBT model is structured to try to favor the Bradley-Terry model as a special case, and this is reflected in our sampler, by explicitly incorporating the completely transitive cluster ${\theta }_0$ as an ever present cluster, even if no pairs are allocated to this cluster at a given iteration of the sampler. To ensure the skill and intransitivity levels both remain ordered, the updates to these levels occur in a transformed space such that no update can lead to a change in order.

The reversible jump algorithm used is a split-merge sampler (Green and Richardson 2001), which is adapted from the work of Ludkin (2020). The sampler comprises three separate moves: a standard Markov chain Monte Carlo Metropolis-Hastings move, which samples parameters $\varphi$, $\theta$, and reallocates clusters $Y,Z$; splitting or merging clusters; and adding or deleting empty clusters. For the construction of the algorithm and how it is implemented, see the supplementary material.

4.4 Model Assessment

The inferences produced by any model are only meaningful if the model itself is accurate. This accuracy is measured here by how well the model fits out of sample. If $\mathcal{C}$ is the set of total observed pairwise comparisons, then let ${\mathcal{C}}_s$ be the set of comparisons on which the model is fitted and ${\mathcal{C}}_t$ be the set of comparisons on which the model performance is analyzed, such that ${\mathcal{C}}_s\cup {\mathcal{C}}_t=\mathcal{C}$ and ${\mathcal{C}}_s\cap {\mathcal{C}}_t=\varnothing$. We use log-loss $l(x^{*})$ of the test dataset $x^{*}$ to measure model performance, which we take to be the average negative log-likelihood per observation in $x^{*}$, that is,(21) $$l(x^{*})=-\frac{1}{\left|{\mathcal{C}}_t\right|}\sum _{c\in {\mathcal{C}}_t}\left[x_c\hspace{0.17em}\log \hspace{0.17em}({\widehat{p}}_{i_c{j}_c})+(1-x_c)\hspace{0.17em}\log \hspace{0.17em}(1-{\widehat{p}}_{i_c{j}_c})\right],$$(21) where ${\widehat{p}}_{ij}$ is the point estimate of $p_{ij}$ based on the training dataset of comparisons ${\mathcal{C}}_s$ and $x^{*}:=\{x_c:c\in {\mathcal{C}}_t\}$ is the set of test data, where the notation is as used in the expression for the likelihood (16) .

4.5 Simulation Study

The model was tested using simulated datasets where the number of objects, the number of round-robin tournaments, and the amount of intransitivity varied between the datasets. The sensitivity of our model to these parameters was then tested by comparing out of sample prediction accuracy with a standard Bradley-Terry model. This provided insights into the amount of data, and the amount of intransitivity, required for our more complex model to outperform the Bradley-Terry model. A full analysis is provided in the supplementary material.

5 Baseball Data

5.1 Data

Baseball was chosen to illustrate the methodology due to the high frequency of games, with accessible data for the American League Baseball obtained from www.retrosheet.org. The data are from the 2010–2018 seasons, with the 2010–2012 seasons involving 14 teams, and the 2013–2018 seasons involving 15 teams due to the Houston Astros moving from the National League to the American League. We analyze each season’s data separately here, and jointly over years in the supplementary material.

The tournament structure is not as simple as the round robin tournament we considered in the simulation study. The American League is split into three divisions based on location: East, Central and West, with five teams in each (since 2013). Within the same division, pairs of teams play each other approximately 20 times, and pairs of teams from different divisions play each other around 5–7 times, as well as any Playoffs and World Series matchups, totaling around 140–160 matches per team every season, depending on the season and the team. Baseball is known to be a highly strategic game, with issues such as player selection, handedness of the of batters, strength and speed of players, and tactics such as “small ball” versus “long ball” all considered of great importance. So we anticipate that the level of intransitivity will be high.

The vast majority (at least $99.5\%$) of all matches are played at the home of one of the two teams competing in the game, with the rest played at neutral venues. Playing at home is well known to have the potential to increase the probability of the home team winning the match across a range of sports (Dixon and Coles 1997). Although prediction and model interpretation could be improved by incorporating this effect, we decided not to address home advantage here. Our reason was that none of the existing intransitivity models have such a feature, as they were developed for applications devoid of home advantage, such as e-sports, so a comparison of the different models would only be fair if we did not include this property. However, in Section 6 we formulate the home advantage adaptation given its potential interest.

If pairs of teams do not play equally home and away, then ignoring home advantage could lead to misinterpretation of the estimated ICBT model parameters, for example, if team i mostly played team k with team i at home, the home advantage would feed into ${\theta }_{ik}$. We do not believe this is problematic due to the near perfect balance of home to away matches per team, and the maximum home percentage within pairs of teams is 70% for 2010–2011 and only 57% subsequently.

5.2 Inference

The baseball data are analyzed using the ICBT model, and its results are compared with those of the Bradley-Terry model and with the existing models of Section 2 except for the model of Duan et al. (2017) due to the subjective choices required for some parameters.

The ICBT model incorporates uncertainty in the choice of model itself, that is, the number of clusters, and therefore, how many parameters. Our prior distributions for number of intransitivity levels K and skill levels A are shown in Figure 1. We used a Poisson(${\lambda }_K=2$) prior for K, with the hyperparameter to give a 95% prior chance that $K\in [0,5]$, as it was thought that there would only be a few different pairwise strategies. Similarly, the prior for A was taken to be Poisson(λ_A = 7) as this hyper-parameter choice gave a 95% prior chance that $A\in [2,13]$. The justification for our choice of the other hyper-parameters $({\gamma }_K,{\gamma }_A,\alpha ,\beta ,{\nu }_A$) and a sensitivity analysis to hyper-parameter choice is reported in the supplementary material.

Figure 1: Posterior distributions of the K intransitivity levels (left) and the A skill levels (right) for the 2018 season: with the associated prior distributions in a lighter color.

Now consider the posterior distributions for K and A based on the 2018 season data, also shown in Figure 1. Despite the prior only providing vague information across values of $K\le 5$, the data clearly favors having a single intransitivity level, meaning three possible clusters for each pair: a positive level, the completely transitive level, and the mirrored negative level. Further, although the prior gave a 14% probability to the Bradley-Terry model $(K=0)$, the posterior probability for that model is estimated to be zero, showing strong evidence of intransitivity in the dataset. For the distinct skill levels the change from prior to posterior is relatively small, with a mean (and 95% credible interval) of $7.94(4,12)$, with the number of distinct skills levels favoring $A\in [6,9]$.

The posterior estimated values of the teams’ skills, and the variance of these values in particular, provides a helpful summary of how competitive the tournament is in a season, with the smaller the variance the more closely contested the tournament. For our model the skills’ variance ranged from 0.027 (2012 season) to 0.32 (2018 season) and for the standard Bradley-Terry model, 0.031 (2015 season) to 0.23 (2018 season) - both models suggesting that 2018 was the least competitive season. In the 2013 season the variance of skill levels according to Bradley-Terry is almost 3 times that of our model. However, the 2013 season was found to contain a particularly large amount of intransitivity, indicating that the large range of skills in the Bradley-Terry model could be the result of compensating for an inability to express the intransitivity. This has perhaps resulted in the Bradley-Terry model concluding that the 2013 season was less competitive than it was in reality.

Now consider the pairwise interactions between teams. These interactions could be inferred from either the intransitivity of the posterior mean ${\widehat{\theta }}_{ij}^{*},\forall i\ne j$, or by the posterior mean of the intransitivity parameter ${\widehat{\theta }}_{ij},\forall i\ne j$. The supplementary material contains a comparison for both and concludes that ${\widehat{\theta }}_{ij}^{*}$ is more meaningful and interpretable here, so we focus on that. For the 2018 season, Figure 2 (left) shows ${\widehat{\theta }}_{ij}^{*}$, for each pair of teams $i>j\in \mathcal{I}$: recall that intransitivity has rotational symmetry, that is, ${\theta }_{ij}^{*}=-{\theta }_{ji}^{*},\forall i\ne j$. The teams are sorted by their rank according to $p_{.}$, given by definition (12), see Figure 2 (right). Reading from the teams on the y-axis to x-axis there is a large positive value of intransitivity from Baltimore (BAL) to Tampa Bay (TBA) of 0.78 with 95% credible interval $(0.36,1.22)$, indicating that Baltimore played better against Tampa Bay than expected, given their overall abilities. This is consistent with the data, with Baltimore winning 11 out of 19 matches between the two teams, despite being ranked lower.

Figure 2: Analysis of 2018 season: the posterior mean of the intransitivity parameter, ${\widehat{\theta }}_{ij}^{*}$ across all pairs of teams $i>j\in \mathcal{I}$ (left); ranking according to definition (12) (black) and Bradley-Terry model (red) for all teams $i\in \mathcal{I}$ (right).

The analysis of these intransitivities between pairs, and that of the skills of each team, can be combined to produce an overall ranking of the teams. As discussed in Section 3.3, with further details in the supplementary material, $p_{.}$ provides a suitable ranking of the teams. For the 2018 season Figure 2 (right) shows the ranks according to $p_{.}$ compared to the Bradley-Terry ranks. Both have been linearly scaled to help with a visual comparison, such that the best and worst teams have abilities 1 and 0, respectively. The two sets of estimated rankings using $p_{.}$ are clearly correlated; however, there is some difference in the ordering of the ranks, indicating that intransitivity may have been masking the true ranks of some teams. For example, consider Tampa Bay, ranked 6th by the Bradley-Terry model. Tampa Bay’s good record against Kansas City (KCA) has a much lower weighting than their poor record against Baltimore in the Bradley-Terry model due to the differing frequency of these match-ups, and therefore, impacts the overall rank of Tampa Bay. The ICBT model, however, recognizes that good or bad records against particular teams could be due to the presence of intransitivity, and therefore, penalizes Tampa Bay less overall, ranking them 5th, thus, illustrating our point in Section 1 that the ICBT model makes adjustments for tournament imbalance. Similar plots and inferences are drawn from the other seasons (2010–2017) but with different team rankings in each year.

5.3 Model Performance

To test the model performance, 70% of games from each season are randomly selected to be training data, on which the model is fitted, with the remaining 30% used as test data, on which the log-loss score is calculated. This random selection is appropriate as none of the models compared take time-dependency into consideration, a feature discussed in the supplementary material. The variation due to this random sampling in the training-test split is accounted for by taking 100 separate random training-test splits for each season. For each replicate of training data the model is fitted separately to each season’s data. Relative log-loss is then calculated by subtracting the log-loss of a baseline coin tossing model.

Table 1 shows these negative relative log-loss scores for all years of baseball data, along with 95% confidence intervals, with this measure evaluated for the ICBT, Bradley-Terry, blade-chest and majority vote models. Since a larger value of negative relative log-loss indicates better model performance, a positive value indicates an improvement on the coin tossing model. So all four models improve on simply using coin tossing, showing that there is information to be exploited for inference and prediction. The Bradley-Terry, blade chest and majority vote models all have somewhat similar performance to each other across the years, with the most improved fit being in 2018. In contrast our model is the best performing out of the four models in terms of out-of-sample prediction on all years of data. When assessed as the cumulative improvement over years, relative to coin tossing, the ICBT model is 2.8 times better than the Bradley-Terry model, showing that we have substantially improved predictive performance. The difference in log-loss scores relative to the Bradley-Terry model is largest for the 2013 season, which in Section 5.2 has been identified as the season with the largest intransitivity.

Table 1: Negative relative log-loss $\times {10}^3$ (compared to a coin-tossing model) for each year of baseball data for the ICBT, Bradley-Terry, blade-chest and majority vote models.

Year	ICBT	BT	Blade-chest	Majority vote
2010	44(38, 46)	17(15, 18)	17(2, 27)	20(13, 25)
2011	46(33, 49)	15(13, 17)	17(-1, 26)	20(13, 26)
2012	49(44, 53)	14(11, 16)	24(8, 33)	22(14, 31)
2013	64(36, 69)	23(21, 25)	33(13, 42)	31(22, 39)
2014	39(29, 45)	9(7, 11)	10(-12, 21)	13(6, 19)
2015	34(12, 45)	5(2, 7)	9(-14, 17)	9(1, 16)
2016	42(32, 55)	10(8, 12)	18(2, 30)	18(10, 28)
2017	36(10, 50)	13(11, 15)	13(-4, 22)	16(9, 22)
2018	73(65, 79)	46(44, 47)	48(19, 56)	51(43, 57)

95% confidence intervals, in parentheses, come from random training-test splits of the data.

6 Conclusions and Discussion

We have proposed a new model and inference structure for paired comparison data. We frame this in the context of sport competitions, baseball in particular, with teams competing against each other, though the potential applications of the model are much broader. Our proposed model, the Intransitive Clustered Bradley-Terry (ICBT) model, extends the standard Bradley Terry model, which is widely considered as the baseline model for such data. The extension allows for intransitivity so that the difference in skill levels between two objects being compared is not the only factor affecting the probabilities of the outcomes. There are a number of models which already allow for intransitivity, but each of these are quite restricted in the parametric form of intransitivity relative to our semi-parametric approach, which recognizes that certain patterns of interaction between pairs of objects can be common over multiple pairs. Our model also allows for objects’ skills to be clustered, a feature that is novel to paired comparisons, with this inducing parsimony and avoiding obtaining distinct rankings for some items when there is no evidence from the data that they are not equally good. We have shown evidence from American League baseball that our model provides a distinct improvement on existing models.

The ICBT model has complete flexibility, in the sense that cluster allocation to skill and intransitivity levels is not predetermined. In order that the data identify the appropriate structure of clustering, and for the inference to account for the uncertainty in this choice, the model is fitted via RJMCMC.

Based on the clusters with the highest posterior probabilities, we anticipate that experts in the particular sport may be able to identify some patterns of clustering that are interpretable, for example, associated with different styles of play. In such cases, these clustering features could be hard wired into the model as the only options, resulting in more efficient inference. A referee made the helpful suggestion that if accounting for clustering uncertainty was not an issue then the inference could be simplified by estimating the ICBT model with group lasso penalties to induce clusters. We feel that our model works sufficiently well for the current applications but agree that it presents an exciting springboard for the consideration of various extensions to the model and its inference. We finish by illustrating a few such possible extensions.

In Section 5 we did not attempt to account for home advantage, which is widely recognized as an important feature in sport, for example, Cattelan, Varin, and Firth (2013) incorporate it in a Bradley-Terry model, though to the best of our knowledge it has not been accounted for in the existing intransitivity models. The most natural way to achieve this is to change $p_{ik}$ given by expression (6) to a probability $p_{ik}^{(i)}$ of the home team i beating the away team k, with(22) $$\begin{array}{c}{p}_{ik}^{(i)}=\frac{1}{1+\hspace{0.17em}\exp [-\left({\theta }_{ik}+\gamma +r_i-r_k\right)]},\\ \text{\hspace{1em}\hspace{1em} and}{p}_{ki}^{(i)}=1-p_{ik}^{(i)}\forall i\ne k\in \mathcal{I},\end{array}$$(22) where $\gamma \in \mathbb{R}$ determines the effect of playing at home, which here is common over all pairs of teams. If $\gamma >0$ ($\gamma <0$) then the probability of a home win is increased (decreased) relative to the other factors of skill and intransitivity. This effect can be extended to vary over teams by replacing $\gamma$ by ${\gamma }_i$ in expression (6). To ensure these ${\gamma }_i$ parameters are all identifiable, we fix ${\gamma }_1=0$, though no additional constraints are needed if there is a common $\gamma$, but that is all that is required under the conditions of Proposition 1 on the other parameters, as we are able to exploit data that distinguishes which team is at home.

This article only considered win-loss scenarios. Extensions of the Bradley-Terry have been proposed for handling draws. Two distinct methods for handling draws are given by Cattelan, Varin, and Firth (2013) and Hankin (2020). The former use ordinal logistic regression, treating win, loss, draw as outcomes of an ordered multinomial random variable, which can then be analyzed via an ordered link model. In contrast, the latter treats the problem as a competition between the two teams and a third theoretical team, such that when the theoretical team wins the outcome of the match corresponds to a draw between the two actual teams. The ICBT model can be adapted similarly, with the use of the clustering strategy extended to pooling teams to account for their similar cautiousness, leading to them drawing more often than would be expected.

We have assumed that all teams play each other. If this is not the case we cannot improve on the prior inference for the ${\theta }_{ik}$ parameters for pairs $(i,k)$ that do not play each other. This is not a restriction for Bradley-Terry or the existing intransitivity models, where the associated $p_{ik}$ are determined by the observed pairs. This raises issues about identifiability of the ICBT model parameters. Our approach, through Proposition 1, is no longer sufficient leaving the open problem of which parameters to fix in order to give the most efficient inference.

Acknowledgments

We thank Dr. Matthew Ludkin for giving detailed insight into the reversible jump algorithm. Spearing gratefully acknowledges funding of the EPSRC funded STOR-i Centre for Doctoral Training (grant number EP/L015692/1), and ATASS Sports. We thank the referees for their helpful comments.

Supplementary Materials

Additional Analyses:

The supplementary material for this article include extra analyses: a simulation study; and an extended analysis of the baseball application, including a pooled analysis of data from all seasons. (ICBT_supplementary.pdf)

Reversible Jump Markov Chain Monte Carlo Sampler:

The supplementary material includes full details of the RJMCMC algorithm, including the initialization procedure and highlighting the main novelties of the algorithm. (ICBT_supplementary.pdf)

R Code:

The analyses can be replicated using the code and data found in the supplementary material. See the README contained in the zip file for more details. (spearing_ICBTcode.zip, zip archive)

Appendix:

The supplementary material contains the Appendix, which provides the proof of Proposition 1. (ICBTappendix.pdf)

Engineering and Physical Sciences Research Council;

Disclosure Statement

The authors report there are no competing interests to declare.