Efficient Bayesian Inference for Generalized Bradley-Terry Models

Introduction

Consider a set of $K$ elements. These elements are repeatedly compared with one another in pairs. For two elements $i$ and $j$ of this set, Bradley and Terry, (1952) suggested the following model

where $\lambda_{l}>0$ is a parameter associated to element $l\in\left\{1,2,\ldots,K\right\}$ that represents its skill rating and we denote $\lambda:=\left\{\lambda_{i}\right\}_{i=1}^{K}$.

This model has found numerous applications. As mentioned in (Hunter, , 2004), as early as 1976, a published bibliography on paired comparisons includes several hundred entries (Davidson and Farquhar, , 1976). For example, it has been adopted by the World Chess Federation and the European Go Federation to rank players and it is a standard approach to build multiclass classifiers based on the output of binary classifiers (Hastie and Tibshirani, , 1998). Various extensions have been proposed to handle home advantage (Agresti, , 1990), draws (Rao and Kupper, , 1967), multiple (Plackett, , 1975; Luce, , 1959) and team comparisons (Huang et al., , 2006). In particular, the popular extension to multiple comparisons, named the Plackett-Luce model (Plackett, , 1975; Luce, , 1959), defines a prior distribution over permutations and has been used for ranking of multiple individuals and for choice models (Luce, , 1977). The monographs of David, (1988) and Diaconis, (1988, Chap. 9) provide detailed discussions on the statistical foundations of these models.

For the basic Bradley-Terry model (1), it is possible to find the maximum likelihood (ML) estimate of the skill ratings $\lambda$ using a simple iterative procedure (Zermelo, , 1929; Hunter, , 2004). Lange et al., (2000) established that this procedure is a specific case of the general class of algorithms referred to as MM algorithms. Generally speaking, MM algorithms use surrogate minimizing functions of the log-likelihood to define an iterative procedure converging to a local maximum. EM algorithms are thus just a special case of MM algorithms. An excellent survey of the MM approach and its applications can be found in (Lange et al., , 2000). Hunter, (2004) further derived MM algorithms for generalized Bradley-Terry models and established sufficient conditions under which these algorithms are guaranteed to converge towards the ML estimate.

Recently several authors have proposed to perform Bayesian inference for (generalized) Bradley-Terry models (Adams, , 2005; Gormley and Murphy, , 2009; Görür et al., , 2006; Guiver and Snelson, , 2009). The resulting posterior density is typically not tractable and needs to be approximated. An Expectation-Propagation method is developed in (Guiver and Snelson, , 2009); this yields an approximation of the posterior which can be computed quickly and might be suitable for very large scale applications. However, it relies on a functional approximation of the posterior and the convergence properties of this algorithm are not well-understood. M-H algorithms have been proposed in (Adams, , 2005; Gormley and Murphy, , 2009; Görür et al., , 2006). Gormley and Murphy, (2009) suggested a carefully designed proposal distribution, though it can perform poorly in some scenarios as demonstrated in section 6.

Our contribution here is three-fold. First, we show that by introducing suitable sets of latent variables, the MM algorithms proposed by Hunter, (2004) for the basic Bradley-Terry model and its generalizations to take into account home advantage, ties and multiple comparisons can be reinterpreted as standard EM algorithms. We believe that this non-trivial reinterpretation is potentially fruitful for statisticians who usually like thinking in terms of latent variables. Note that the latent variables introduced here differ from the ones introduced in the standard Thurstonian interpretation of the Bradley-Terry model (Diaconis, , 1988, Chap. 9) and lead to more efficient algorithms as discussed in section 2. Second, using similar ideas, we propose original EM algorithms for some recent generalizations of the Bradley-Terry model including group comparisons and random graphs. Third, based on the sets of latent variables introduced to derive these EM algorithms, we propose Gibbs samplers to perform Bayesian inference in this important class of models. To the best of our knowledge, no Gibbs sampler has ever been proposed in this context. These algorithms have the great advantage of allowing us to bypass the design of proposal distributions for M-H updates and we demonstrate experimentally that they perform very well.

