¿Cómo adaptar el modelo Bradley – Terry – Luce en R, sin una fórmula complicada?

El modelo Bradley – Terry – Luce (BTL) establece que , donde es la probabilidad de que el objeto se considere "mejor", más pesado, etc., que el objeto , y , y son parámetros.$p_{ji} = logit^{-1}(\delta_j - \delta_i)$$p_{ij}$$j$$i$$\delta_i$$\delta_j$

Esto parece ser un candidato para la función glm, con family = binomial. Sin embargo, la fórmula sería algo así como "Éxito ~ S1 + S2 + S3 + S4 + ...", donde Sn es una variable ficticia, es decir 1 si el objeto n es el primer objeto en la comparación, -1 si es el segundo, y 0 en caso contrario. Entonces el coeficiente de Sn sería el correspondiente .$delta_n$

Esto sería bastante fácil de administrar con solo unos pocos objetos, pero podría conducir a una fórmula muy larga y la necesidad de crear una variable ficticia para cada objeto. Solo me pregunto si hay un método más simple. Suponga que el nombre o el número de los dos objetos que se comparan son variables (¿factores?) Object1 y Object2, y Success es 1 si el objeto 1 se juzga mejor, y 0 si el objeto 2 lo es.

Creo que el mejor paquete para datos de Comparación Pareada (PC) en R es el paquete prefmod , que permite preparar convenientemente los datos para que se ajusten (log lineales) a modelos BTL en R. Utiliza un Poisson GLM (más exactamente, un logit multinomial en Poisson formulación ver, por ejemplo, esta discusión ).

Lo bueno es que tiene una función prefmod::llbt.designque convierte automáticamente sus datos en el formato necesario y la matriz de diseño necesaria.

Por ejemplo, supongamos que tiene 6 objetos todos comparados por pares. Entonces

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

construirá la matriz de diseño a partir de una matriz de datos que se ve así:

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

con filas que indican personas, columnas que indican comparaciones y 0 significa indeciso 1 significa objeto 1 preferido y 2 significa objeto 2 preferido. Se permiten valores perdidos. Editar : como esto probablemente no sea algo que se infiera simplemente de los datos anteriores, lo deletreo aquí. Las comparaciones deben ordenarse de la siguiente manera ((12) objeto de comparación de medias 1 con objeto 2):

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

El ajuste se realiza más convenientemente con la gnm::gnmfunción, ya que le permite realizar modelos estadísticos. (Editar: también puede usar la prefmod::llbt.fitfunción, que es un poco más simple ya que solo toma los recuentos y la matriz de diseño).

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

Tenga en cuenta que el término eliminar omitirá los parámetros molestos del resumen. Luego puede obtener los parámetros de valor (sus 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

Y puedes trazarlos con

R> plotworth(wmat)

Si tiene muchos objetos y desea escribir un objeto de fórmula o1+o2+...+onrápidamente, puede 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 generar la fórmula para gnm(que no necesitarías llbt.fit).

Hay un artículo de JSS , consulte también https://r-forge.r-project.org/projects/prefmod/ y la documentación a través de ?llbt.design.