Introducción

Este cuaderno busca abordar el modelamiento de preferencias a través de comparaciones pareadas, utilizando como referencia el modelo de Bradley-Terry. El conjunto de datos utilizado proviene del paquete MPsychoR, y presenta los resultados de una encuesta en la que 200 personas compararon pares de bandas musicales.

1. ¿Qué es el modelo de Bradley-Terry?

El modelo de Bradley-Terry asigna una “habilidad” a cada ítem (por ejemplo, una banda) en una escala numérica. Esta habilidad nos permite predecir la probabilidad de que un ítem sea preferido sobre otro cuando se comparan de a pares.

2. ¿Cómo se calcula paso a paso?

A continuación desglosamos cada fase, desde los datos originales hasta la obtención de los parámetros.

2.1. Datos observados (formato binomial)

library(MPsychoR)
data("bandpref")
kable(head(bandpref), caption = "Primeras filas de bandpref")
Primeras filas de bandpref
Band1 Band2 Win1 Win2
Slayer Rush 142 58
Slayer Death 54 146
Slayer Emperor 158 42
Slayer Scorpions 34 166
Rush Death 121 79
Rush Emperor 147 53
  • Para cada par de bandas (i, j), contamos:
    • Win1: número de veces que i venció a j.
    • Win2: número de veces que j venció a i.

2.2. Transformación interna a matriz de diseño y vector de respuesta

# Construcción interna de X e Y usando BTm.setup
setup <- BradleyTerry2:::BTm.setup(
  outcome = cbind(bandpref$Win1, bandpref$Win2),
  player1 = bandpref$Band1,
  player2 = bandpref$Band2,
  data    = bandpref
)

# Mostrar las primeras 6 filas de la matriz de diseño X
knitr::kable(setup$X, caption = "Matriz de diseño X")
Matriz de diseño X
..Rush ..Death ..Emperor ..Scorpions
-1 0 0 0
0 -1 0 0
0 0 -1 0
0 0 0 -1
1 -1 0 0
1 0 -1 0
1 0 0 -1
0 1 -1 0
0 1 0 -1
0 0 1 -1 
# Mostrar las primeras 6 entradas del vector de respuesta Y (Wins vs Losses)
knitr::kable(data.frame(Win1 = bandpref$Win1, Win2 = bandpref$Win2), 
             caption = "Vector de respuesta Y: Win1 vs Win2")
Vector de respuesta Y: Win1 vs Win2
Win1 Win2
142 58
54 146
158 42
34 166
121 79
147 53
155 45
72 128
140 60
159 41
2.2.1 ¿Cómo se explica este diseño en el modelo?

Cada fila de setup$X indica qué bandas se enfrentan: +1 para la banda en player1 (ganadora en el diseño) y -1 para la banda en player2 (perdedora), con 0 en las demás columnas.

El vector de respuesta Y se construye a partir de dos columnas:

  • Win1: número de veces que player1 ganó.
  • Win2: número de veces que player2 ganó.

Cuando usamos en glm() la fórmula:

cbind(Win1, Win2) ~ X - 1

le decimos a R que ajuste un modelo binomial para los conteos de éxitos y fracasos en cada comparación:

  • La parte izquierda (cbind(Win1, Win2)) entrega un par de números que representa el número de éxitos (Win1) y fracasos (Win2).
  • La parte derecha (X - 1) utiliza la matriz de diseño sin intercepto, de modo que cada coeficiente estimado corresponde directamente a la habilidad logarítmica λᵢ de cada banda.

En otras palabras, R modela la proporción de éxito en cada fila (Win1 / (Win1+Win2)) como función lineal en logit de los valores +1/−1 de X. Esto produce coeficientes que son exactamente los λᵢ del modelo Bradley–Terry.

# Preparar data frame para glm
# ... (resto del chunk permanece igual)
# Preparar data frame para glm
df_manual <- cbind(
  data.frame(
    Win1 = bandpref$Win1,
    Win2 = bandpref$Win2
  ),
  as.data.frame(setup$X)
)

# Ajuste manual del modelo: cbind(Win1, Win2) ~ X - 1
fit_manual <- glm(
  cbind(Win1, Win2) ~ . -1,
  data   = df_manual,
  family = binomial(link = "logit")
)

# Coeficientes estimados: equivalen a las habilidades logarítmicas λᵢ
knitr::kable(
  data.frame(
    Banda  = names(coef(fit_manual)),
    Lambda = round(coef(fit_manual), 3)
  ),
  caption = "Estimación manual de λᵢ con glm(binomial)"
)
Estimación manual de λᵢ con glm(binomial)
Banda Lambda
..Rush ..Rush 0.380
..Death ..Death 0.199
..Emperor ..Emperor -0.024
..Scorpions ..Scorpions -0.312

2.3. Ajuste del modelo con BTm. Ajuste del modelo con BTm Ajuste del modelo con BTm

# En lugar de armar manualmente la regresión, usamos directamente BTm:
bandsBT <- BTm(
  cbind(Win1, Win2),
  Band1,
  Band2,
  data = bandpref
)
summary(bandsBT)
## 
## Call:
## BTm(outcome = cbind(Win1, Win2), player1 = Band1, player2 = Band2, 
##     data = bandpref)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## ..Rush       0.37988    0.09099   4.175 2.98e-05 ***
## ..Death      0.19891    0.09027   2.204 0.027556 *  
## ..Emperor   -0.02434    0.09008  -0.270 0.786976    
## ..Scorpions -0.31179    0.09100  -3.426 0.000612 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 488.94  on 10  degrees of freedom
## Residual deviance: 423.67  on  6  degrees of freedom
## AIC: 486.45
## 
## Number of Fisher Scoring iterations: 4