The rest of this paper is organized as follows. In section 2, we consider the basic Bradley-Terry model (1). Based on the introduction of a suitable set of latent variables, we present an EM reinterpretation of the MM algorithm presented by Hunter, (2004) for Maximum a Posteriori (MAP) parameter estimation and an original data augmentation algorithm to sample from the posterior. In section 3, various standard extensions of the Bradley-Terry model allowing for home advantage, ties and competition between teams are described. EM algorithms and original Gibbs sampling schemes are proposed. The Plackett-Luce model (Plackett, , 1975; Luce, , 1959), a very popular generalization of the Bradley-Terry model for multiple comparisons, is presented in section 4. Algorithms applicable to further extensions of the Bradley-Terry model to choice models, random graphs and classification are presented in section 5. In section 6, these algorithms are applied to the NASCAR 2002 dataset and to chess competition data.

Bradley-Terry model

Suppose we have observed a number of statistically independent pairwise comparisons among $K$ individuals. We denote by $D$ the associated data. Let also $w_{ij}$ denote the number of comparisons where $i$ beats $j$ and $n_{ij}=w_{ij}+w_{ji}$ the total number of comparisons between $i$ and $j$. Based on the Bradley-Terry model (1), the log-likelihood function is given by

where the notation $1\leq i\neq j\leq K$ is an abusive notation to denote the set $\left\{\left(i,j\right)\in\left\{1,...,K\right\}^{2}\right.$ $\left.\text{ such that }i\neq j\right\}$ and $1\leq i<j\leq K$ stands for $\left\{\left(i,j\right)\in\left\{1,...,K\right\}^{2}\text{ such that }i<j\right\}$.

We seek to introduce latent variables which are such that the resulting complete log-likelihood admits a simple form. It is well-known that the Bradley-Terry model enjoys the following Thurstonian interpretation (Diaconis, , 1988, Chap. 9): for each pair $1\leq i<j\leq K$ and for each associated pair comparison $k=1,\ldots,n_{ij}$, let $Y_{ki}\sim\mathcal{E}(\lambda_{i})$ and $Y_{kj}\sim\mathcal{E}(\lambda_{j})$ where $\mathcal{E}\left(\varsigma\right)$ is the exponential distribution of rate parameter $\varsigma$ then

These latent variables can be interpreted as arrival times and the individual with the lowest arrival time wins. These latent variables would allow us to define EM and data augmentation algorithms. However, if we only introduce for each pair $i,j$ the sum of the lowest arrival times

instead of $\left\{Y_{ki},Y_{kj}\right\}$ for $k=1,\ldots,n_{ij}$, then the resulting complete log-likelihood remains simple. Here $\mathcal{G}\left(\alpha,\beta\right)$ denote the Gamma distribution of shape $\alpha$ and inverse scale $\beta$. As the fraction of missing information is reduced, this leads to faster rates of convergence for the resulting EM and data augmentation algorithms (Liu, , 2001, Chap. 6). We will essentially proceed similarly to introduce latent variables for generalized Bradley-Terry models.

To summarize, for $1\leq i<j\leq K$ such that $n_{ij}>0$, we introduce the latent variables $Z=\left\{Z_{ij}\right\}$ which are such that

The resulting complete data log-likelihood is given by

If we assign additionally a prior to $\lambda$ such that

as in (Gormley and Murphy, , 2009; Guiver and Snelson, , 2009) then we can maximize the resulting log-posterior using the EM algorithm which proceeds as follows at iteration $t$:

with “$\equiv$” meaning “equal up to terms independent of the first argument of the $Q$ function” and $w_{i}=\sum_{j=1,j\neq i}^{K}w_{ij}$ is the number of wins of element $i$. Using (5), it follows that

For $a=1$ and $b=0$, the MAP and ML estimates coincide. In this case (6) is exactly the minorizing function of the MM algorithm proposed in (Hunter, , 2004, Eq. (10)) and thus the MM algorithm is given by (7).

For $1\leq i<j\leq K$ s.t. $n_{ij}>0$, sample

Generalized Bradley-Terry models

Consider now that the pairwise comparisons are modeled using the Bradley-Terry model with “home-field advantage” (Agresti, , 1990) where

The parameter $\theta$, $\theta>0$, measures the strength of the home-field advantage ($\theta>1$) or disadvantage ($\theta<1$). Let $a_{ij}$ be the number of times that $i$ is at home and beats $j$ and $b_{ij}$ is the number of times that $i$ is at home and loses to $j$.

The log-likelihood of the skill ratings $\lambda$ and $\theta$ is given by

where $n_{ij}=a_{ij}+b_{ij}$ is the number of times $i$ plays at home against $j$, $c=\sum_{1\leq i\neq j\leq K}a_{ij}$ is the total number of home-field wins and $w_{i}$ is the total number of wins of element $i$.

