Bagaimana agar sesuai dengan model Bradley-Terry-Luce dalam R, tanpa formula yang rumit?

9

Model Bradley – Terry – Luce (BTL) menyatakan bahwa , di mana adalah probabilitas objek yang dinilai "lebih baik", lebih berat, dll, daripada objek , dan , dan adalah parameter.p i j j i δ i δ jpji=logit1(δjδi)pijjiδiδj

Ini tampaknya menjadi kandidat untuk fungsi glm, dengan family = binomial. Namun, rumusnya akan seperti "Sukses ~ S1 + S2 + S3 + S4 + ...", di mana Sn adalah variabel dummy, yaitu 1 jika objek n adalah objek pertama dalam perbandingan, -1 jika itu adalah yang kedua, dan 0 sebaliknya. Maka koefisien Sn akan menjadi sesuai .deltan

Ini akan cukup mudah untuk dikelola dengan hanya beberapa objek, tetapi dapat menyebabkan formula yang sangat panjang, dan kebutuhan untuk membuat variabel dummy untuk setiap objek. Saya hanya ingin tahu apakah ada metode yang lebih sederhana. Misalkan nama atau jumlah dari dua objek yang dibandingkan adalah variabel (faktor?) Object1, dan Object2, dan Sukses adalah 1 jika objek 1 dinilai lebih baik, dan 0 jika objek 2 adalah.

Gegat
sumber
3
Ada paket R untuk model Bradley-Terry. Lihatlah Rseek.
kardinal
Saya juga memberikan beberapa tautan pada pertanyaan terkait: stats.stackexchange.com/a/10741/930
chl
Paket @ cardinal disebutkan, btw: BradleyTerry2
conjugateprior

Jawaban:

17

Saya pikir paket terbaik untuk data Paired Comparison (PC) dalam R adalah paket prefmod , yang memungkinkan untuk dengan mudah menyiapkan data agar sesuai (log linear) model BTL di R. Ia menggunakan Poisson GLM (lebih tepatnya, multinomial logit di Poisson formulasi lihat misalnya diskusi ini ).

Yang menyenangkan adalah ia memiliki fungsi prefmod::llbt.designyang secara otomatis mengubah data Anda ke dalam format yang diperlukan dan matriks desain yang diperlukan.

Misalnya, Anda memiliki 6 objek yang semuanya berpasangan dibandingkan. Kemudian

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

akan membangun matriks desain dari matriks data yang terlihat seperti ini:

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

dengan baris yang menunjukkan orang, kolom yang menunjukkan perbandingan dan 0 berarti ragu-ragu 1 berarti objek 1 disukai dan 2 berarti objek 2 disukai. Nilai yang hilang diizinkan. Sunting : Karena ini mungkin bukan sesuatu yang dapat disimpulkan hanya dari data di atas, saya mengejanya di sini. Perbandingan harus dipesan dengan cara berikut ((12) berarti objek perbandingan 1 dengan objek 2):

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

Pemasangan lebih mudah dilakukan dengan gnm::gnmfungsi, karena memungkinkan Anda untuk melakukan pemodelan statistik. (Sunting: Anda juga dapat menggunakan prefmod::llbt.fitfungsi, yang sedikit lebih sederhana karena hanya menggunakan hitungan dan matriks desain.)

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

Harap perhatikan bahwa istilah eliminasi akan menghilangkan parameter gangguan dari ringkasan. Anda kemudian bisa mendapatkan parameter layak (delta Anda) sebagai

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

Dan Anda dapat merencanakannya

R> plotworth(wmat)

Jika Anda memiliki banyak objek dan ingin menulis objek rumus dengan o1+o2+...+oncepat, Anda dapat menggunakannya

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

untuk menghasilkan rumus untuk gnm(yang tidak Anda perlukan llbt.fit).

Ada artikel JSS , lihat juga https://r-forge.r-project.org/projects/prefmod/ dan dokumentasinya via ?llbt.design.

Momo
sumber
1
Itu respons yang sangat menyeluruh. Terima kasih. Sepertinya prefmod akan menjadi paket yang baik untuk digunakan. Saya tertarik menggunakan model untuk mencoba memprediksi hasil pertandingan olahraga.
Silverfish
Tidak masalah, senang jika itu membantu. Saya tidak tahu persis bagaimana Anda bermaksud memprediksi, tetapi Leitner et al. telah menggunakan model ini untuk memprediksi acara olahraga. Lihat tesisnya epubdev.wu.ac.at/2925 . Semoga berhasil.
Momo
Mungkin tautan ini lebih baik epubdev.wu.ac.at/view/creators/…
Momo
Apakah mungkin untuk menghitung signifikansi untuk perbedaan antara pasangan individu (misalnya, o1 dan o2) dari data ini? Atau apakah Anda harus mengatur ulang rumus, menggunakan o2 sebagai faktor terakhir dan hidup tanpa estimasi Std.error dalam kasus itu?
TNT
1
Sudah lama, jadi saya tidak ingat apakah Anda dapat dengan mudah menggunakan batasan linear tetapi apa yang dapat Anda lakukan dalam kasus Anda adalah menggunakan satu sebagai tingkat referensi, katakanlah o1, dan gunakan nilai t yang lain, katakanlah o2, dari ringkasan - ini secara efektif merupakan tes apakah perbedaan antara o1 dan o2 adalah nol.
Momo