Ladies’ figure skating in the 2002 Winter Olympics
R. K. S. Hankin
(takes about twenty minutes to run without cache)
To cite the hyper2 package in publications, please use
Hankin (2017). This document analyses Bradley-Terry strengths of
the competitors in the Ladies’ Figure Skating at the 2002 Winter
Olympics. In the package, dataframe skating_table is copied from
Lock and Lock -Lock, Robin and Lock, Kari Frazer (2003). It also creates file data/skating.rda
which contains R objects skating, skating_table, and
skating_maxp.
set.seed(0)
skating_table <- ordertable(read.table("skating.txt", header=TRUE))
skating_table## An ordertable:
## J1 J2 J3 J4 J5 J6 J7 J8 J9
## hughes 1 4 3 4 1 2 1 1 1
## slutskaya 3 1 1 1 4 1 2 3 2
## kwan 2 3 2 2 2 3 3 2 3
## cohen 5 2 4 3 3 4 4 4 4
## suguri 4 8 5 5 5 7 5 5 5
## butyrskaya 6 5 8 7 12 5 8 7 6
## robinson 7 7 7 9 6 8 10 6 7
## sebestyen 8 10 12 8 7 6 12 8 8
## kettunen 9 9 13 6 13 10 7 11 14
## volchkova 10 6 14 11 10 12 6 9 15
## maniachenko 13 12 11 12 16 11 11 10 9
## fontana 14 11 18 16 9 15 9 12 10
## liashenko 15 13 6 10 8 14 13 14 16
## onda 11 14 10 15 15 13 15 13 11
## hubert 12 17 17 13 11 16 14 15 13
## meier 16 16 9 14 14 9 16 16 12
## gusmeroli 17 15 15 17 17 18 17 17 17
## soldatova 19 18 22 20 21 17 18 18 19
## hegel 20 21 16 22 18 19 21 19 18
## giunchi 18 19 20 21 19 20 20 20 20
## babiakova 22 20 19 19 20 21 19 22 22
## kopac 21 22 23 18 22 22 22 21 21
## luca 23 23 21 23 23 23 23 23 23(we note in passing that the Plackett-Luce strengths of every competitors is strictly positive. We can see this by observing that:
\[\begin{alignat*}{4} \textbf{luka} &\succ_3\mbox{kopak} &&\succ_1\mbox{babiakova} &&\succ_3\mbox{giunchi} &&\succ_1\mbox{hegel} \\ &\succ_3\mbox{soldatova} &&\succ_6\mbox{gusmeroli} && \succ_2\mbox{meier} &&\succ_3\mbox{hubert} \\ &\succ_4\mbox{onda} &&\succ_1\mbox{kiashenko} &&\succ_3\mbox{fontana} &&\succ_2\mbox{maniachenko}\\ &\succ_3\mbox{volchkova} &&\succ_2\mbox{kettunen} &&\succ_2\mbox{sebestyen} &&\succ_4\mbox{robinson} \\ &\succ_3\mbox{butyrskaya} &&\succ_2\mbox{suguri} &&\succ_1\mbox{cohen} &&\succ_2\mbox{kwan} \\ &\succ_1\mbox{slutskaya} &&\succ_2\mbox{hughes} \\ &\succ_1\textbf{luka} \end{alignat*}\]
and if any competitor has zero BT strength, at least one of the above successors is false with probability one).
Object skating_table is an order table. It is structured so that
each competitor is a row, and each judge is a column. One good thing
to do is to present the dataset as a rank table:
as.ranktable(skating_table)## A ranktable:
## c1 c2 c3 c4 c5 c6 c7 c8 c9
## J1 hughes kwan slutskaya suguri cohen butyrskaya robinson sebestyen kettunen
## J2 slutskaya cohen kwan hughes butyrskaya volchkova robinson suguri kettunen
## J3 slutskaya kwan hughes cohen suguri liashenko robinson butyrskaya meier
## J4 slutskaya kwan cohen hughes suguri kettunen butyrskaya sebestyen robinson
## J5 hughes kwan cohen slutskaya suguri robinson sebestyen liashenko fontana
## J6 slutskaya hughes kwan cohen butyrskaya sebestyen suguri robinson meier
## J7 hughes slutskaya kwan cohen suguri volchkova kettunen butyrskaya fontana
## J8 hughes kwan slutskaya cohen suguri robinson butyrskaya sebestyen volchkova
## J9 hughes slutskaya kwan cohen suguri butyrskaya robinson sebestyen maniachenko
## c10 c11 c12 c13 c14 c15 c16 c17
## J1 volchkova onda hubert maniachenko fontana liashenko meier gusmeroli
## J2 sebestyen fontana maniachenko liashenko onda gusmeroli meier hubert
## J3 onda maniachenko sebestyen kettunen volchkova gusmeroli hegel hubert
## J4 liashenko volchkova maniachenko hubert meier onda fontana gusmeroli
## J5 volchkova hubert butyrskaya kettunen meier onda maniachenko gusmeroli
## J6 kettunen maniachenko volchkova onda liashenko fontana hubert soldatova
## J7 robinson maniachenko sebestyen liashenko hubert onda meier gusmeroli
## J8 maniachenko kettunen fontana onda liashenko hubert meier gusmeroli
## J9 fontana onda meier hubert kettunen volchkova liashenko gusmeroli
## c18 c19 c20 c21 c22 c23
## J1 giunchi soldatova hegel kopac babiakova luca
## J2 soldatova giunchi babiakova hegel kopac luca
## J3 fontana babiakova giunchi luca soldatova kopac
## J4 kopac babiakova soldatova giunchi hegel luca
## J5 hegel giunchi babiakova soldatova kopac luca
## J6 gusmeroli hegel giunchi babiakova kopac luca
## J7 soldatova babiakova giunchi hegel kopac luca
## J8 soldatova hegel giunchi kopac babiakova luca
## J9 hegel soldatova giunchi kopac babiakova lucaIn the above, the column names should be read “came first” (c1),
“came second” (c2) and so on to “came 23rd” (c23). From the first
column, we see that Hughes came first 5 times [according to judges
1,5,7,8,9] and Slutskya came first four times [judges 2,3,4,6]. From
the last column we see that all judges except J3 awarded Luca last
place. File man/rrank.Rd has a detailed discussion of the
differences between rank and order.
We can convert an order table to a support function:
summary(suppfun(skating_table))## A hyper2 object of size 23.
## pnames: babiakova butyrskaya cohen fontana giunchi gusmeroli hegel hubert hughes kettunen kopac kwan liashenko luca maniachenko meier onda robinson sebestyen slutskaya soldatova suguri volchkova
## Number of brackets: 135
## Sum of powers: 0
##
## Table of bracket lengths:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## 23 4 4 6 4 3 3 7 7 8 7 9 8 8 6 7 7 2 2 3 4 2 1
##
## Table of powers:
## -9 -8 -7 -6 -5 -4 -3 -2 -1 1 8 9
## 1 1 3 1 1 5 8 13 79 1 1 21Further, we might ask how many judges ranked each competitor first, second and so on:
oo <- sapply(seq_len(23),function(i){rowSums(skating_table==i)})
colnames(oo) <- rep(" ",23) # nicer print
oo##
## hughes 5 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## slutskaya 4 2 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## kwan 0 5 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## cohen 0 1 2 5 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## suguri 0 0 0 1 6 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## butyrskaya 0 0 0 0 2 2 2 2 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## robinson 0 0 0 0 0 2 4 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## sebestyen 0 0 0 0 0 1 1 4 0 1 0 2 0 0 0 0 0 0 0 0 0 0 0
## kettunen 0 0 0 0 0 1 1 0 2 1 1 0 2 1 0 0 0 0 0 0 0 0 0
## volchkova 0 0 0 0 0 2 0 0 1 2 1 1 0 1 1 0 0 0 0 0 0 0 0
## maniachenko 0 0 0 0 0 0 0 0 1 1 3 2 1 0 0 1 0 0 0 0 0 0 0
## fontana 0 0 0 0 0 0 0 0 2 1 1 1 0 1 1 1 0 1 0 0 0 0 0
## liashenko 0 0 0 0 0 1 0 1 0 1 0 0 2 2 1 1 0 0 0 0 0 0 0
## onda 0 0 0 0 0 0 0 0 0 1 2 0 2 1 3 0 0 0 0 0 0 0 0
## hubert 0 0 0 0 0 0 0 0 0 0 1 1 2 1 1 1 2 0 0 0 0 0 0
## meier 0 0 0 0 0 0 0 0 2 0 0 1 0 2 0 4 0 0 0 0 0 0 0
## gusmeroli 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 6 1 0 0 0 0 0
## soldatova 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 3 2 1 1 1 0
## hegel 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 2 2 1 2 1 0
## giunchi 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 5 1 0 0
## babiakova 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 2 1 3 0
## kopac 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 3 4 1
## luca 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 8In the above, each row is a competitor and each column corresponds to a ranking. Thus the first column corresponds to “first place” and shows that five judges ranked Hughes first and four judges ranked Slutskaya first. The second column corresponds to “second place” and shows that one judge gave Hughes second place, two judges placed Slutskaya second, five gave Kawn second, and one gave Cohen second. The first row corresponds to Hughes and we see that, of \(5+1+1+2=9\) judges, 5 placed Hughes first, one placed her second, one placed her third, and two placed her fourth.
The formal ordering used in competition is, according to Lock and Lock, given by the median ordinal:
apply(skating_table,1,median)## hughes slutskaya kwan cohen suguri butyrskaya robinson sebestyen
## 1 2 2 4 5 7 7 8
## kettunen volchkova maniachenko fontana liashenko onda hubert meier
## 10 10 11 12 13 13 14 14
## gusmeroli soldatova hegel giunchi babiakova kopac luca
## 17 19 19 20 20 22 23but with ties broken using a complicated hierarchy, the first of which is the “size of the majority”:
rowSums(sweep(skating_table,1,seq_len(23))>=0)## hughes slutskaya kwan cohen suguri butyrskaya robinson sebestyen
## 9 5 4 6 8 7 7 7
## kettunen volchkova maniachenko fontana liashenko onda hubert meier
## 7 6 7 5 6 4 4 4
## gusmeroli soldatova hegel giunchi babiakova kopac luca
## 7 8 6 6 4 5 8which would suggest that Slutskaya beats Kwan (5-4), Butyrskaya draws
with Robinson (7-7), and so on. The rows of skating_table are in
the order given by this system.
1 Data visualization
Now some data visualization. First the MLE for the strengths:
skating <- suppfun(skating_table)
options("use_alabama" = TRUE)
skating_maxp <- maxp(skating,startp=
c(5.0493913251e-06, 0.0053037393083, 0.090240695999,
0.00036893464900, 1.1397045995e-05, 8.3424027249e-05,
9.8129903007e-06, 0.00026960143514, 0.29831585124,
0.001376333090, 2.9434081237e-06, 0.26550514587,
0.00051946513092, 1.7e-06, 0.00080745068338,
0.00029358481485, 0.00048852120227, 0.0053075302458,
0.0025877438737, 0.30889220253, 1.2869403350e-05,
0.0184563624305147)) # values from previous repeated runs of maxp()
m <- skating_maxp # for ease of typing
m## babiakova butyrskaya cohen fontana giunchi gusmeroli hegel hubert
## 5.0472e-06 5.3045e-03 9.0241e-02 3.6833e-04 1.1373e-05 8.3169e-05 9.7856e-06 2.6898e-04
## hughes kettunen kopac kwan liashenko luca maniachenko meier
## 2.9832e-01 1.3764e-03 2.9967e-06 2.6551e-01 5.1886e-04 1.0000e-06 8.0688e-04 2.9281e-04
## onda robinson sebestyen slutskaya soldatova suguri volchkova
## 4.8802e-04 5.3083e-03 2.5883e-03 3.0889e-01 1.2693e-05 1.8457e-02 1.1384e-03equalp.test(skating)##
## Constrained support maximization
##
## data: skating
## null hypothesis: babiakova = butyrskaya = cohen = fontana = giunchi = gusmeroli = hegel = hubert = hughes = kettunen = kopac = kwan = liashenko = luca = maniachenko = meier = onda = robinson = sebestyen = slutskaya = soldatova = suguri = volchkova
## null estimate:
## babiakova butyrskaya cohen fontana giunchi gusmeroli hegel hubert
## 0.043478 0.043478 0.043478 0.043478 0.043478 0.043478 0.043478 0.043478
## hughes kettunen kopac kwan liashenko luca maniachenko meier
## 0.043478 0.043478 0.043478 0.043478 0.043478 0.043478 0.043478 0.043478
## onda robinson sebestyen slutskaya soldatova suguri volchkova
## 0.043478 0.043478 0.043478 0.043478 0.043478 0.043478 0.043478
## (argmax, constrained optimization)
## Support for null: -464.46 + K
##
## alternative hypothesis: sum p_i=1
## alternative estimate:
## babiakova butyrskaya cohen fontana giunchi gusmeroli hegel hubert
## 8.5593e-06 1.0413e-02 1.0465e-01 9.9774e-04 1.9198e-05 1.8303e-04 1.7491e-05 6.9266e-04
## hughes kettunen kopac kwan liashenko luca maniachenko meier
## 2.5423e-01 3.7718e-03 4.9234e-06 2.5893e-01 1.3404e-03 1.0000e-06 3.8385e-03 7.5582e-04
## onda robinson sebestyen slutskaya soldatova suguri volchkova
## 1.3209e-03 1.0957e-02 6.3675e-03 3.0368e-01 2.0746e-05 3.4651e-02 3.1466e-03
## (argmax, free optimization)
## Support for alternative: -233.35 + K
##
## degrees of freedom: 22
## support difference = 231.11
## p-value: 5.3208e-84dotchart(m)
dotchart(log(m))
1.1 Evidence for medal positions
Looking at the dotcharts, it seems that the medallists—Hughes (gold), Slutskya (silver), Kwan (bronze)—were considerably higher in strength than the rest of the field. Here I will test the hypothesis that the medallists were in fact the strongest three competitors. Technically you need to optimize over the union of the possibilities that one of the three medallists did not come in the top three; but this is hard. We will do something much easier but numerically equivalent: optimize over the union of outcomes where either Cohen or Suguri (who placed fourth and fifth respectively) had higher strength than any of the medallists.
f <- function(smaller,bigger){
jj <- rep(0,22)
names(jj) <- pnames(skating)[-23] # no fillup
jj[smaller] <- -1
jj[bigger ] <- +1
return(jj)
}
problem_constraints <- NULL
problem_constraints %<>% rbind(f("hughes" ,"cohen" ))
problem_constraints %<>% rbind(f("slutskaya","cohen" ))
problem_constraints %<>% rbind(f("kwan" ,"cohen" ))
problem_constraints %<>% rbind(f("hughes" ,"suguri"))
problem_constraints %<>% rbind(f("slutskaya","suguri"))
problem_constraints %<>% rbind(f("kwan" ,"suguri"))
rownames(problem_constraints) <-
c("hug < coh", "slu < coh", "kwa < coh", "hug < sug", "slu < sug", "kwa < sug" )
colnames(problem_constraints) %<>% substr(1,3)
problem_constraints## bab but coh fon giu gus heg hub hug ket kop kwa lia luc man mei ond rob seb slu sol sug
## hug < coh 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0
## slu < coh 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0
## kwa < coh 0 0 1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0
## hug < sug 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1
## slu < sug 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 1
## kwa < sug 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 1small <- 0.01
big <- 0.1
start_vector <- rep(small,22)
names(start_vector) <- pnames(skating)[-23] # no fillup
start_vector["cohen"] <- big
start_vector["suguri"] <- bigout <- rep(0,nrow(problem_constraints))
fullout <- list()
mle_skaters <- NULL
for(i in seq_len(nrow(problem_constraints))){
jj <- maxp(skating, startp=start_vector, give=TRUE,fcm=problem_constraints[i,], fcv=0,n=10,SMALL=1e-5)
fullout[[i]] <- jj
out[i] <- jj$value
mle_skaters <- cbind(mle_skaters,jj$par)
}
out## [1] -240.49 -240.87 -242.51 -245.34 -245.76 -246.36out-max(out)## [1] 0.00000 -0.37607 -2.02439 -4.84482 -5.27450 -5.86863colnames(mle_skaters) <- rownames(problem_constraints)
mle_skaters## hug < coh slu < coh kwa < coh hug < sug slu < sug kwa < sug
## babiakova 3.2667e-05 2.5014e-05 4.8079e-05 2.9086e-05 3.0673e-05 5.2544e-05
## butyrskaya 1.2137e-02 1.0329e-02 1.5217e-02 1.7104e-02 2.0783e-02 1.8437e-02
## cohen 1.8480e-01 1.9422e-01 1.8330e-01 1.4611e-01 1.2358e-01 1.3285e-01
## fontana 1.1025e-03 9.4827e-04 1.3901e-03 1.4003e-03 1.5262e-03 1.5581e-03
## giunchi 5.4870e-05 4.6039e-05 9.4139e-05 5.7172e-05 6.3842e-05 8.6316e-05
## gusmeroli 3.2183e-04 2.4757e-04 3.8966e-04 3.5637e-04 3.5552e-04 4.6524e-04
## hegel 5.6206e-05 4.3580e-05 8.7857e-05 5.2264e-05 5.8444e-05 8.2770e-05
## hubert 8.2030e-04 7.0930e-04 9.8510e-04 1.0355e-03 1.0860e-03 1.1334e-03
## hughes 1.8332e-01 3.2606e-01 1.7143e-01 9.5729e-02 3.2605e-01 2.7023e-01
## kettunen 3.6403e-03 3.2187e-03 4.7842e-03 4.9077e-03 5.3165e-03 5.2343e-03
## kopac 2.0396e-05 1.5873e-05 2.5018e-05 1.8586e-05 1.9063e-05 4.1954e-05
## kwan 2.2298e-01 2.1696e-01 1.8270e-01 2.4359e-01 2.8976e-01 8.9326e-02
## liashenko 1.4785e-03 1.5748e-03 1.4781e-03 1.9413e-03 2.1835e-03 2.2420e-03
## luca 1.0000e-05 1.0000e-05 1.0000e-05 1.0000e-05 1.0000e-05 1.0000e-05
## maniachenko 2.1493e-03 1.8427e-03 2.6640e-03 2.9463e-03 3.2808e-03 3.2110e-03
## meier 8.8519e-04 7.1397e-04 1.1080e-03 1.1292e-03 1.2017e-03 1.3137e-03
## onda 1.3968e-03 1.3830e-03 1.5750e-03 1.8229e-03 1.9951e-03 1.9881e-03
## robinson 1.2452e-02 1.1396e-02 1.4620e-02 1.7747e-02 1.9257e-02 1.9155e-02
## sebestyen 6.9877e-03 5.2787e-03 8.9093e-03 9.0378e-03 1.1345e-02 9.2573e-03
## slutskaya 3.3115e-01 1.9374e-01 3.6834e-01 3.5499e-01 9.3861e-02 3.4931e-01
## soldatova 5.9899e-05 5.5970e-05 9.6161e-05 6.5944e-05 6.9721e-05 9.5251e-05
## suguri 3.1103e-02 2.8599e-02 3.6279e-02 9.5833e-02 9.3879e-02 8.9330e-02Now compare these values with the unconstrained maximum likelihoods:
options(use_alabama = TRUE)
mgv <- maxp(skating, give=TRUE, n=1)$value
mgv## [1] -233.35mgv - out## [1] 7.1446 7.5207 9.1690 11.9894 12.4191 13.0132mgv - max(out)## [1] 7.1446that is, a little over 7 units of support, clearly exceeding the two
units suggested by Edwards. We have strong evidence to suggest that
the medallists were indeed the strongest three competitors. Observe
that the maximum likelihood among the six alternative hypotheses is
that of number 3, in which the maximization was constrained to obey
Kwan < Cohen. If we have to change the medallists, then
exchanging the order of Kwan and Cohen incurs the lowest likelihood
penalty. If, instead, we ask whether there is evidence that Suguri
[who was much weaker than Cohen] should not have been a medallist, we
find
mgv - max(out[4:6])## [1] 11.989Much stronger.
plot(out,ylab="log likelihood",axes=FALSE)
axis(2)
axis(side=1,at=1:6,srt=45,labels=c(
"Hug<Coh", "Slu<Coh", "Kwa<Coh",
"Hug<Sug", "Slu<Sug", "Kwa<Sug"
))
In the plot above, the vertical axis shows the support. The six
points on the x-axis correspond to the six rows of
problem_constraints; names have been abbreviated to the first three
letters. Thus the first three points are maximum likelihoods for
Hughes < Cohen, Slutskya < Cohen, and Kwan < Cohen respectively.
We may test the hypothesis that Slutskaya == Kwan == Hughes:
samep.test(skating,c("slutskaya","kwan","hughes"))##
## Constrained support maximization
##
## data: skating
## null hypothesis: slutskaya = kwan = hughes
## null estimate:
## babiakova butyrskaya cohen fontana giunchi gusmeroli hegel hubert
## 7.7128e-06 1.1710e-02 1.0659e-01 7.0844e-04 1.7800e-05 1.5591e-04 1.5842e-05 6.4838e-04
## hughes kettunen kopac kwan liashenko luca maniachenko meier
## 2.7228e-01 2.7147e-03 4.5348e-06 2.7228e-01 1.0225e-03 1.0108e-06 1.7857e-03 6.3340e-04
## onda robinson sebestyen slutskaya soldatova suguri volchkova
## 1.0171e-03 1.1618e-02 5.8138e-03 2.7228e-01 2.1518e-05 3.6249e-02 2.4317e-03
## (argmax, constrained optimization)
## Support for null: -231.2 + K
##
## alternative hypothesis: sum p_i=1
## alternative estimate:
## babiakova butyrskaya cohen fontana giunchi gusmeroli hegel hubert
## 8.5593e-06 1.0413e-02 1.0465e-01 9.9774e-04 1.9198e-05 1.8303e-04 1.7491e-05 6.9266e-04
## hughes kettunen kopac kwan liashenko luca maniachenko meier
## 2.5423e-01 3.7718e-03 4.9234e-06 2.5893e-01 1.3404e-03 1.0000e-06 3.8385e-03 7.5582e-04
## onda robinson sebestyen slutskaya soldatova suguri volchkova
## 1.3209e-03 1.0957e-02 6.3675e-03 3.0368e-01 2.0746e-05 3.4651e-02 3.1466e-03
## (argmax, free optimization)
## Support for alternative: -233.35 + K
##
## degrees of freedom: 2
## support difference = -2.1476
## p-value: 11.2 Ranking methodologies.
maximum likelihood estimated strengths furnish an ordering for the competitors. We can compare ranking by strengths with the official point-tallying method:
rL <- sort(skating_maxp,decreasing=TRUE)
rL[] <- seq_along(rL)
rO <- seq_len(nrow(skating_table))
names(rO) <- rownames(skating_table)
ordertransplot(rO,rL,xlab="offical rank",ylab="likelihood rank",main="Ladies free skating, 2002 Winter Olympics")
Figure 1.1: Official ranks vs likelihood rank, Ladies free skating, 2002 Winter Olympics
In figure 1.1, there is close agreement between the two methods in the large, but differences in detail. For example, the official ordering is hughes first, then slutskaya second, then kwan third; likelihood winner is slutskya, followed by hughes and then kwan. We can gain some understanding of this result by looking at the raw judging results for the three medallists.
2 Extraction methods for ordertables
The entire ordertable is given here again for convenience:
skating_table## An ordertable:
## J1 J2 J3 J4 J5 J6 J7 J8 J9
## hughes 1 4 3 4 1 2 1 1 1
## slutskaya 3 1 1 1 4 1 2 3 2
## kwan 2 3 2 2 2 3 3 2 3
## cohen 5 2 4 3 3 4 4 4 4
## suguri 4 8 5 5 5 7 5 5 5
## butyrskaya 6 5 8 7 12 5 8 7 6
## robinson 7 7 7 9 6 8 10 6 7
## sebestyen 8 10 12 8 7 6 12 8 8
## kettunen 9 9 13 6 13 10 7 11 14
## volchkova 10 6 14 11 10 12 6 9 15
## maniachenko 13 12 11 12 16 11 11 10 9
## fontana 14 11 18 16 9 15 9 12 10
## liashenko 15 13 6 10 8 14 13 14 16
## onda 11 14 10 15 15 13 15 13 11
## hubert 12 17 17 13 11 16 14 15 13
## meier 16 16 9 14 14 9 16 16 12
## gusmeroli 17 15 15 17 17 18 17 17 17
## soldatova 19 18 22 20 21 17 18 18 19
## hegel 20 21 16 22 18 19 21 19 18
## giunchi 18 19 20 21 19 20 20 20 20
## babiakova 22 20 19 19 20 21 19 22 22
## kopac 21 22 23 18 22 22 22 21 21
## luca 23 23 21 23 23 23 23 23 23Suppose we are interested only in hughes, slutskaya, and kwan;
we thus focus on the first three lines of skating_table. If we were
to consider the entries of the first three lines, see how hughes
came fourth according to J2. This cannot be converted to a
ranktable because the judges’ observations are not sequential. Each
column should be a permutation of \(\left\lbrace 1,2,3\right\rbrace\).
To do this, we use the extraction method:
skating_table[1:3,]## An ordertable:
## J1 J2 J3 J4 J5 J6 J7 J8 J9
## hughes 1 3 3 3 1 2 1 1 1
## slutskaya 3 1 1 1 3 1 2 3 2
## kwan 2 2 2 2 2 3 3 2 3Above we see just the order for these three competitors. This may be converted to a ranktable, or indeed a likelihood function:
as.ranktable(skating_table[1:3,])## A ranktable:
## c1 c2 c3
## J1 hughes kwan slutskaya
## J2 slutskaya kwan hughes
## J3 slutskaya kwan hughes
## J4 slutskaya kwan hughes
## J5 hughes kwan slutskaya
## J6 slutskaya hughes kwan
## J7 hughes slutskaya kwan
## J8 hughes kwan slutskaya
## J9 hughes slutskaya kwansuppfun(skating_table[1:3,])## log(hughes^6 * (hughes + kwan)^-4 * (hughes + kwan + slutskaya)^-9 * kwan^6 * (kwan +
## slutskaya)^-5 * slutskaya^6)and then calculating a table of the results:
hughes <- unclass(skating_table)[1,]
slutskaya <- unclass(skating_table)[2,]
table(hughes)## hughes
## 1 2 3 4
## 5 1 1 2table(slutskaya)## slutskaya
## 1 2 3 4
## 4 2 2 1Looking at the above, we see that Hughes has five judges who gave her
“1”s, while Slutskaya has only 4, which is why she was considered to
be the top. It seems that Hughes’s two “4”s (compared with
Slutskaya’s one) have cost her more likelihood-based strength than her
extra “1” gave her; note also that the two judges who placed Hughes
fourth (viz J2 and J4) placed Slutskaya first. Also:
table(hughes < slutskaya)##
## FALSE TRUE
## 4 5(note that equality is disallowed) showing that five judges preferred Hughes to Slutskaya and four preferred Slutsakaya to Hughes.
3 Just the top seven competitors
Now we consider only the top seven competitors (that is, the top seven according to Judge 1).
(o <- skating_table[1:7,])## An ordertable:
## J1 J2 J3 J4 J5 J6 J7 J8 J9
## hughes 1 4 3 4 1 2 1 1 1
## slutskaya 3 1 1 1 4 1 2 3 2
## kwan 2 3 2 2 2 3 3 2 3
## cohen 5 2 4 3 3 4 4 4 4
## suguri 4 7 5 5 5 6 5 5 5
## butyrskaya 6 5 7 6 7 5 6 7 6
## robinson 7 6 6 7 6 7 7 6 7(o <- as.ranktable(o))## A ranktable:
## c1 c2 c3 c4 c5 c6 c7
## J1 hughes kwan slutskaya suguri cohen butyrskaya robinson
## J2 slutskaya cohen kwan hughes butyrskaya robinson suguri
## J3 slutskaya kwan hughes cohen suguri robinson butyrskaya
## J4 slutskaya kwan cohen hughes suguri butyrskaya robinson
## J5 hughes kwan cohen slutskaya suguri robinson butyrskaya
## J6 slutskaya hughes kwan cohen butyrskaya suguri robinson
## J7 hughes slutskaya kwan cohen suguri butyrskaya robinson
## J8 hughes kwan slutskaya cohen suguri robinson butyrskaya
## J9 hughes slutskaya kwan cohen suguri butyrskaya robinsonNow create a likelihood function:
o <- suppfun(o)
o## log(butyrskaya^6 * (butyrskaya + cohen + hughes + kwan + robinson + slutskaya +
## suguri)^-9 * (butyrskaya + cohen + hughes + kwan + robinson + suguri)^-4 * (butyrskaya +
## cohen + hughes + robinson + suguri)^-2 * (butyrskaya + cohen + kwan + robinson +
## slutskaya + suguri)^-5 * (butyrskaya + cohen + kwan + robinson + suguri)^-3 * (butyrskaya
## + cohen + robinson)^-1 * (butyrskaya + cohen + robinson + slutskaya + suguri)^-3 *
## (butyrskaya + cohen + robinson + suguri)^-6 * (butyrskaya + hughes + kwan + robinson +
## suguri)^-1 * (butyrskaya + hughes + robinson + suguri)^-2 * (butyrskaya + robinson)^-7 *
## (butyrskaya + robinson + slutskaya + suguri)^-1 * (butyrskaya + robinson + suguri)^-8 *
## cohen^9 * hughes^9 * kwan^9 * robinson^4 * (robinson + suguri)^-2 * slutskaya^9 *
## suguri^8)and have some fun with it:
(o <- maxp(o))## butyrskaya cohen hughes kwan robinson slutskaya suguri
## 0.00124147 0.06015656 0.32275047 0.27376975 0.00076939 0.33661023 0.00470213pie(o)
Package dataset
Following lines create skating.rda, residing in the data/
directory of the package.
save(skating_table, skating, skating_maxp, file="skating.rda")References
hyper2 Package: Likelihood Functions for Generalized Bradley-Terry Models.” The R Journal 9 (2): 429–39.