For $1\leq i\neq j\leq K$ such that $n_{ij}>0$, let us introduce the latent variables $Z=\left\{Z_{ij}\right\}$ which are such that

The associated complete data log-likelihood is given by

Using independent priors for $\lambda$ and $\theta$, i.e. $p\left(\lambda,\theta\right)=p\left(\lambda\right)p\left(\theta\right)$, where $p\left(\lambda\right)$ is defined as (4) and

then we can maximize the resulting posterior using the EM algorithm

For $a=a_{\theta}=1$ and $b=b_{\theta}=0$, i.e. if we use flat priors, this EM algorithm is similar to the MM algorithm proposed in (Hunter, , 2004, pp. 389).

Using the same latent variables, we can sample from the posterior distribution of $\left(\lambda,\theta,Z\right)$ using the Gibbs sampler which updates iteratively $Z,$ $\lambda$ and $\theta$ as follows at iteration $t$:

For $1\leq i\neq j\leq K$ s.t. $n_{ij}>0$, sample

2 Model with ties

If we now want to allow for ties in pairwise comparisons, we can use the following model proposed by Rao and Kupper, (1967)

where $\theta>1$. The log-likelihood function for $\left(\lambda,\theta\right)$ is given by

where $t_{ij}=t_{ji}$ is the number of ties between $i$ and $j$ and $s_{ij}=w_{ij}+t_{ij}$.

For $1\leq i\neq j\leq K$ such that $s_{ij}>0$, let us introduce the latent variables $Z=\left\{Z_{ij}\right\}$ which are such that

which yields the following complete log-likelihood

where $T=\frac{1}{2}\sum_{1\leq i\neq j\leq K}t_{ij}$ is the total number of ties. If we adopt for $\theta$ a flat improper prior on $\left[1,\infty\right)$ and we select $p\left(\lambda\right)$ as (4) then we obtain

and we recover once again the minorizing function in (Hunter, , 2004, pp. 389-390) for $a=1$ and $b=0$. This can be maximized using the following procedure

Using the same latent variables, we can sample from the posterior distribution of $\left(\lambda,\theta,Z\right)$ using the following Gibbs sampler which updates iteratively $Z,$ $\lambda$ and $\theta$ as follows at iteration $t$:

For $1\leq i\neq j\leq K$ s.t. $s_{ij}>0$, sample

It is possible to sample from (10) exactly. By performing a change of variable $\overline{\theta}=\theta-1$, we obtain

which is a mixture of Gamma distributions.

3 Group comparisons

Consider now that we have $n$ pairwise comparisons betweem teams. For each comparison $i=1,\ldots,n$, let $T_{i}^{+}\subset\{1,\ldots,K\}$ be the winning team and $T_{i}^{-}\subset\{1,\ldots,K\}$ the losing team where $T_{i}^{+}\cap T_{i}^{-}=\varnothing$ and $T_{i}=T_{i}^{+}\cup T_{i}^{-}$. Recently Huang et al., (2006) have proposed the following model

The log-likelihood function for $\lambda$ is thus given by

For $i=1,...,n$ we introduce the latent variables $Z=\left\{Z_{i}\right\}$ and $C=\left\{C_{i}\right\}$ such that

where $\mathcal{E}\left(x;\alpha\right)$ is the exponential density of argument $x$ and parameter $\alpha$. It follows that the complete log-likelihood is given by

The $Q$ function associated to the EM algorithm is given by

where $\alpha_{ik}=1$ if $k\in T_{i}^{+}$ and otherwise and $\gamma_{ik}=1$ if $k\in T_{i}$ and otherwise. It follows that the EM update is given by

Using the same latent variables, we obtain a data augmentation sampler to sample from $p\left(\left.\lambda,z,c\right|D\right)$ by iteratively sampling $\left(Z,C\right)$ and $\lambda.$ This proceeds as follows at iteration $t$:

where $\delta_{u,v}=1$ if $u=v$ and otherwise.

Multiple comparisons

We now consider the popular Plackett-Luce model (Luce, , 1959; Plackett, , 1975) which is an extension of the Bradley-Terry model to comparisons involving more than two elements. Assume that $p_{i}\leq K$ individuals are ranked for comparison $i$ where $i=1,...,n$. We write $\rho_{i}=(\rho_{i1},\ldots,\rho_{ip_{i}})$ where $\rho_{i1}$ is the first individual, $\rho_{i2}$, the second, etc. The Plackett-Luce model assumes

