Como encaixar o modelo de Bradley-Terry-Luce em R, sem fórmula complicada?

9

O modelo de Bradley-Terry-Luce (BTL) afirma que , em que é a probabilidade de o objeto ser considerado "melhor", mais pesado, etc, do que o objecto i , e \ delta_i , e \ delta_j são parâmetros.p i j j ipjEu=euogEut-1 1(δj-δEu)pEujjEuδEuδj

Parece ser um candidato para a função glm, com família = binomial. No entanto, a fórmula seria algo como "Sucesso ~ S1 + S2 + S3 + S4 + ...", em que Sn é uma variável dummy, ou seja, 1 se o objeto n for o primeiro objeto na comparação, -1 se for o segundo e 0 caso contrário. Então o coeficiente de Sn seria o delta_n correspondente deeutuman.

Isso seria bastante fácil de gerenciar com apenas alguns objetos, mas poderia levar a uma fórmula muito longa e à necessidade de criar uma variável dummy para cada objeto. Eu só me pergunto se existe um método mais simples. Suponha que o nome ou o número dos dois objetos que estão sendo comparados sejam variáveis ​​(fatores?) Objeto1 e Objeto2, e Sucesso será 1 se o objeto 1 for julgado melhor e 0 se o objeto 2 for.

Silverfish
fonte
3
Existe um pacote R para o modelo Bradley-Terry. Olhe para Rseek.
cardeal
Também forneci alguns links para uma pergunta relacionada: stats.stackexchange.com/a/10741/930
chl
O pacote @ cardinal mencionado, btw: BradleyTerry2
conjugateprior

Respostas:

17

Eu acho que o melhor pacote para dados de Comparação emparelhada (PC) em R é o pacote prefmod , que permite preparar convenientemente dados para ajustar modelos BTL (log linear) em R. Ele usa um Poisson GLM (mais precisamente, um logit multinomial em Poisson formulação ver, por exemplo, esta discussão ).

O bom é que ele possui uma função prefmod::llbt.designque converte automaticamente seus dados no formato e na matriz de design necessários.

Por exemplo, digamos que você tenha 6 objetos, todos emparelhados. Então

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

criará a matriz de design a partir de uma matriz de dados que se parece com isso:

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

com linhas indicando pessoas, colunas indicando comparações e 0 significa indeciso 1 significa objeto 1 preferido e 2 significa objeto 2 preferido. Valores ausentes são permitidos. Edit : Como isso provavelmente não é algo para inferir simplesmente a partir dos dados acima, eu soletrar aqui. As comparações devem ser ordenadas da seguinte maneira ((12): objeto de comparação médio 1 com o objeto 2):

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

O ajuste é mais convenientemente realizado com a gnm::gnmfunção, pois permite a modelagem estatística. (Editar: você também pode usar a prefmod::llbt.fitfunção, que é um pouco mais simples, pois requer apenas as contagens e a matriz de design.)

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

Observe que o termo de eliminação omitirá os parâmetros incômodos do resumo. Você pode obter os parâmetros de valor (seus deltas) como

## 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

E você pode plotá-los com

R> plotworth(wmat)

Se você possui muitos objetos e deseja escrever o1+o2+...+onrapidamente um objeto de fórmula , pode usar

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

para gerar a fórmula gnm(da qual você não precisaria llbt.fit).

Há um artigo JSS , consulte também https://r-forge.r-project.org/projects/prefmod/ e a documentação via ?llbt.design.

Momo
fonte
11
Essa é uma resposta muito completa. Obrigado. Parece que prefmod seria um bom pacote para usar. Estou interessado em usar o modelo para tentar prever os resultados de jogos esportivos, a propósito.
quer
Não tem problema, feliz se ajudou. Não sei exatamente como você pretende prever, mas Leitner et al. usaram esses modelos para prever eventos esportivos. Veja sua tese epubdev.wu.ac.at/2925 . Boa sorte.
Momo
Talvez este link seja melhor epubdev.wu.ac.at/view/creators/…
Momo
É possível calcular significados para as diferenças entre pares individuais (por exemplo, o1 e o2) a partir desses dados? Ou você precisa reorganizar a fórmula, usar o2 como último fator e viver sem uma estimativa de erro padrão nesse caso?
TNT
11
Já faz um tempo, então não me lembro se você pode usar restrições lineares convenientemente, mas o que você pode fazer no seu caso é usar um como nível de referência, digamos o1, e usar o valor t do outro, digamos o2, do resumo - efetivamente constitui um teste para determinar se a diferença entre o1 e o2 é zero.
Momo