2.4. Extracción de parámetros (λᵢ) con BTabilities. Extracción de parámetros (λᵢ) con BTabilities

bandsBT  <- BTm(cbind(Win1, Win2), Band1, Band2, data = bandpref)
bandsAbil <- BTabilities(bandsBT)
lambda    <- round(bandsAbil[, "ability"],  3)
kable(data.frame(Banda = names(lambda), Lambda = lambda),
      caption = "Habilidades logarítmicas (λ)")
Habilidades logarítmicas (λ)
Banda Lambda
Slayer Slayer 0.000
Rush Rush 0.380
Death Death 0.199
Emperor Emperor -0.024
Scorpions Scorpions -0.312

2.5. Conversión a escala de odds y probabilidad directa

alpha       <- exp(lambda)  # odds relativos
prob_matrix <- outer(alpha, alpha, function(a,b) round(a/(a+b),2))

Para un par específico:

alpha["Rush"] / (alpha["Rush"] + alpha["Death"])  # P(Rush > Death)
##      Rush 
## 0.5451269

Interpretación de betas (λᵢ), odds y probabilidades

  • Cada coeficiente λᵢ (beta) es el log‑odds de la habilidad de la banda i.
  • Al aplicar la función exponencial obtenemos αᵢ = e^{λᵢ}, que son las odds relativas:
    • Por ejemplo, si αᵢ = 2, significa que la banda i tiene el doble de chances de ganar frente a una banda con α = 1.
  • Para convertir las odds de un solo ítem en una probabilidad unaria (frente a un oponente con odds = 1), usamos:
$$ P_i = \frac{\alpha_i}{1 + \alpha_i} $$
  • Ejemplo: si \(α_i = 2\), entonces \(P_i = 2/(1+2) = 0.667\), lo que equivale al 66.7% de probabilidad de ganar.
  • Para comparar dos bandas i y j, la probabilidad de que i gane se calcula como:
$$ P(i > j) = \frac{\alpha_i}{\alpha_i + \alpha_j} $$

Esta fórmula extiende la interpretación de las odds al caso de enfrentamientos entre pares.

Interpretación de P(Rush > Death):

La probabilidad estimada de que Rush sea preferido por encima de Death es de 0.545 (54.5 %).
En otras palabras, de cada 100 enfrentamientos teóricos entre ambas bandas, Rush ganaría aproximadamente 55 veces y Death 45 veces.

En términos de odds, la razón de probabilidades es:

\[ \frac{\alpha_{\text{Rush}}}{\alpha_{\text{Death}}} \;=\; \frac{e^{\lambda_{\text{Rush}}}}{e^{\lambda_{\text{Death}}}} \;\approx\; \frac{1.463}{1.220} \;\approx\; 1.20 \]

Esto indica que Rush tiene un 20 % más de oportunidades de ser elegido frente a Death e

3. ¿Por qué escala logarítmica?.

  • Convierte comparaciones multiplicativas (odds) en diferencias aditivas.
  • Garantiza probabilidades válidas (0–1).

Resultados Completos

Tabla de habilidades y odds

result <- data.frame(
  Banda  = rownames(bandsAbil),
  Lambda = round(bandsAbil[,"ability"],3),
  Odds   = round(alpha, 3)
)
knitr::kable(result, caption = "Habilidades (λ) y Odds relativos")
Habilidades (λ) y Odds relativos
Banda Lambda Odds
Slayer Slayer 0.000 1.000
Rush Rush 0.380 1.462
Death Death 0.199 1.220
Emperor Emperor -0.024 0.976
Scorpions Scorpions -0.312 0.732

Heatmap de probabilidades

library(reshape2); library(ggplot2)
prob_df <- melt(prob_matrix)
colnames(prob_df) <- c("Ganador","Perdedor","Probabilidad")

ggplot(prob_df, aes(Perdedor, Ganador, fill=Probabilidad)) +
  geom_tile(color="white") +
  geom_text(aes(label=round(Probabilidad,2))) +
  scale_fill_gradient(low="white", high="steelblue") +
  labs(title="Probabilidad de victoria entre pares") +
  theme_minimal()

Top 5 comparaciones más y menos probables

prob_rank   <- subset(prob_df, Ganador!=Perdedor)
prob_top    <- head(prob_rank[order(-prob_rank$Probabilidad),],5)
prob_bottom <- head(prob_rank[order(prob_rank$Probabilidad),],5)

knitr::kable(prob_top,    caption="Top 5 más probables")
Top 5 más probables
Ganador Perdedor Probabilidad
22 Rush Scorpions 0.6663372
23 Death Scorpions 0.6249705
17 Rush Emperor 0.5997013
2 Rush Slayer 0.5938432
21 Slayer Scorpions 0.5773219 
knitr::kable(prob_bottom, caption="Top 5 menos probables")
Top 5 menos probables
Ganador Perdedor Probabilidad
10 Scorpions Rush 0.3336628
15 Scorpions Death 0.3750295
9 Emperor Rush 0.4002987
6 Slayer Rush 0.4061568
5 Scorpions Slayer 0.4226781

Gráficos de barras

prob_top$Comparación    <- paste(prob_top$Ganador, "vs", prob_top$Perdedor)
prob_bottom$Comparación <- paste(prob_bottom$Ganador, "vs", prob_bottom$Perdedor)

ggplot(prob_top, aes(reorder(Comparación, -Probabilidad), Probabilidad)) +
  geom_col(fill="darkgreen") + coord_flip() +
  labs(title="Top 5 más probables", x="", y="Probabilidad") +
  theme_minimal()

ggplot(prob_bottom, aes(reorder(Comparación, Probabilidad), Probabilidad)) +
  geom_col(fill="firebrick") + coord_flip() +
  labs(title="Top 5 menos probables", x="", y="Probabilidad") +
  theme_minimal()