For $i=1,\ldots,n$ and $j=1,\ldots,p_{i}-1$, we introduce the following latent variables $Z=\left\{Z_{ij}\right\}$

which leads to the complete log-likelihood

The $Q$ function associated to the EM algorithm is given by

which is once again equivalent to the majorizing function in (Hunter, , 2004, pp. 398) for $a=1,$ $b=0$. It follows that the EM algorithm is given at iteration $t$ by

where $w_{k}$ is the number of rankings where the $k^{\text{th}}$ individual is not in the last ranking position and $\delta_{ijk}$ is defined as

i.e. $\delta_{ijk}$ is the indicator of the event that individual $k$ receives a rank no better than $j$ in the $i^{\text{th}}$ ranking.

To sample from $p\left(\left.\lambda,z\right|D\right)$, we can use the following data augmentation sampler. At iteration $t$, this proceeds as follows:

For $i=1,...,n$, for $j=1,\ldots,p_{i}-1$, sample

Using exactly the same augmentation, EM and Gibbs samplers can be defined for further extensions of these models such as mixtures of Plackett-Luce models (Gormley and Murphy, , 2008).

Discussion

Consider the basic Bradley-Terry model and its extensions to group comparisons and multiple comparisons. Let us define

and write $\pi:=\left\{\pi_{i}\right\}_{i=1}^{K}$. The likelihood is invariant to a rescaling of the vector $\lambda$ so the parameter $\Lambda$ is not likelihood-identifiable and

From (4), it follows that $\pi\sim\mathcal{D}(a,\ldots,a)$ where $\mathcal{D}$ is the Dirichlet distribution and $\Lambda\sim\mathcal{G}(Ka,b)$, hence

To improve the mixing of the MCMC algorithms in this context, an additional sampling step can be added where we normalize the current parameter estimate $\lambda^{\left(t\right)}$ and then rescale them randomly using a prior draw for $\Lambda$.

For $i=1,\ldots,K$, set $\lambda_{i}^{\ast(t)}=\frac{\lambda_{i}^{(t)}}{\sum_{j=1}^{K}\lambda_{j}^{(t)}}\Lambda^{(t)}$ where $\Lambda^{(t)}\sim\mathcal{G}(Ka,b)$.

This step can drastically improve the mixing of the Markov chain. However, if we are only interested in the normalized values $\pi$ of $\lambda$ then this additional step is useless.

As an alternative, it is also possible to consider an EM algorithm for the basic Bradley-Terry model which does not require the introduction of a scale parameter. Assume $\pi\sim\mathcal{D}(a,\ldots,a)$ and let us introduce latent variables $M_{ij}$, $C_{ij}=(C_{ij1},\ldots,C_{ijM_{ij}})$ for $1\leq i\neq j\leq K$ such that $n_{ij}>0$

where $NB(r,p)$ is the negative binomial distribution. The complete log-likelihood is given by

where $r_{ijk}$ is the number of $c_{ijl},l=1,\ldots,m_{ij}$ that take value $k$. Omitting the terms independent of $\pi$, the $Q$ function is given by

where $C$ is a term independent of $\pi$. It follows that the EM update is given by

with $\sum_{k=1}^{K}\pi_{k}^{\left(t\right)}=1.$ Although the above EM algorithm does not rely on unidentifiable scaling parameters, it suffers from a slow convergence rate. When $\pi_{k}$ takes small values, $\sum_{i\neq k}\sum_{j\neq k}\frac{n_{ij}}{\pi_{i}+\pi_{j}}$ is large and it slows down the convergence of the EM algorithm. The same augmentation can be used to define a Gibbs sampler, but the same slow mixing issues arise for the Markov chain.

2 Hyperparameter estimation

The prior (4) is specified by the hyperparameters $a$ and $b$. However, the inverse scale parameter $b$ is not likelihood identifiable so there is no point assigning a prior to it. However it might be interesting to set a prior $p(a)$ on $a$ and estimate it from the data. Given $\lambda$, we have

It is possible to sample from this density using auxiliary variables $U_{1},U_{2}$ defined on $\left(0,\infty\right)$ as described in (Damien et al., , 1999). We introduce

A Gibbs sampler can now be implemented to sample from $p\left(\left.a,u_{1},u_{2}\right|\lambda\right)$. We can directly sample from the full conditionals of $U_{1}$ and $U_{2}$

where $\mathcal{U}\left(\alpha,\beta\right)$ is the uniform distribution on $\left(\alpha,\beta\right)$. The full conditional of $a$ given $u_{1},u_{2}$ is given by

