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

Question

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.

Answer 1

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.