Как вписать модель Брэдли – Терри – Люса в R без сложной формулы?

9

Модель Брэдли – Терри – Люса (BTL) утверждает, что , где - вероятность того, что объект j будет оценен как «лучший», тяжелее, и т. д., чем объект i , и \ delta_i , и \ delta_j являются параметрами.pji=logit1(δjδi)pijjiδiδj

Кажется, это кандидат на функцию glm с family = binomial. Однако формула будет выглядеть примерно так: «Success ~ S1 + S2 + S3 + S4 + ...», где Sn - фиктивная переменная, то есть 1, если объект n является первым объектом в сравнении, -1, если это второе и 0 в противном случае. Тогда коэффициент Sn будет соответствующим deltan .

Это было бы довольно легко управлять только несколькими объектами, но могло бы привести к очень длинной формуле и необходимости создания фиктивной переменной для каждого объекта. Мне просто интересно, есть ли более простой метод. Предположим, что имя или номер двух сравниваемых объектов являются переменными (факторами?) Object1 и Object2, и Success равен 1, если объект 1 оценен лучше, и 0, если объект 2 является.

тарпон
источник
3
Существует пакет R для модели Брэдли-Терри. Посмотри на Rseek.
кардинал
Я также предоставил несколько ссылок на связанный вопрос: stats.stackexchange.com/a/10741/930
chl
Упомянутый пакет @cardinal, кстати: BradleyTerry2
сопряженный

Ответы:

17

Я думаю, что лучшим пакетом для данных Paired Comparison (PC) в R является пакет prefmod , который позволяет удобно подготовить данные для подгонки (логарифмировать) BTL-моделей в R. Он использует Poisson GLM (точнее, многочленный логит в Poisson). Формулировка см., например, это обсуждение ).

Приятно то, что у него есть функция, prefmod::llbt.designкоторая автоматически конвертирует ваши данные в нужный формат и необходимую матрицу дизайна.

Например, скажем, у вас есть 6 объектов для сравнения. затем

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

создаст матрицу дизайна из матрицы данных, которая выглядит следующим образом:

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

со строками, обозначающими людей, столбцы, обозначающие сравнения, и 0 означает неопределенность, 1 означает предпочтительный объект 1, а 2 означает предпочтительный объект 2. Пропущенные значения допускаются. Изменить : Поскольку это, вероятно, не то, что выводить просто из данных выше, я изложил это здесь. Сравнения должны быть упорядочены следующим образом ((12) означает сравнение объекта 1 с объектом 2):

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

Подгонка удобнее всего выполнять с помощью gnm::gnmфункции, поскольку она позволяет выполнять статистическое моделирование. (Редактировать: Вы также можете использовать prefmod::llbt.fitфункцию, которая немного проще, так как она требует только количества и матрицы дизайна.)

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

Обратите внимание, что термин «исключение» пропустит параметры помех в сводке. Затем вы можете получить параметры стоимости (ваши дельты) как

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

И вы можете построить их с

R> plotworth(wmat)

Если у вас много объектов и вы хотите o1+o2+...+onбыстро написать объект формулы , вы можете использовать

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

создать формулу gnm( для которой вам не нужно llbt.fit).

Существует статья JSS , см. Также https://r-forge.r-project.org/projects/prefmod/ и документацию через ?llbt.design.

Момо
источник
1
Это очень тщательный ответ. Спасибо. Кажется, что prefmod будет хорошим пакетом для использования. Кстати, мне интересно использовать модель, чтобы попытаться предсказать результаты спортивных матчей.
Серебряная рыба
Нет проблем, рад, если это помогло. Я точно не знаю, как вы хотите предсказать, но Leitner et al. использовали эти модели для прогнозирования спортивных событий. Смотрите его диссертацию epubdev.wu.ac.at/2925 . Удачи.
Момо
Может быть, эта ссылка лучше epubdev.wu.ac.at/view/creators/…
Momo
Можно ли рассчитать значения для различий между отдельными парами (например, o1 и o2) из ​​этих данных? Или вам нужно изменить формулу, использовать o2 как последний фактор и жить без оценки Std.error в этом случае?
TNT
1
Это было какое - то время, так что я не помню , можно ли удобно использовать линейные ограничения , но то , что вы можете сделать в вашем случае использовать один в качестве опорного уровня, скажем , o1, и использовать т значение другого, скажем , o2, Исходя из резюме, он фактически представляет собой тест, равный нулю разницы между o1 и o2.
Момо