Alternatively we can update $a$ using a M-H random walk on $\log(a)$. We can propose $a^{\star}=\exp(\sigma_{a}^{2}z)a$ where $z\sim\mathcal{N}(0,1)$ and $a^{\star}$ is accepted with probability

3 Further extensions

A model closely related to Bradley-Terry has been proposed for undirected random graphs with $K$ vertices (Holland and Leinhardt, , 1981; Chatterjee et al., , 2010; Park and Newman, , 2004). In this model, the degree sequence $(d_{1},\ldots,d_{K})$ of a given graph, where $d_{i}$ is the degree of node $i$, is supposed to capture all the information in the graph. It can be formalized by saying that the degree sequence is a sufficient statistic for a probability distribution on graphs (Chatterjee et al., , 2010).

In this model an edge is inserted between vertices $i$ and $j$ for $1\leq i<j\leq K$ with probability

where $\lambda_{k}>0$ for $k\in\{1,\ldots,K\}$. Let $r_{ij}=1$ if there is an edge between $i$ and $j$ and otherwise. Given the observations $D=\left\{r_{ij}\right\}_{1\leq i<j\leq K}$, the log-likelihood function for $\lambda$ is given by

We introduce the following latent variables $Z=\left\{Z_{ij}\right\}_{1\leq i<j\leq K}$ such that

The $Q$ function associated to the EM algorithm is given by

Solving $\partial Q(\lambda,\lambda^{\ast})/\partial\lambda_{i}=0$ requires solving a quadratic equation. For sake of brevity, we do not present these details here.

Once again, we can define a data augmentation sampler to sample from $p\left(\left.\lambda,z\right|D\right)$ by iteratively sampling $Z$ and $\lambda.$ This proceeds as follows at iteration $t$:

Here $\mathcal{GIG}\left(\alpha,\beta,\gamma\right)$ denotes the generalized inverse Gaussian distribution (see e.g. (Barndorff-Nielsen and Shephard, , 2001)) whose density for an argument $x$ is proportional to

Algorithms to sample exactly from this distribution are available.

3.2 Choice models

Other extensions of the Bradley-Terry model are the choice models introduced by Restle, (1961) and Tversky, 1972a ; Tversky, 1972b in psychology; see also (Wickelmaier and Shmidt, , 2004; Görür et al., , 2006). In these models, we are given a set of $n$ elements. To each element $i$ is associated a set of $K$ features represented by a binary vector $f_{i}\in\{0,1\}^{K}$. The probability that element $i$ is chosen over element $j$ is given by

where $\lambda_{k}>0$ is a weight representing the importance of feature $k$. The term $\sum_{k=1}^{K}\lambda_{k}f_{ik}(1-f_{jk})$ corresponds to the sum of the weights of features possessed by object $i$ but not object $j$. EM and Gibbs algorithms can be derived by following the same construction as with group comparisons.

3.3 Classification model

Let consider the following original model for categorical data analysis

Experimental results

In all the above models, the parameter $b$ is just a scaling parameter on $\lambda_{k}$. As the likelihood is invariant to a rescaling of the $\lambda_{k}$, this parameter does not have any influence on inference. Hence to ensure that the MAP estimate satisfies $\sum_{k=1}^{K}\widehat{\lambda}_{k}=1$, we set $b=Ka-1$ henceforth as explained in section 5. We demonstrate our algorithms on one synthetic and two real-world data sets.

We first study the Plackett-Luce model, comparing experimentally the mixing properties of the Gibbs sampler relative to a slightly modified version of the M-H algorithm proposed by Gormley and Murphy, (2009). In this latter paper, the authors propose to update the skill parameters simultaneously using the following proposal distributionThe authors actually use a normal approximation of the gamma distribution, and work with normalized data. To obtain similar algorithms, we consider unnormalized data. at iteration $t$

We simulated $500$ dataset of $n$ rankings of $K=4$ individuals, for various values of $n$ with $a=5$. For each dataset, 10,000 iterations of the Gibbs sampler presented in section 4 were run. The sample lag-1 autocorrelation was then computed for the four skill parameters. For a given sample size $n$, the mean value over skill parameters and simulated data is reported on Figure 1 together with 90% confidence bounds. The algorithm of Gormley and Murphy, (2009) performs reasonably well when the sample size is large, which is the case for the voting data they considered, but poorly for small sample sizes.

2 Nascar 2002 dataset

