# How to fit Bradley–Terry–Luce model in R, without complicated formula?

The Bradley–Terry–Luce(BTL) model states that $p_{ji} = logit^{-1}(\delta_j - \delta_i)$, where $p_{ij}$ is the probability that object $j$ is judged to be "better", heavier, etc, than object $i$, and $\delta_i$, and $\delta_j$ are parameters.

This seems to be a candidate for the glm function, with family = binomial. However, the formula would be something like "Success ~ S1 + S2 + S3 + S4 +...", where Sn is a dummy variable, that is 1 if object n is the first object in the comparison, -1 if it is the second, and 0 otherwise. Then the coefficient of Sn would be the corresponding $delta_n$.

This would be fairly easy to manage with only a few objects, but could lead to a very long formula, and the need to create a dummy variable for each object. I just wonder if there is a simpler method. Suppose that the name or number of the two objects being compared are variables (factors?) Object1, and Object2, and Success is 1 if object 1 is judged better, and 0 if object 2 is.

• There is an R package for the Bradley-Terry model. Look on Rseek. Apr 24, 2012 at 2:05
• I also provided some links on a related question: stats.stackexchange.com/a/10741/930
– chl
Apr 24, 2012 at 6:58
• The package @cardinal mentioned, btw: BradleyTerry2 Jul 20, 2012 at 13:55

I think the best package for Paired Comparison (PC) data in R is the prefmod package, which allows to conveniently prepare data to fit (log linear) BTL models in R. It uses a Poisson GLM (more accurately, a multinomial logit in Poisson formulation see e.g. this discussion).

The nice thing is that it has a function prefmod::llbt.design that automatically converts your data into the necessary format and necessary design matrix.

For example, say you have 6 objects all pairwise compared. Then

R> library(prefmod)
R> des<-llbt.design(data, nitems=6)

will build the design matrix from a data matrix that looks like this:

P1  0  0 NA  2  2  2  0  0  1   0   0   0   1   0   1   1   2
P2  0  0 NA  0  2  2  0  2  2   2   0   2   2   0   2   1   1
P3  1  0 NA  0  0  2  0  0  1   0   0   0   1   0   1   1   2
P4  0  0 NA  0  2  0  0  0  0   0   0   0   0   0   2   1   1
P5  0  0 NA  2  2  2  2  2  2   0   0   0   0   0   2   2   2
P6  2  2 NA  0  0  0  2  2  2   2   0   0   0   0   2   1   2

with rows denoting people, columns denoting comparisons and 0 means undecided 1 means object 1 preferred and 2 means object 2 preferred. Missing values are allowed. Edit: As this is probably not something to infer simply from the data above, I spell it out here. The comparisons must be ordered the following way((12) mean comparison object 1 with object 2):

(12) (13) (23) (14) (24) (34) (15) (25) etc.

Fitting is than most conveniently carried out with the gnm::gnm function, as it allows you to do statistical modelling. (Edit: You can also use the prefmod::llbt.fit function, which is a bit simpler as it takes only the counts and the design matrix.)

R> res<-gnm(y~o1+o2+o3+o4+o5+o6, eliminate=mu, family=poisson, data=des)
R> summary(res)
Call:
gnm(formula = y ~ o1 + o2 + o3 + o4 + o5 + o6, eliminate = mu,
family = poisson, data = des)

Deviance Residuals:
Min      1Q  Median      3Q     Max
-7.669  -4.484  -2.234   4.625  10.353

Coefficients of interest:
Estimate Std. Error z value Pr(>|z|)
o1  1.05368    0.04665  22.586  < 2e-16 ***
o2  0.52833    0.04360  12.118  < 2e-16 ***
o3  0.13888    0.04297   3.232  0.00123 **
o4  0.24185    0.04238   5.707 1.15e-08 ***
o5  0.10699    0.04245   2.521  0.01171 *
o6  0.00000         NA      NA       NA
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Std. Error is NA where coefficient has been constrained or is unidentified

Residual deviance: 2212.7 on 70 degrees of freedom
AIC: 2735.3

Please note that the eliminate term will omit the nuisance parameters from the summary. You then can get the worth parameters (your deltas) as

## calculating and plotting worth parameters
R> wmat<-llbt.worth(res)
worth
o1 0.50518407
o2 0.17666128
o3 0.08107183
o4 0.09961109
o5 0.07606193
o6 0.06140979

And you can plot them with

R> plotworth(wmat)

If you have many objects and want to write a formula object o1+o2+...+on fast, you can use

R> n<-30
R> objnam<-paste("o",1:n,sep="")
R> fmla<-as.formula(paste("y~",paste(objnam, collapse= "+")))
R> fmla
y ~ o1 + o2 + o3 + o4 + o5 + o6 + o7 + o8 + o9 + o10 + o11 +
o12 + o13 + o14 + o15 + o16 + o17 + o18 + o19 + o20 + o21 +
o22 + o23 + o24 + o25 + o26 + o27 + o28 + o29 + o30

to generate the formula for gnm (which you wouldn't need for llbt.fit).

There is a JSS article, see also https://r-forge.r-project.org/projects/prefmod/ and the documentation via ?llbt.design.

• That's a very thorough response. Thank you. It seems like prefmod would be a good package to use. I'm interested in using the model to try to predict the results of sports matches, by the way. Apr 25, 2012 at 17:08
• No problem, glad if it helped. I don't exactly know how you mean to predict, but Leitner et al. have used these models to predict sport events. See his thesis epubdev.wu.ac.at/2925. Good luck.
– Momo
Apr 25, 2012 at 19:18
• Maybe this link is better epubdev.wu.ac.at/view/creators/…
– Momo
Apr 25, 2012 at 19:21
• Is it possible to calculate significances for the differences between individual pairs (e.g. o1 and o2) from this data? Or do you have to rearrange the formula, use o2 as last factor and live without a Std.error estimation in that case?
– TNT
Dec 12, 2016 at 7:01
• It's been a while, so I don't recall whether you can conveniently use linear restrictions but what you can do in your case is to use one as the reference level, say o1, and use the t value of the other, say o2, from the summary - it effectively constitutes a test whether the difference between o1 and o2 is zero.
– Momo
Dec 12, 2016 at 11:50