Analysis of pudding following Davidson 1970
R. K. S. Hankin
To cite the hyper2 package in publications, please use Hankin (2017).
This short document creates the pudding suite of R objects and
presents some preliminary analysis. It is not particularly
well-structured or terse and is essentially a sequence of notes. It
follows Davidson 1970 and the PlackettLuce package.
Here is the pudding result:
pudding_table <- pudding # standard hyper2 terminology
rownames(pudding_table) <- apply(allbinom(6,2),2,paste,collapse="v")
pudding_table## i j r_ij w_ij w_ji t_ij
## 1v2 1 2 57 19 22 16
## 1v3 1 3 47 16 19 12
## 1v4 2 3 48 19 19 10
## 1v5 1 4 54 18 23 13
## 1v6 2 4 51 23 19 9
## 2v3 3 4 54 19 20 15
## 2v4 1 5 50 13 19 18
## 2v5 2 5 48 16 20 12
## 2v6 3 5 48 16 15 17
## 3v4 4 5 47 17 14 16
## 3v5 1 6 51 18 21 12
## 3v6 2 6 54 22 20 12
## 4v5 3 6 41 13 18 10
## 4v6 4 6 51 14 19 18
## 5v6 5 6 44 11 21 12Here is a method to convert the pudding matrix to a hyper2
likelihood function:
numbers <- as.character(1:6)
puddings <- paste("pudding",numbers,sep="_")
tyers <- paste("tie",numbers,sep="_")
`convert` <- function(M){
H <- hyper2()
I <- hyper2()
J <- hyper2()
for(i in seq_len(nrow(M))){
r_ij <- M[i,3]
w_ij <- M[i,4]
w_ji <- M[i,5]
t_ij <- M[i,6]
H[c(puddings[c(M[i,1] )] )] %<>% inc(w_ij)
H[c(puddings[c( M[i,2])] )] %<>% inc(w_ji)
H[c( "tie")] %<>% inc(t_ij)
H[c(puddings[c(M[i,1],M[i,2])],"tie")] %<>% dec(r_ij)
I[c(puddings[c(M[i,1] )] )] %<>% inc(w_ij)
I[c(puddings[c( M[i,2])] )] %<>% inc(w_ji)
I[c(tyers [c(M[i,1],M[i,2])] )] %<>% inc(t_ij)
I[c(puddings[c(M[i,1],M[i,2])],tyers[c(M[i,1],M[i,2])])] %<>% dec(r_ij)
J[c(puddings[c(M[i,1] )],"first")] %<>% inc(w_ij)
J[c(puddings[c( M[i,2])] )] %<>% inc(w_ji)
J[c(tyers [c(M[i,1],M[i,2])] )] %<>% inc(t_ij)
J[c(puddings[c(M[i,1],M[i,2])],tyers[c(M[i,1],M[i,2])],"first")] %<>% dec(r_ij)
}
return(list(H,I,J))
}Function convert() gives a list of three hyper2 objects. The
first assumes that each pudding has a strength \(p_i\), \(1\leq i\leq 6\)
and
\[ \operatorname{Prob}(i\succ j) = \frac{\pi_i}{\pi_i+\pi_j+D}\\ \operatorname{Prob}(j\succ i) = \frac{\pi_j}{\pi_i+\pi_j+D}\\ \operatorname{Prob}(i\sim j) = \frac{D}{\pi_i+\pi_j+D} \]
where \(p_i,D\geq 0\) and \(\sum p_i+D=1\) (here \(i\succ j\) means \(i\) is preferred to \(j\) and \(i\sim j\) means no preference). The second element is a generalization of the first, allowing the tying propensities to differ:
\[ \operatorname{Prob}(i\succ j) = \frac{\pi_i}{\pi_i+\pi_j+D_i+D_j}\\ \operatorname{Prob}(j\succ i) = \frac{\pi_j}{\pi_i+\pi_j+D_i+D_j}\\ \operatorname{Prob}(i\sim j) = \frac{D}{\pi_i+\pi_j+D_i+D_j} \]
where \(\sum p_i+D_i=1\). The third element is a generalization of the second, giving the first pudding a “first-mover” advantage analogous to the home-team advantage in football, or playing white in chess:
\[ \operatorname{Prob}(i\succ j) = \frac{\pi_i+F}{\pi_i+\pi_j+D_i+D_j+F}\\ \operatorname{Prob}(j\succ i) = \frac{\pi_j}{\pi_i+\pi_j+D_i+D_j+F}\\ \operatorname{Prob}(i\sim j) = \frac{D}{\pi_i+\pi_j+D_i+D_j+F} \]
Above, the first-mover advantage is given to pudding \(i\). We can use these now:
lf <- convert(pudding_table)
pudding <- lf[[1]] # standard hyper2 terminology
pudding2 <- lf[[2]]
pudding3 <- lf[[3]]
pudding_maxp <- maxp(pudding)
pudding2_maxp <- maxp(pudding2)
pudding3_maxp <- maxp(pudding3)pudding_maxp## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.12598 0.16609 0.14549 0.14614 0.13023 0.17619 0.10988pudding2_maxp## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## 0.101583 0.126039 0.117246 0.119040 0.115048 0.147511 0.041762 0.024338 0.041120 0.044663 0.068661 0.052990pudding3_maxp## first pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3
## 1.1075e-06 1.0162e-01 1.2597e-01 1.1725e-01 1.1905e-01 1.1504e-01 1.4756e-01 4.1777e-02 2.4305e-02 4.1121e-02
## tie_4 tie_5 tie_6
## 4.4675e-02 6.8674e-02 5.2955e-021 Pudding 1
We will analyse pudding first:
summary(pudding)## A hyper2 object of size 7.
## pnames: pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## Number of brackets: 22
## Sum of powers: 0
##
## Table of bracket lengths:
## 1 3
## 7 15
##
## Table of powers:
## -57 -54 -51 -50 -48 -47 -44 -41 79 84 86 93 99 102 202
## 1 3 3 1 3 2 1 1 1 1 1 1 1 1 1pie(pudding_maxp)
First we use pudding to test a null of equal pudding strengths:
samep.test(pudding,puddings)## Warning in optim(theta.old, fun, gradient, control = control, method = method, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly##
## Constrained support maximization
##
## data: pudding
## null hypothesis: pudding_1 = pudding_2 = pudding_3 = pudding_4 = pudding_5 = pudding_6
## null estimate:
## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.14828 0.14828 0.14828 0.14828 0.14828 0.14828 0.11031
## (argmax, constrained optimization)
## Support for null: -811.75 + K
##
## alternative hypothesis: sum p_i=1
## alternative estimate:
## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.12598 0.16609 0.14549 0.14614 0.13023 0.17619 0.10988
## (argmax, free optimization)
## Support for alternative: -809 + K
##
## degrees of freedom: 5
## support difference = 2.7499
## p-value: 0.35797so according to samep.test() the puddings are all the same strength.
We note that the propensity to tie cannot possibly be zero but we can
give a profile likelihood curve:
ptie <- seq(from=0.08,len=20,to=0.15)
l <- profsupp(pudding,"tie",ptie)plot(ptie,l,type='b',lwd=2)
abline(h=c(0,-2))
#abline(v=c(0.095,0.127))
grid()
and we see that a credible interval for the proclivity to tie would be \((0.095,0.127)\).
2 Pudding 2
Object pudding2 includes a pudding-specific propensity to tie. We
can test different nulls:
samep.test(pudding2,puddings)##
## Constrained support maximization
##
## data: pudding2
## null hypothesis: pudding_1 = pudding_2 = pudding_3 = pudding_4 = pudding_5 = pudding_6
## null estimate:
## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## 0.120920 0.120920 0.120920 0.120920 0.120920 0.120920 0.048670 0.022195 0.042442 0.046405 0.070684 0.044083
## (argmax, constrained optimization)
## Support for null: -808.78 + K
##
## alternative hypothesis: sum p_i=1
## alternative estimate:
## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## 0.101583 0.126039 0.117246 0.119040 0.115048 0.147511 0.041762 0.024338 0.041120 0.044663 0.068661 0.052990
## (argmax, free optimization)
## Support for alternative: -806.74 + K
##
## degrees of freedom: 5
## support difference = 2.0416
## p-value: 0.53749samep.test(pudding2,tyers)##
## Constrained support maximization
##
## data: pudding2
## null hypothesis: tie_1 = tie_2 = tie_3 = tie_4 = tie_5 = tie_6
## null estimate:
## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## 0.103292 0.136139 0.119258 0.119863 0.106796 0.144432 0.045037 0.045037 0.045037 0.045037 0.045037 0.045037
## (argmax, constrained optimization)
## Support for null: -809 + K
##
## alternative hypothesis: sum p_i=1
## alternative estimate:
## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## 0.101583 0.126039 0.117246 0.119040 0.115048 0.147511 0.041762 0.024338 0.041120 0.044663 0.068661 0.052990
## (argmax, free optimization)
## Support for alternative: -806.74 + K
##
## degrees of freedom: 5
## support difference = 2.2594
## p-value: 0.477363 Pudding 3
Object pudding3 includes a first-mover advantage:
summary(pudding3)## A hyper2 object of size 13.
## pnames: first pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## Number of brackets: 40
## Sum of powers: 0
##
## Table of bracket lengths:
## 1 2 5
## 5 20 15
##
## Table of powers:
## -57 -54 -51 -50 -48 -47 -44 -41 9 10 11 12 13 15 16 17 18 22 31 38 48 62 68 80 84 99
## 1 3 3 1 3 2 1 1 1 2 1 5 1 1 2 1 2 1 1 1 1 1 1 1 1 1pudding3_maxp## first pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3
## 1.1075e-06 1.0162e-01 1.2597e-01 1.1725e-01 1.1905e-01 1.1504e-01 1.4756e-01 4.1777e-02 2.4305e-02 4.1121e-02
## tie_4 tie_5 tie_6
## 4.4675e-02 6.8674e-02 5.2955e-02We can see above that the first-mover advantage is small (and indeed a “second-mover advantage”) is also very small.
4 Davidson’s analysis
Davidson (1970) considers the following probability model:
\[ \operatorname{Prob}(i\succ j) = \frac{\pi_i}{\pi_i+\pi_j+\nu\sqrt{\pi_i\pi_j}}\\ \operatorname{Prob}(j\succ i) = \frac{\pi_j}{\pi_i+\pi_j+\nu\sqrt{\pi_i\pi_j}}\\ \operatorname{Prob}(i\sim j) = \frac{\nu\sqrt{\pi_i\pi_j}}{\pi_i+\pi_j+\nu\sqrt{\pi_i\pi_j}} \]
and gives the following estimate for the pudding data:
davidson <- c(0.1388005,0.1729985, 0.1617420, 0.1653930, 0.1586805, 0.2023855, 0.7468147)
names(davidson) <- c(puddings,"tie")
davidson## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.13880 0.17300 0.16174 0.16539 0.15868 0.20239 0.74681par(pty='s')
plot(davidson[1:6],pudding_maxp[1:6]/(1-pudding_maxp[7]),asp=1,xlim=c(0.13,0.2),ylim=c(0.13,0.2),pch=16)
abline(0,1)
Although the probability model is different, we may compare the probabilities of \(i\succ j\), \(j\succ i\), and \(i\sim j\) for \(i<j\) with the two models. This gives us a total of \(3{6\choose 2}=45\) (nonindependent) probabilities and we may plot a scattergraph:
ph <- pudding_maxp # ph == 'pudding hankin'; saves typing
pd <- davidson # pd == 'pudding davidson'
print(ph)## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.12598 0.16609 0.14549 0.14614 0.13023 0.17619 0.10988print(pd)## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.13880 0.17300 0.16174 0.16539 0.15868 0.20239 0.74681`hankin` <- function(ph){
## ph a vector of length n+1; first n entries are strengths, last one
## propensity to tie
n <- length(ph)-1
m <- allbinom(n,2)
out <- matrix(NA,nrow=ncol(m),ncol=3)
colnames(out) <- c("i","j","tie")
rownames(out) <- apply(m,2,paste,collapse="v")
th <- ph[n+1]
for(o in seq_len(ncol(m))){
i <- m[1,o]
j <- m[2,o]
dh <- ph[i] + ph[j] + th # dh == 'denominator hankin'
out[o,1] <- ph[i]/dh
out[o,2] <- ph[j]/dh
out[o,3] <- th/dh
}
return(out)
}
`davidson` <- function(pd){
## pd a vector of length n+1; first n entries are strengths, last one
## nu
n <- length(pd)-1 # how many puddings
nu <- pd[n+1]
m <- allbinom(n,2)
out <- matrix(NA,nrow=ncol(m),ncol=3)
colnames(out) <- c("i","j","tie")
rownames(out) <- apply(m,2,paste,collapse="v")
nu <- pd[n+1]
for(o in seq_len(ncol(m))){
i <- m[1,o]
j <- m[2,o]
dd <- pd[i] + pd[j] + nu*sqrt(pd[i]*pd[j]) # dd == 'denominator davidson'
out[o,1] <- pd[i]/dd
out[o,2] <- pd[j]/dd
out[o,3] <- nu*sqrt(pd[i]*pd[j])/dd
}
return(out)
}
plotter <- function(ph,pd,...){
h <- hankin(ph)
d <- davidson(pd)
nh <- nrow(h)
plot(c(h),c(d),pch=16,asp=1,col=c(rep("red",nh),rep("blue",nh),rep("green",nh)),...)
}print(ph)## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.12598 0.16609 0.14549 0.14614 0.13023 0.17619 0.10988print(hankin(ph))## i j tie
## 1v2 0.31342 0.41322 0.27336
## 1v3 0.33036 0.38151 0.28813
## 1v4 0.32979 0.38258 0.28763
## 1v5 0.34412 0.35574 0.30013
## 1v6 0.30574 0.42759 0.26666
## 2v3 0.39409 0.34520 0.26071
## 2v4 0.39348 0.34622 0.26030
## 2v5 0.40889 0.32061 0.27050
## 2v6 0.36733 0.38966 0.24300
## 3v4 0.36235 0.36399 0.27366
## 3v5 0.37730 0.33775 0.28495
## 3v6 0.33712 0.40827 0.25461
## 4v5 0.37836 0.33717 0.28447
## 4v6 0.33813 0.40765 0.25422
## 5v6 0.31284 0.42323 0.26394print(pd)## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.13880 0.17300 0.16174 0.16539 0.15868 0.20239 0.74681print(davidson(pd))## i j tie
## 1v2 0.32466 0.40465 0.27069
## 1v3 0.33654 0.39216 0.27131
## 1v4 0.33258 0.39630 0.27113
## 1v5 0.33994 0.38862 0.27144
## 1v6 0.29763 0.43397 0.26840
## 2v3 0.37636 0.35187 0.27177
## 2v4 0.37227 0.35590 0.27183
## 2v5 0.37987 0.34843 0.27170
## 2v6 0.33584 0.39289 0.27128
## 3v4 0.36000 0.36813 0.27187
## 3v5 0.36754 0.36058 0.27187
## 3v6 0.32397 0.40538 0.27064
## 4v5 0.37162 0.35654 0.27184
## 4v6 0.32789 0.40123 0.27088
## 5v6 0.32063 0.40894 0.27043jj <- c(2,5)/10
par(pty='s')
plotter(ph,pd,xlim=jj,ylim=jj)
grid()
Now calculate likelihoods:
n <- pudding_table[,-(1:3)]
h <- hankin(ph)
d <- davidson(pd)
showsome(n)## w_ij w_ji t_ij
## 1v2 19 22 16
## 1v3 16 19 12
## 1v4 19 19 10
## ... ... ... ...
## 5v6 11 21 12showsome(h)## i j tie
## 1v2 0.3134 0.4132 0.2734
## 1v3 0.3304 0.3815 0.2881
## 1v4 0.3298 0.3826 0.2876
## ... ... ... ...
## 5v6 0.3128 0.4232 0.2639showsome(d)## i j tie
## 1v2 0.3247 0.4047 0.2707
## 1v3 0.3365 0.3922 0.2713
## 1v4 0.3326 0.3963 0.2711
## ... ... ... ...
## 5v6 0.3206 0.4089 0.2704sum(n*log(h))## [1] -815.78sum(n*log(d))## [1] -812.78(suppdiff <- sum(n*log(d))-sum(n*log(h)))## [1] 2.9959Above we see that Davidson’s method has almost 3 units of support more than RBT. However, we note that Davidson’s model has seven degrees of freedom while RBT has only six (because of the the unit sum constraint). Here Wilks’s theorem is not applicable because the two models are not nested but even so the likelihood ratio is a meaningful statistic.
pchisq(suppdiff,df=1,lower.tail=FALSE)## [1] 0.0834775 Appendix: synthetic datasets
Above we see that, according to this model, there is no evidence to suggest that the puddings’ strengths, or their tying ability, differ. But we can adduce a dataset for which the tying ability does differ. Consider this:
p <- pudding_table
p[1 ,4:6] <- c(10,10,37) # 1v2
p[2 ,4:6] <- c(10,10,27) # 1v3
p[4 ,4:6] <- c(10,10,34) # 1v4
p[7 ,4:6] <- c(20,20,10) # 1v5
p[11,4:6] <- c(10,10,31) # 1v6
p## i j r_ij w_ij w_ji t_ij
## 1v2 1 2 57 10 10 37
## 1v3 1 3 47 10 10 27
## 1v4 2 3 48 19 19 10
## 1v5 1 4 54 10 10 34
## 1v6 2 4 51 23 19 9
## 2v3 3 4 54 19 20 15
## 2v4 1 5 50 20 20 10
## 2v5 2 5 48 16 20 12
## 2v6 3 5 48 16 15 17
## 3v4 4 5 47 17 14 16
## 3v5 1 6 51 10 10 31
## 3v6 2 6 54 22 20 12
## 4v5 3 6 41 13 18 10
## 4v6 4 6 51 14 19 18
## 5v6 5 6 44 11 21 12(above, we doctor the dataset so that any comparison involving pudding 1 has a high probability of tying). Then:
Hp <- convert(p)[[2]]
options(use_alabama = TRUE)
mHp <- maxp(Hp)
mHpt <- mHp[7:12] #just the tyersvisualise:
mHpt## tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## 0.202526 0.014232 0.038001 0.042226 0.043104 0.050527plot(mHpt)
(above we see that pudding 1 does indeed have a higher tie probability). Test the null of equal tying probabilities:
(tHp <- samep.test(Hp,tyers))##
## Constrained support maximization
##
## data: Hp
## null hypothesis: tie_1 = tie_2 = tie_3 = tie_4 = tie_5 = tie_6
## null estimate:
## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## 0.070468 0.119096 0.107462 0.102020 0.113834 0.127296 0.059971 0.059971 0.059971 0.059971 0.059971 0.059971
## (argmax, constrained optimization)
## Support for null: -810.93 + K
##
## alternative hypothesis: sum p_i=1
## alternative estimate:
## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie_1 tie_2 tie_3 tie_4 tie_5 tie_6
## 0.102552 0.098235 0.094684 0.092287 0.103188 0.118439 0.202526 0.014232 0.038001 0.042226 0.043104 0.050527
## (argmax, free optimization)
## Support for alternative: -788.54 + K
##
## degrees of freedom: 5
## support difference = 22.39
## p-value: 1.6082e-08Above analysis shows that we can reject the null that all puddings have the same tying propensity.
5.1 Test of Davidson vs Hankin
Although the likelihood for Davidson’s model was higher than mine for the original dataset, we can create a synthetic dataset which is a random observation drawn from a reified Bradley-Terry probability model. I will then test the hypothesis that the data is in fact drawn from an RBT distribution by calculating a likelihood ratio.
set.seed(0)
f <- function(i){rmultinom(1,50,prob=h[i,])} # f(i) = sample from row i of h
synth <- pudding_table
synth[,-(1:3)] <- t(sapply(1:15,f)) # random sample
synth[,3] <- 50
rownames(synth) <- rownames(p)
colnames(synth) <- colnames(p)
showsome(synth)## i j r_ij w_ij w_ji t_ij
## 1v2 1 2 50 20 20 10
## 1v3 1 3 50 15 19 16
## 1v4 2 3 50 21 19 10
## ... ... ... ... ... ... ...
## 5v6 5 6 50 15 18 17Above, synth is a dataset drawn from an RBT model with maximum
likelihood parameters obtained from the original dataset. We now use
both methods. First reified Bradley-Terry:
(H_synth <- convert(synth)[[1]])## log(pudding_1^94 * (pudding_1 + pudding_2 + tie)^-50 * (pudding_1 + pudding_3 + tie)^-50 * (pudding_1 +
## pudding_4 + tie)^-50 * (pudding_1 + pudding_5 + tie)^-50 * (pudding_1 + pudding_6 + tie)^-50 * pudding_2^91
## * (pudding_2 + pudding_3 + tie)^-50 * (pudding_2 + pudding_4 + tie)^-50 * (pudding_2 + pudding_5 + tie)^-50
## * (pudding_2 + pudding_6 + tie)^-50 * pudding_3^85 * (pudding_3 + pudding_4 + tie)^-50 * (pudding_3 +
## pudding_5 + tie)^-50 * (pudding_3 + pudding_6 + tie)^-50 * pudding_4^95 * (pudding_4 + pudding_5 + tie)^-50
## * (pudding_4 + pudding_6 + tie)^-50 * pudding_5^80 * (pudding_5 + pudding_6 + tie)^-50 * pudding_6^99 *
## tie^206)ph_synth <- maxp(H_synth)ph_synth## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.15534 0.14832 0.13492 0.15772 0.12436 0.16754 0.11180The above strengths are different from (but not too different from)
the value obtained from maxp(H) above. Now analyse the same dataset
but using the method given in the PlackettLuce vignette:
PlackettLuce_coef <- function(pudding_table){
i_wins <- data.frame(Winner = pudding_table$i, Loser = pudding_table$j)
j_wins <- data.frame(Winner = pudding_table$j, Loser = pudding_table$i)
ties <- data.frame(Winner = asplit(pudding_table[c("i", "j")], 1),
Loser = rep(NA, 15))
R <- as.rankings(rbind(i_wins, j_wins, ties),
input = "orderings")
w <- unlist(pudding_table[c("w_ij", "w_ji", "t_ij")])
mod <- PlackettLuce(R, weights = w, npseudo = 0, maxit = 70)
return(coef(mod,log=FALSE))
} As a consistency check, try verifying the original result:
PlackettLuce_coef(pudding_table)## 1 2 3 4 5 6 tie2
## 0.13880 0.17300 0.16175 0.16537 0.15869 0.20239 0.74682Above, we see values that match those in the PlackettLuce vignette.
Now apply the same method to the synthetic dataset synth:
(pd_synth <- PlackettLuce_coef(synth))## 1 2 3 4 5 6 tie2
## 0.18832 0.16098 0.14679 0.16703 0.14146 0.19542 0.76054Now we need to calculate the probabilities:
showsome(h_synth <- hankin(ph_synth))## i j tie
## 1v2 0.3739 0.357 0.2691
## 1v3 0.3864 0.3356 0.2781
## 1v4 0.3656 0.3712 0.2631
## ... ... ... ...
## 5v6 0.308 0.415 0.2769showsome(d_synth <- davidson(pd_synth))## i j tie
## 1v2 0.3909 0.3342 0.2749
## 1v3 0.408 0.318 0.274
## 1v4 0.3841 0.3407 0.2751
## ... ... ... ...
## 5v6 0.3053 0.4218 0.2729And calculate a likelihood ratio for the two hypotheses, RBT and Davidson:
showsome(n <- synth[,-(1:3),])## w_ij w_ji t_ij
## 1v2 20 20 10
## 1v3 15 19 16
## 1v4 21 19 10
## ... ... ... ...
## 5v6 15 18 17(lh_synth <- sum(n*log(h_synth))) # Hankin; RBT## [1] -816.95(ld_synth <- sum(n*log(d_synth))) # Davidson## [1] -817.89lh_synth - ld_synth## [1] 0.93563So reified Bradley-Terry is better than Davidson, but (in this case at least) not by much. The difference is small because the puddings have very similar strengths. We can repeat the analysis but with puddings of different strengths:
dpudstrengths <- c(zipf(6)*0.9,tie=0.1) # different pudding strengths
names(dpudstrengths) <- c(puddings,"tie")
dpudstrengths## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.367347 0.183673 0.122449 0.091837 0.073469 0.061224 0.100000Above we see that pudding_1 is six times stronger than pudding_6
in accordance with Zipf’s law. Same analysis as before:
showsome(h2 <- hankin(dpudstrengths))## i j tie
## 1v2 0.5643 0.2821 0.1536
## 1v3 0.6228 0.2076 0.1696
## 1v4 0.6569 0.1642 0.1788
## ... ... ... ...
## 5v6 0.313 0.2609 0.4261f <- function(i){rmultinom(1,50,prob=h2[i,])}
synth2 <- pudding_table
synth2[,-(1:3)] <- t(sapply(1:15,f)) # random sample
synth2[,3] <- 50
rownames(synth2) <- rownames(p)
colnames(synth2) <- colnames(p)
showsome(synth2)## i j r_ij w_ij w_ji t_ij
## 1v2 1 2 50 30 13 7
## 1v3 1 3 50 30 11 9
## 1v4 2 3 50 36 8 6
## ... ... ... ... ... ... ...
## 5v6 5 6 50 16 14 20OK so synth2 is a different synthetic dataset but this time with
pudding strengths following Zipf. As before, create a hyper2
object and find the evaluate:
H_synth2 <- convert(synth2)[[1]] # a hyper2 likelihood function
ph_synth2 <- maxp(H_synth2)ph_synth2## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## 0.281754 0.269806 0.112006 0.077407 0.080582 0.074286 0.104159As a consistency check compare estimate with true value:
ph_synth2-dpudstrengths # difference is small## pudding_1 pudding_2 pudding_3 pudding_4 pudding_5 pudding_6 tie
## -0.0855925 0.0861323 -0.0104431 -0.0144299 0.0071127 0.0130612 0.0041593Now Davidson estimate using synth2:
(pd_synth2 <- PlackettLuce_coef(synth2))## 1 2 3 4 5 6 tie2
## 0.324795 0.297882 0.119001 0.073000 0.087016 0.098306 0.861888h_synth2 <- hankin(ph_synth2)
d_synth2 <- davidson(pd_synth2)
showsome(h_synth2) # same as before## i j tie
## 1v2 0.4297 0.4115 0.1588
## 1v3 0.5659 0.2249 0.2092
## 1v4 0.6081 0.1671 0.2248
## ... ... ... ...
## 5v6 0.3111 0.2868 0.4021And calculate a support function:
n <- synth2[,-(1:3),] # as before, extract just the counts
(lh_synth2 <- sum(n*log(h_synth2)))## [1] -754.96(ld_synth2 <- sum(n*log(d_synth2)))## [1] -770.12lh_synth2 - ld_synth2## [1] 15.161So this time the difference in support is larger, due to the puddings having different strengths.
Package dataset
Following lines create pudding.rda, residing in the data/ directory of the package.
save(pudding_table, pudding, pudding_maxp, file = "pudding.rda")References
- R. D. Davidson 1970. “On extending the Bradley-Terry model to accommodate ties in paired comparison experiments”. Journal of the American Statistical Association, Volume 65, Number 329, pp317–328
hyper2 Package: Likelihood Functions for Generalized Bradley-Terry Models.” The R Journal 9 (2): 429–39.