NASCAR is the primary sanctioning body for stock car auto racing in the United States. Each race involves 43 drivers. During the 2002 season, 87 different drivers participated in 36 races. Some drivers participated in all of the races while others participated in only one. We propose to apply the Plackett-Luce model with gamma prior on the parameters. The NASCAR datasetThe data can be downloaded from http://www.stat.psu.edu/ dhunter/code/btmatlab/ has been studied by Hunter, (2004), who noted that the MLE cannot be found for the original data set as four drivers placed last in each race they entered, and therefore had to be removed. This does not need to be done if we follow a Bayesian approach. We focus here on predicting the outcome of the next race based on the previous ones, starting from race 5; i.e. we predict the results of race 6 based on the MAP estimates obtained with the first 5 races, then the results of race 7 based on the MAP estimated obtained with the first 6 races, etc. For each race, we compute the test log-likelihood using the MAP estimates. The mean value and 90% confidence bounds are represented in Figure 2 w.r.t. the value of $a$. The EM algorithm was initialized using $(\lambda_{1}^{\left(0\right)},\ldots,\lambda_{83}^{\left(0\right)})=(\frac{1}{83},\ldots,\frac{1}{83})$.

The Gibbs sampler was also applied to the same dataset. The skill parameters were initialized at the same value, and the parameter $a$ was assigned a flat improper prior and sampled as described in section 5. We ran 50,000 iterations with 2,000 burn-in. As detailed in Section 5, only the normalized weights $\pi_{i}$ are likelihood identifiable. Skill ratings are usually represented on the real line, and we use the following one-to-one mapping $\beta_{i}=\log\pi_{i}-\log 1/83$. The marginal posterior densities of the reparameterized skill ratings for the first four drivers according to their average place are reported in Figure 3. The Bayesian approach can effectively capture the uncertainty in the skill ratings of the drivers. ML and MMSE (minimum mean squared error) estimates together with standard deviations are reported in Table 1 for the first ten and last ten drivers according to average place.

3 Chess data

Rating the skills of chess players is of major practical interest. It allows organizers of a tournament to avoid having strong players playing against each other at early stages, or to restrict the tournament to players with skills above a given threshold. The international chess federation adopted the so-called “Elo” system which is based on the Bradley-Terry model (Elo, , 1978). For historical considerations about the rating system in chess, the reader should refer to Glickman, (1995).

We consider here game-by-game chess results over 100 months, consisting of 65,053 matches between 8631 playersChess data can be downloaded from http://kaggle.com/chess. The outcome of each game is either win (+1), tie (+0.5) or loss (0). We estimate the parameters of the Bradley-Terry model with ties presented in section 3.2 on the first 95 months and then predict the outcome of the games of the last 5 months. The hyperparameters for the tie parameter $\theta$ are set to $a_{\theta}=1,$ $b_{\theta}=0$. Given the large sample size, it is not possible to sample from Eq. (10) as the number of elements in the mixture is very large. We therefore use a M-H step with a normal random walk proposal of standard deviation $0.1$. The mean squared error is reported for predictions based on MAP estimates and full Bayesian predictive based on the Gibbs sampler outcomes, for different values of the hyperparameter $a$. EM and Gibbs samplers were initialized at $(\lambda_{1}^{\left(0\right)},\ldots,\lambda_{8631}^{\left(0\right)})=(\frac{1}{8631},\ldots,\frac{1}{8631})$ and $\theta^{\left(0\right)}=1,5$. The Gibbs samplers were run with 10,000 iterations and 1,000 burn-in iterations. The results are reported in Figure 4. The results demonstrate the benefit of penalizing the skill rating parameters and the improvement brought up by a full Bayesian analysis. In Figure 5 we also report the autocorrelation function associated to the parameter $\theta$ and the skill parameters with largest mean values. The Markov chain displays good mixing properties.

Conclusion

The Bradley-Terry model and its generalizations arise in numerous applications. We have shown here that most of the MM algorithms proposed in Hunter, (2004) can be reinterpreted as special cases of EM algorithms. Additionally we have proposed original EM algorithms for some recent generalizations of the Bradley-Terry models. Finally we have shown how the latent variables introduced to derive these EM algorithms lead straightforwardly to Gibbs sampling algorithms. These elegant MCMC algorithms mix experimentally well and outperform a recently proposed M-H algorithm.

Acknowledgment. The authors are grateful to Persi Diaconis for helpful discussions and pointers to references on the Plackett-Luce and random graph models and to Luke Bornn for helpful comments.

References