In RobinHankin/hyper2: The Hyperdirichlet Distribution, Mark 2
Introduction
Formula 1 racing is an important and prestigious motor sport
\citep{codling2017,jenkins2010}. Season ranking is based on a points
allocation system wherein competitors are awarded points based on race
finishing order; points accumulate additively. The overall
competition winner is the competitor who accumulates the most points
after the final race.
However, in the case of Formula 1 Motor Racing, the points system is
the subject of much controversy, having changed often since the
competition's inauguration in 1950 when the points allocation, used
from 1950 to 1959 was $(8,6,4,3,2)$, thus crediting only the first
five finishers. As of 2020, the current points system of
$(25,18,15,12,10,8,6,4,2,1)$ credits the first 10. Arguably these two
systems could introduce different rational behaviour under zero-sum
assumptions: if, in a race, a driver knows he will place seventh under
a low-risk strategy but may place sixth by dint of driving more
aggressively, the low-risk strategy might be rational under the first
points system (which does not reward the extra ranking), but not under
the second, which does. Still, it is reasonable to assume that each
driver strives to maximise his rank and this will be done here. If
this is so, then changing the points system might change the
competitors' rankings but not their behaviour: surely a defect of
using points to rank the competitors.
Given that racing is a zero-sum game---and that points are
monotonically decreasing---each player will try to get as high a rank
as possible regardless of the actual points system used. However,
there are other consistent interpretations. \cite{bakhrankova2011},
for example, considers the possibility of inter-driver collusion, a
phenomenon not pursued here; and authors such as
\citeauthor{mastromarco2009} suggest that the frequency of rule
changes is driven by factors such as driver safety and revenue
optimization.
Points systems similar to that of Formula 1 are common in other sports
and examples would include modern pentathlon, curling, chess, rowing,
and figure skating; all have the common feature of translating ranks
into scores which combine additively to generate an overall ranking.
Further, we see points systems frequently in the wider context of
competitive situations such as the Eurovision Song Contest, Australian
MasterChef, University rankings, educational league tables, car safety
ratings, and many others.
Bradley-Terry and generalizations for rank statistics
The Bradley-Terry model \citep{bradley1952} assigns non-negative
strengths $p_1,\ldots, p_n$ to each of $n$ competitors in such a way
that the probability of $i$ beating $j\neq i$ in pairwise competition
is $\frac{p_i}{p_i+p_j}$; it is conventional to normalize so that
$\sum p_i=1$. Further, we use a generalization due to
\citet{luce1959}, in which the probability of competitor $i$ winning
in a field of $1,\ldots, n$ is $\frac{p_i}{p_1+\cdots +p_n}$. Noting
that there is information in the whole of the finishing order, and not
just the first across the line, we can follow \cite{plackett1975} and
consider the runner-up to be the winner among the remaining
competitors, and so on down the finishing order. Without loss of
generality, if the order of finishing were $1,2,3,4,5$, then a
suitable \citeauthor{plackett1975}-\citeauthor{luce1959} likelihood
function would be
\begin{equation}\label{competitors_1_to_5_likelihood}
\frac{p_1}{p_1+p_2+p_3+p_4+p_5}\cdot
\frac{p_2}{p_2+p_3+p_4+p_5}\cdot
\frac{p_3}{p_3+p_4+p_5}\cdot
\frac{p_4}{p_4+p_5}\cdot
\frac{p_5}{p_5}
\end{equation}
and this would be a forward ranking Plackett-Luce model in the
terminology of \cite{johnson2020}. A slight generalization allows the
incorporation of non-finishers (DNF etc). If, say, competitors 4 and
5 did not finish, we would have
\begin{equation}\label{competitors_1_to_3_only_finished}
\frac{p_1}{p_1+p_2+p_3+p_4+p_5}\cdot
\frac{p_2}{p_2+p_3+p_4+p_5}\cdot
\frac{p_3}{p_3+p_4+p_5}
\end{equation}
(observe how this likelihood function, while informative about
$p_4+p_5$, is uninformative about $p_4\left|p_4+p_5\right.$). We now
use a technique due to \cite{hankin2010,hankin2017} and introduce
fictional (reified) entities whose nonzero Bradley-Terry strength
helps certain competitors or sets of competitors under certain
conditions. The canonical example would be the home-ground advantage
in association football. If players (teams) $1,2$ with strengths
$p_1,p_2$ compete, and if our observation were $a$ home wins and $b$
away wins for team $1$, and $c$ home wins and $d$ away wins for team
$2$, then a suitable likelihood function would be
[
\left(\frac{p_1+p_H}{p_1+p_2+p_H}\right)^a
\left(\frac{p_1}{p_1+p_2+p_H}\right)^b
\left(\frac{p_2+p_H}{p_1+p_2+p_H}\right)^c
\left(\frac{p_2+p_H}{p_1+p_2+p_H}\right)^d
]
where $p_H$ is a quantification of the beneficial home ground effect.
Similar techniques have been used to account for the first-move
advantage in chess, and effective coordination between members of
doubles tennis teams; we may use a similar device to account for pole
position in Formula 1. Here I analyse seasons 2016-2019 using the
hyper2 package \citep{hankin2017} which implements the Plackett-Luce
likelihood function with additional reified entities.
Formula 1 dataset
Taking 2017 as an example, we see the rank table as follows:
\small
\begin{CodeChunk}
\begin{CodeOutput}
AUS CHN BHR RUS ESP MON CAN AZE AUT GBR HUN BEL ITA SIN MAL JPN USA MEX BRA ABU
Hamilton 2 1 2 4 1 7 1 5 4 1 4 1 1 1 2 1 1 9 4 2
Vettel 1 2 1 2 2 1 4 4 2 7 1 2 3 Ret 4 Ret 2 4 1 3
Bottas 3 6 3 1 Ret 4 2 2 1 2 3 5 2 3 5 4 5 2 2 1
Raikkonen 4 5 4 3 Ret 2 7 14 5 3 2 4 5 Ret DNS 5 3 3 3 4
Ricciardo Ret 4 5 Ret 3 3 3 1 3 5 Ret 3 4 2 3 3 Ret Ret 6 Ret
...
Hartley 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 13 Ret Ret 15
Button 0 0 0 0 0 Ret 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Resta 0 0 0 0 0 0 0 0 0 0 Ret 0 0 0 0 0 0 0 0 0
\end{CodeOutput}
\end{CodeChunk}
top <- 11
Each row is a driver and each column (after the first) a venue. We
see that Hamilton, the first row, came second in Australia, first in
China, second in Bahrain, fourth in Russia, and so on (Hartley,
Button, and Resta placed last). In the first column, we see the
result from Australia (AUS) in which Hamilton came second, Vettel
first, Bottas third, and so on. Notation used also includes six
different classes of no-score such as Ret for retired, WD for
withdrawn, and so on; a zero entry means "did not finish". It is
straightforward to translate this table into a Plackett-Luce
log-likelihood function using the hyper2 package; for simplicity we
will consider only the r top top-ranked drivers\footnote{In the
Plackett-Luce likelihood function, the performance of lower-ranked
players can be weakly informative about higher-ranked players'
strengths. For example, we see that Vettel retired twice (in
Singapore and Japan), so any player who placed in those venues will
effectively steal'' strength from Vettel, and generallygive'' it
to Hamilton or Bottas.}. Although it has many terms, the overall loglikelihood
expression is of the general form
library("hyper2",quietly=TRUE)
library("magrittr",quietly=TRUE)
top <- 11
H2016 <- ordertable2supp(F1_table_2016[seq_len(top),])
H2017 <- ordertable2supp(F1_table_2017[seq_len(top),])
H2018 <- ordertable2supp(F1_table_2018[seq_len(top),])
H2019 <- ordertable2supp(F1_table_2019[seq_len(top),])
\newcommand{\pham}{p_\mathrm{Ham}}
\newcommand{\pvet}{p_\mathrm{Vet}}
\newcommand{\pbot}{p_\mathrm{Bot}}
\newcommand{\prai}{p_\mathrm{Rai}}
\newcommand{\pric}{p_\mathrm{Ric}}
\newcommand{\pver}{p_\mathrm{Ver}}
\newcommand{\pper}{p_\mathrm{Per}}
\newcommand{\poco}{p_\mathrm{Oco}}
\newcommand{\psai}{p_\mathrm{Sai}}
\newcommand{\phul}{p_\mathrm{Hul}}
\newcommand{\pmas}{p_\mathrm{Mas}}
\begin{equation}
\frac{\pham^{20}\,\pmas^{16}\,\pbot^{19}\ldots
}{\parbox{4in}{$(\pper+\poco)(\pver+\pper+\psai+\phul+\pmas)\
\rule{10mm}{0mm}(\pric+\pper+\pmas)(\psai+\pmas)\ \rule{20mm}{0mm}
(\pver + \pper+ \poco + \psai+ \phul)\ldots
$}}
\end{equation}
Finding the maximum likelihood estimate for the players' strengths is
straightforward. The hyper2 package includes a suite of numerical
optimization routines, and because they have access to derivatives,
convergence is rapid. A graphical diagram of the strengths is given
in figure \ref{piechartstrength}
mH2016 <- maxp(H2016)
mH2017 <- maxp(H2017)
mH2018 <- maxp(H2018)
mH2019 <- maxp(H2019)
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, seasons 2016-2019"}
par(mfrow=c(2,2))
par(mar = c(2,2,2,2)/2)
pie(mH2016[c(1,6,11,2,7,3,8,4,9,5,10)],cex=0.7,main="2016")
pie(mH2017[c(1,6,11,2,7,3,8,4,9,5,10)],cex=0.7,main="2017")
pie(mH2018[c(1,6,11,2,7,3,8,4,9,5,10)],cex=0.7,main="2018")
pie(mH2019[c(1,6,11,2,7,3,8,4,9,5,10)],cex=0.7,main="2019")
```r
library("qcc")
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, 2016"}
pareto.chart(mH2016[c(1,6,11,2,7,3,8,4,9,5,10)],main="2016",ylab="strength")
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, 2017"}
pareto.chart(mH2017[c(1,6,11,2,7,3,8,4,9,5,10)],main="2017",ylab="strength")
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, 2018"}
pareto.chart(mH2018[c(1,6,11,2,7,3,8,4,9,5,10)],main="2018",ylab="strength")
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, 2019"}
pareto.chart(mH2019[c(1,6,11,2,7,3,8,4,9,5,10)],main="2019",ylab="strength")
We see in figure \ref{piechartstrength} that the driver with the
largest estimated strength was Rosberg in 2016 at about 30\%, and in
years 2017-8-9 was Hamilton (at about 29\%, 37\%, and 42\%
respectively). As an illustration of the value of likelihood methods
(as opposed to points-based methods), a likelihood ratio test
[samep.test(), supplied with the package] rejects the null that
Hamilton and Vettel have the same strength in 2018 ($H_0\colon \pham=
\pvet$) with a likelihood ratio of $e^{2.76}\simeq 15.8$,
corresponding (by Wilks's theorem) to an asymptotic $p$-value of about
$0.02$. We may also use the reified entity concept to test the null
hypothesis that Hamilton's strength remains constant 2016, where
Rosberg had the highest estimated strength, and 2017-9, where Hamilton
did; we fail to reject this null.
H2016a <- H2016 %>% psubs("Hamilton","Hamilton_pre")
H2017a <- H2017 %>% psubs("Hamilton","Hamilton_post")
H2018a <- H2018 %>% psubs("Hamilton","Hamilton_post")
H2019a <- H2019 %>% psubs("Hamilton","Hamilton_post")
Ha <- H2016a + H2017a + H2018a + H2019a
Do a statistical test:
ignorethis <- samep.test(Ha,c("Hamilton_pre","Hamilton_post"))
ignorethis
Do another statistical test:
ignore <- samep.test(H2018,c("Hamilton","Vettel"))
Likelihood scoring vs points scoring
Applying the current points system, for example, to the 2017 results table we
would rank the drivers as follows:
real <- ordertable2points(ranktable_to_ordertable(wikitable_to_ranktable(F1_table_2017)),formula1_points_systems(top)$real)
real <- round(real[seq_len(top)],2)
real <- sort(real,decreasing=TRUE)
real[] <- seq_along(real)
\begin{multline}
\mbox{Hamilton}\succ
\mbox{Vettel}\succ
\mbox{Bottas}\succ
\mbox{Raikkonen}\succ\mbox{YES I KNOW}
\mbox{Ricciardo}\succ
\mbox{Verstappen}\succ\
\mbox{Perez1}\succ
\mbox{Ocon}\succ
\mbox{Sainz}\succ
\mbox{Hulkenberg}\mbox{YES I KNOW}\succ
\mbox{Massa}
\end{multline}
(this happens to be identical to the ranking after the extra point was
awarded for fastest lap). However, if we were award points according
to Zipf's law we would have:
zipf <- ordertable2points(ranktable_to_ordertable(wikitable_to_ranktable(F1_table_2017)),formula1_points_systems(top)$zipf)
zipf <- round(zipf[seq_len(top)],2)
zipf <- sort(zipf,decreasing=TRUE)
zipf[] <- seq_along(zipf)
\begin{multline}
\mbox{Hamilton}\succ
\mbox{Vettel}\succ
\mbox{Bottas}\succ
\mbox{Ricciardo}\succ
\mbox{Verstappen}\succ
\mbox{Raikkonen}\succ\mbox{YES I KNOW}\
\mbox{Perez}\succ
\mbox{Ocon}\succ
\mbox{Massa}\succ
\mbox{Sainz}\succ
\mbox{Hulkenberg}
\end{multline}
Thus these two systems agree on the first three places but fourth is
awarded to Raikkonen under the F1 system and Ricciardo under Zipf.
Compare the likelihood ranking:
jj <- mH2017
jj <- sort(jj,decreasing=TRUE)
jj[] <- seq_along(jj)
jj
\begin{multline}
\mbox{Hamilton}\succ
\mbox{Vettel}\succ
\mbox{Bottas}\succ
\mbox{Raikkonen}\mbox{YES I KNOW}\succ
\mbox{Ocon}\succ
\mbox{Ricciardo}\succ\
\mbox{Verstappen}\succ
\mbox{Perez}\succ
\mbox{Massa}\succ
\mbox{Sainz}\succ
\mbox{Hulkenberg}\mbox{YES I KNOW}
\end{multline}
Thus, for 2017 at least, we see that the current points system agrees
with likelihood ranking for the top four places, while Zipf's law
agrees to three.
p <- formula1_points_systems()[[1]] # real points system
dp2016 <- ordertable2points(as.ordertable(F1_table_2016[seq_len(top),]),p)
dp2016 <- sort(dp2016,decreasing=TRUE)
dp2016[] <- seq_along(dp2016)
m2016 <- maxp(H2016)
m2016 <- sort(m2016,decreasing=TRUE)
m2016[] <- seq_along(m2016)
dp2017 <- ordertable2points(as.ordertable(F1_table_2017[seq_len(top),]),p)
dp2017 <- sort(dp2017,decreasing=TRUE)
dp2017[] <- seq_along(dp2017)
m2017 <- maxp(H2017)
m2017 <- sort(m2017,decreasing=TRUE)
m2017[] <- seq_along(m2017)
dp2018 <- ordertable2points(as.ordertable(F1_table_2018[seq_len(top),]),p)
dp2018 <- sort(dp2018,decreasing=TRUE)
dp2018[] <- seq_along(dp2018)
m2018 <- maxp(H2018)
m2018 <- sort(m2018,decreasing=TRUE)
m2018[] <- seq_along(m2018)
dp2019 <- ordertable2points(as.ordertable(F1_table_2019[seq_len(top),]),p)
dp2019 <- sort(dp2019,decreasing=TRUE)
dp2019[] <- seq_along(dp2019)
m2019 <- maxp(H2019)
m2019 <- sort(m2019,decreasing=TRUE)
m2019[] <- seq_along(m2019)
```r top-ranked 11 drivers in Formula 1 seasons 2016-2019. Horizontal axis gives official (points-based) order, and the vertical axis gives the likelihood horder. Thus, taking 2017 as an example, the points-based ordering would be Hamilton first, then Vettel, then R\"{a}ikk\"{o}nen; YES I KNOW IT IS WRONG: RMARKDOWN IS DEFECTIVE while the likelihood ordering is (reading vertically) Hamilton, Vettel, Bottas ALSO R PLOTTING SUCKS HERE AND I WILL HAVE TO CHANGE THE PDF BY HAND WITH INKSCAPE TO MAKE IT LOOK GOOD"}
par(mfrow=c(2,2))
par(mar = c(3,3,3,3))
par(pty='s') # square plot
ox <- dp2016
print("m2016")
dput(m2016)
print("dp2016")
dput(dp2016)
oy <- ordertrans(m2016,names(dp2016))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2016 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2017
oy <- ordertrans(m2017,names(dp2017))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2017 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2018
oy <- ordertrans(m2018,names(dp2018))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2018 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2019
oy <- ordertrans(m2019,names(dp2019))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2019 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
```r
par(pty="s")
ox <- dp2016
oy <- ordertrans(m2016,names(dp2016))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2016 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty="s")
ox <- dp2017
oy <- ordertrans(m2017,names(dp2017))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2017 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty="s")
ox <- dp2018
oy <- ordertrans(m2018,names(dp2018))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2018 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty="s")
ox <- dp2019
oy <- ordertrans(m2019,names(dp2019))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2019 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
We can plot one ranking against the other\footnote{The historically
correct points awarded to the drivers differs from that calculated
here. That is for two reasons: firstly, ``fastest lap'' points are
not included here, and also the truncation of the order table to the
first 11 drivers can increase the rank of a driver if non-first-11
drivers are placed}, shown in Figure
\ref{orderingbypointsandlikelihood}. We thus see a comparison between
the two ordering systems. Taking 2017 as an example, drivers Hamilton
and Vettel are are respectively first and second according to both
ranking ranking procedures; but R\"{a}ikk\"{o}nen YES I KNOW IT IS
WRONG: RMARKDOWN IS DEFECTIVE and Bottas are third and fourth, and
fourth and third, according to points and likelihood respectively. We
define the \emph{degree of agreement} between the two ranking systems
as the maximal value of $r$ such that places $1,2,\ldots r$ all match.
Thus, from figure \ref{orderingbypointsandlikelihood}, the degree of
agreement between likelihood ranking and points ranking for the years
2016-2019 would be 8,2,6, and 2 respectively.
However, the points system used is essentially arbitrary. We could
use, for example, Zipf's law to rank the drivers: award one point to
the winner, half a point to second place, one third of a point to
third, and so on; see figure \ref{orderingbypointsandlikelihood_zipf}
in which we see generally poorer agreement between points-based ranks
and likelihood-based ranks, with a degree of agreement of 4,2,2,2 for
the years 2016-2019 respectively. This might be an indication that
using a Zipf law points allocation is objectively worse than the
current points system. There are a number of plausible points systems
that might be used:
pz <- formula1_points_systems()$zipf
dp2016z <- ordertable2points(as.ordertable(F1_table_2016[seq_len(top),]),pz)
dp2016z <- sort(dp2016z,decreasing=TRUE)
dp2016z[] <- seq_along(dp2016z)
dp2017z <- ordertable2points(as.ordertable(F1_table_2017[seq_len(top),]),pz)
dp2017z <- sort(dp2017z,decreasing=TRUE)
dp2017z[] <- seq_along(dp2017z)
dp2018z <- ordertable2points(as.ordertable(F1_table_2018[seq_len(top),]),pz)
dp2018z <- sort(dp2018z,decreasing=TRUE)
dp2018z[] <- seq_along(dp2018z)
dp2019z <- ordertable2points(as.ordertable(F1_table_2019[seq_len(top),]),pz)
dp2019z <- sort(dp2019z,decreasing=TRUE)
dp2019z[] <- seq_along(dp2019z)
```r \label{orderingbypointsandlikelihood_zipf} but points ranking calculated by Zipf's law. REF DOES NOT WORK HERE Note the generally poorer agreement between the two ranking systems ALSO R PLOTTING SUCKS HERE AND I WILL HAVE TO CHANGE THE PDF BY HAND WITH INKSCAPE TO MAKE IT LOOK GOOD"}
par(mfrow=c(2,2))
par(mar = c(3,3,3,3))
par(pty='s') # square plot
ox <- dp2016z
oy <- ordertrans(m2016,names(dp2016z))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2016 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2017z
oy <- ordertrans(m2017,names(dp2017z))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2017 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2018z
oy <- ordertrans(m2018,names(dp2018z))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2018 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2019z
oy <- ordertrans(m2019,names(dp2019z))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order",ylab="likelihood order",main='2019 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
```r
par(pty='s') # square plot
ox <- dp2016z
oy <- ordertrans(m2016,names(dp2016z))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order (Zipf)",ylab="likelihood order",main='2016 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2017z
oy <- ordertrans(m2017,names(dp2017z))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order (Zipf)",ylab="likelihood order",main='2017 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2018z
oy <- ordertrans(m2018,names(dp2018z))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order (Zipf)",ylab="likelihood order",main='2018 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
par(pty='s') # square plot
ox <- dp2019z
oy <- ordertrans(m2019,names(dp2019z))
plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16,
xlab="points order (Zipf)",ylab="likelihood order",main='2019 season')
par(xpd=TRUE) # allow drivers' names to appear outside plotting region
for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) }
par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region
abline(0,1)
\begin{itemize}
\item The current Formula 1 system $(25,18,15,12,10,8,6,4,2,1)$
\item The inaugural Formula 1 system $(8,6,4,3,2)$
\item Zipf's law $(1,\frac{1}{2},\frac{1}{3},\frac{1}{4},\ldots)$
\item Linear decrease $(n,n-1,n-2,\ldots,3,2,1,0)$
\item Halving system: $(1,\frac{1}{2},\frac{1}{4},\frac{1}{8},\ldots)$
\item A ``winner takes all'' system $(1,0,\ldots)$
\end{itemize}
We note that some of these may be generalized; we might consider a
more general linear points system $(r,r-1,r-2,\ldots,3,2,1,0)$ for
fixed integer $0<r<n$; the halving system can be generalized to a
geometric distribution; and the winner takes all system can be
replaced by giving equal points to the top $r$ competitors
$(\underbrace{1,1,\ldots ,1}_{r},0,0,\ldots)$, for some integer $r<n$.
agreement2016 <- rep(NA,length(formula1_points_systems(top)))
agreement2017 <- rep(NA,length(formula1_points_systems(top)))
agreement2018 <- rep(NA,length(formula1_points_systems(top)))
agreement2019 <- rep(NA,length(formula1_points_systems(top)))
for(i in seq_along(formula1_points_systems(top))){
pa <- formula1_points_systems(top)[[i]]
dp2016a <- ordertable2points(as.ordertable(F1_table_2016[seq_len(top),]),pa)
dp2016a <- sort(dp2016a,decreasing=TRUE)
dp2016a[] <- seq_along(dp2016a)
agreement2016[i] <- sum(cumprod(names(dp2016a) == names(m2016)))
dp2017a <- ordertable2points(as.ordertable(F1_table_2017[seq_len(top),]),pa)
dp2017a <- sort(dp2017a,decreasing=TRUE)
dp2017a[] <- seq_along(dp2017a)
agreement2017[i] <- sum(cumprod(names(dp2017a) == names(m2017)))
dp2018a <- ordertable2points(as.ordertable(F1_table_2018[seq_len(top),]),pa)
dp2018a <- sort(dp2018a,decreasing=TRUE)
dp2018a[] <- seq_along(dp2018a)
agreement2018[i] <- sum(cumprod(names(dp2018a) == names(m2018)))
dp2019a <- ordertable2points(as.ordertable(F1_table_2019[seq_len(top),]),pa)
dp2019a <- sort(dp2019a,decreasing=TRUE)
dp2019a[] <- seq_along(dp2019a)
agreement2019[i] <- sum(cumprod(names(dp2019a) == names(m2019)))
} # i loop closes
agreement4yearsallp <- cbind(
agreement2016,
agreement2017,
agreement2018,
agreement2019
)
rownames(agreement4yearsallp) <- names(formula1_points_systems(top))
colnames(agreement4yearsallp) <- 2016:2019
wanted <- c(1,2,24,23,3,29)
It is straightforward to calculate the degree of agreement for the
observed rank table for years 2016-2019, shown in the following table:
agreement4yearsallp[wanted,]
It is clear that there is no points system that is the best for all
four years. However, observing that both the likelihood ranking and
the points ranking are random variables in this paradigm suggests a
method whereby we can objectively assess a given points system. Using
sampling techniques we can repeatedly generate an order table \emph{in
silico}, using estimated driver strengths from the observed tables.
For each of, say, 1000 such synthetic tables, calculate drivers'
maximum likelihood \citeauthor{plackett1975} strengths, and also their
points awarded according to any given points system. We then compare
rankings generated by the \citeauthor{plackett1975} strengths and the
points awarded and note the degree of agreement between the two, as
measured by the number of rankings correctly predicted. This
furnishes an objective assessment of the points system used.
Numerical results
load("formula1_results_2016.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist
rownames(OO)[3] <- "wta"
OO2016 <- OO
load("formula1_results_2017.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist
rownames(OO)[3] <- "wta"
OO2017 <- OO
load("formula1_results_2018.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist
rownames(OO)[3] <- "wta"
OO2018 <- OO
load("formula1_results_2019.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist
rownames(OO)[3] <- "wta"
OO2019 <- OO
summarytablemaker <- function(OO){
out <- data.frame(
means = apply(OO,1,mean),
winner_correct = rowSums(OO>1)/ncol(OO),
all_correct = rowSums(OO==top)/ncol(OO)
)
rownames(out) <- names(pointslist)
return(out)
}
s2016 <- summarytablemaker(OO2016)
s2017 <- summarytablemaker(OO2017)
s2018 <- summarytablemaker(OO2018)
s2019 <- summarytablemaker(OO2019)
```r systems when compared with a Placket-Luce strength ordering"}
jj <- cbind(s2016[,1],s2017[,1],s2018[,1],s2019[,1]) # mean
typecol <- c(rep("red",2),rep("orange",10),rep("yellow",11),rep("green",1),rep("blue",7),rep("purple",7))
matplot(jj,type="n",xlab='points system',ylab='mean correct places')
points(s2016[,1],col=typecol,pch=1)
points(s2017[,1],col=typecol,pch=2)
points(s2018[,1],col=typecol,pch=3)
points(s2019[,1],col=typecol,pch=4)
legend("topright",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
matplot(jj[wanted,],type="n",xlab='points system',ylab='mean correct places',axes=FALSE)
jjlabels <- names(formula1_points_systems(top)[wanted])
jjlabels[5] <- "wta"
axis(1,at=1:6,labels=jjlabels)
axis(2,at=1:4)
points(s2016[wanted,1],pch=1,type='b',col='red')
points(s2017[wanted,1],pch=2,type='b',col='orange')
points(s2018[wanted,1],pch=3,type='b',col='green')
points(s2019[wanted,1],pch=4,type='b',col='blue')
legend("topright",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
par(xpd=TRUE)
text(x=6,y=s2016[wanted[6],1],pos=4,"2016",col='red')
text(x=6,y=s2017[wanted[6],1],pos=4,"2017",col='orange')
text(x=6,y=s2018[wanted[6],1],pos=4,"2018",col='green')
text(x=6,y=s2019[wanted[6],1],pos=4,"2019",col='blue')
```r systems"}
jj <- cbind(s2016[,2],s2017[,2],s2018[,2],s2019[,2]) # prob of winner
#matplot(jj,type="n",xlab='points system',ylab='probability of winner correct')
#points(s2016[,2],col=typecol,pch=1)
#points(s2017[,2],col=typecol,pch=2)
#points(s2018[,2],col=typecol,pch=3)
#points(s2019[,2],col=typecol,pch=4)
#legend("topright",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
matplot(jj[wanted,],type="n",xlab='points system',ylab='probability of winner correct',axes=FALSE)
axis(1,at=1:6,labels=jjlabels)
axis(2)
points(s2016[wanted,2],pch=1,type='b',col='red')
points(s2017[wanted,2],pch=2,type='b',col='orange')
points(s2018[wanted,2],pch=3,type='b',col='green')
points(s2019[wanted,2],pch=4,type='b',col='blue')
#legend(x=3.4,y=0.67,lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
par(xpd=TRUE)
text(x=6,y=s2016[wanted[6],2],pos=4,"2016",col='red')
text(x=6,y=s2017[wanted[6],2],pos=4,"2017",col='orange')
text(x=6,y=s2018[wanted[6],2],pos=4,"2018",col='green')
text(x=6,y=s2019[wanted[6],2],pos=4,"2019",col='blue')
```r systems"}
jj <- cbind(s2016[,3],s2017[,3],s2018[,3],s2019[,3]) # total rank order correct
matplot(jj,type="n",xlab='points system',ylab='probability of total rank order correct')
points(s2016[,3],col=typecol,pch=1)
points(s2017[,3],col=typecol,pch=2)
points(s2018[,3],col=typecol,pch=3)
points(s2019[,3],col=typecol,pch=4)
legend("topright",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
matplot(jj[wanted,],type="n",xlab='points system',ylab='probability of total rank order correct',axes=FALSE)
axis(1,at=1:6,labels=jjlabels)
axis(2)
points(s2016[wanted,3],pch=1,type='b',col='red')
points(s2017[wanted,3],pch=2,type='b',col='orange')
points(s2018[wanted,3],pch=3,type='b',col='green')
points(s2019[wanted,3],pch=4,type='b',col='blue')
legend("topleft",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
par(xpd=TRUE)
text(x=5,y=s2016[wanted[5],3],pos=3,"2016",col='red')
text(x=5,y=s2017[wanted[5],3],pos=1,"2017",col='orange')
text(x=5,y=s2018[wanted[5],3],pos=3,"2018",col='green')
text(x=5,y=s2019[wanted[5],3],pos=3,"2019",col='blue')
```
We now assess the six points systems using the methodology outlined
above, using 1000 in silico trials. For each of the six points
systems, each of the 1000 trials results in a single non-negative
integer: the degree of agreement between the Plackett-Luce ranking and
the rankings according to the points system considered. There are
three measures that might be used to assess the distribution of degree
of agreement: (1), the mean degree of agreement $\overline{d}$; (2),
the probability of correctly predicting the winner,
$\operatorname{Pr}(d\geqslant 1)$; and (3), the probability of
correctly predicting the complete order statistic
$\operatorname{Pr}(D=n)$. These summaries are shown for each of the
six points systems in figures \ref{simulated_mean},
\ref{simulated_winner_correct} and \ref{simulated_complete_agreement}
respectively.
We see that the linear points system $(n,n-1,n-2,\ldots,3,2,1)$ gives
the highest mean number of places, and also the highest probability of
correctly predicting the winner. However, the probability of
predicting the complete order statistic is maximized using the winnner
takes all system $(1,0,0,\ldots)$.
Discussion
Many competitive situations involve a ranking of the participants and
one common method for combining repeated observations is to assign
points for the winner, second place, third place, etc. The overall
ranking is decided on the basis of accumulated points. However,
points systems are unavoidably arbitrary. Assuming that participants
attempt to optimize their placing means that the arbitrary
nature of any points system.
Maximum likelihood estimation of participants' Plackett-Luce strengths
furnishes a ranking system that does not suffer from the arbitrariness
of a points system.
The points system used in Formula 1 motor racing is a source of lively
debate from many perspectives, with changes being controversial. Here
I treat total points scored as a random variable and compare different
point allocation schemes against objective Plackett-Luce ranks. Of
the six points systems considered here, a linear system or a
winner-takes-all-system appear to be closest to objective
Plackett-Luce, depending on the exact definition of ``closest''.
\bibliography{chess}
RobinHankin/hyper2 documentation built on April 21, 2024, 11:38 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Introduction
Formula 1 racing is an important and prestigious motor sport \citep{codling2017,jenkins2010}. Season ranking is based on a points allocation system wherein competitors are awarded points based on race finishing order; points accumulate additively. The overall competition winner is the competitor who accumulates the most points after the final race.
However, in the case of Formula 1 Motor Racing, the points system is the subject of much controversy, having changed often since the competition's inauguration in 1950 when the points allocation, used from 1950 to 1959 was $(8,6,4,3,2)$, thus crediting only the first five finishers. As of 2020, the current points system of $(25,18,15,12,10,8,6,4,2,1)$ credits the first 10. Arguably these two systems could introduce different rational behaviour under zero-sum assumptions: if, in a race, a driver knows he will place seventh under a low-risk strategy but may place sixth by dint of driving more aggressively, the low-risk strategy might be rational under the first points system (which does not reward the extra ranking), but not under the second, which does. Still, it is reasonable to assume that each driver strives to maximise his rank and this will be done here. If this is so, then changing the points system might change the competitors' rankings but not their behaviour: surely a defect of using points to rank the competitors.
Given that racing is a zero-sum game---and that points are monotonically decreasing---each player will try to get as high a rank as possible regardless of the actual points system used. However, there are other consistent interpretations. \cite{bakhrankova2011}, for example, considers the possibility of inter-driver collusion, a phenomenon not pursued here; and authors such as \citeauthor{mastromarco2009} suggest that the frequency of rule changes is driven by factors such as driver safety and revenue optimization.
Points systems similar to that of Formula 1 are common in other sports and examples would include modern pentathlon, curling, chess, rowing, and figure skating; all have the common feature of translating ranks into scores which combine additively to generate an overall ranking. Further, we see points systems frequently in the wider context of competitive situations such as the Eurovision Song Contest, Australian MasterChef, University rankings, educational league tables, car safety ratings, and many others.
Bradley-Terry and generalizations for rank statistics
The Bradley-Terry model \citep{bradley1952} assigns non-negative strengths $p_1,\ldots, p_n$ to each of $n$ competitors in such a way that the probability of $i$ beating $j\neq i$ in pairwise competition is $\frac{p_i}{p_i+p_j}$; it is conventional to normalize so that $\sum p_i=1$. Further, we use a generalization due to \citet{luce1959}, in which the probability of competitor $i$ winning in a field of $1,\ldots, n$ is $\frac{p_i}{p_1+\cdots +p_n}$. Noting that there is information in the whole of the finishing order, and not just the first across the line, we can follow \cite{plackett1975} and consider the runner-up to be the winner among the remaining competitors, and so on down the finishing order. Without loss of generality, if the order of finishing were $1,2,3,4,5$, then a suitable \citeauthor{plackett1975}-\citeauthor{luce1959} likelihood function would be
\begin{equation}\label{competitors_1_to_5_likelihood} \frac{p_1}{p_1+p_2+p_3+p_4+p_5}\cdot \frac{p_2}{p_2+p_3+p_4+p_5}\cdot \frac{p_3}{p_3+p_4+p_5}\cdot \frac{p_4}{p_4+p_5}\cdot \frac{p_5}{p_5} \end{equation}
and this would be a forward ranking Plackett-Luce model in the terminology of \cite{johnson2020}. A slight generalization allows the incorporation of non-finishers (DNF etc). If, say, competitors 4 and 5 did not finish, we would have
\begin{equation}\label{competitors_1_to_3_only_finished} \frac{p_1}{p_1+p_2+p_3+p_4+p_5}\cdot \frac{p_2}{p_2+p_3+p_4+p_5}\cdot \frac{p_3}{p_3+p_4+p_5} \end{equation}
(observe how this likelihood function, while informative about $p_4+p_5$, is uninformative about $p_4\left|p_4+p_5\right.$). We now use a technique due to \cite{hankin2010,hankin2017} and introduce fictional (reified) entities whose nonzero Bradley-Terry strength helps certain competitors or sets of competitors under certain conditions. The canonical example would be the home-ground advantage in association football. If players (teams) $1,2$ with strengths $p_1,p_2$ compete, and if our observation were $a$ home wins and $b$ away wins for team $1$, and $c$ home wins and $d$ away wins for team $2$, then a suitable likelihood function would be
[ \left(\frac{p_1+p_H}{p_1+p_2+p_H}\right)^a \left(\frac{p_1}{p_1+p_2+p_H}\right)^b \left(\frac{p_2+p_H}{p_1+p_2+p_H}\right)^c \left(\frac{p_2+p_H}{p_1+p_2+p_H}\right)^d ]
where $p_H$ is a quantification of the beneficial home ground effect.
Similar techniques have been used to account for the first-move
advantage in chess, and effective coordination between members of
doubles tennis teams; we may use a similar device to account for pole
position in Formula 1. Here I analyse seasons 2016-2019 using the
hyper2 package \citep{hankin2017} which implements the Plackett-Luce
likelihood function with additional reified entities.
Formula 1 dataset
Taking 2017 as an example, we see the rank table as follows:
\small \begin{CodeChunk} \begin{CodeOutput} AUS CHN BHR RUS ESP MON CAN AZE AUT GBR HUN BEL ITA SIN MAL JPN USA MEX BRA ABU Hamilton 2 1 2 4 1 7 1 5 4 1 4 1 1 1 2 1 1 9 4 2 Vettel 1 2 1 2 2 1 4 4 2 7 1 2 3 Ret 4 Ret 2 4 1 3 Bottas 3 6 3 1 Ret 4 2 2 1 2 3 5 2 3 5 4 5 2 2 1 Raikkonen 4 5 4 3 Ret 2 7 14 5 3 2 4 5 Ret DNS 5 3 3 3 4 Ricciardo Ret 4 5 Ret 3 3 3 1 3 5 Ret 3 4 2 3 3 Ret Ret 6 Ret ... Hartley 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 13 Ret Ret 15 Button 0 0 0 0 0 Ret 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Resta 0 0 0 0 0 0 0 0 0 0 Ret 0 0 0 0 0 0 0 0 0 \end{CodeOutput} \end{CodeChunk}
top <- 11
Each row is a driver and each column (after the first) a venue. We
see that Hamilton, the first row, came second in Australia, first in
China, second in Bahrain, fourth in Russia, and so on (Hartley,
Button, and Resta placed last). In the first column, we see the
result from Australia (AUS) in which Hamilton came second, Vettel
first, Bottas third, and so on. Notation used also includes six
different classes of no-score such as Ret for retired, WD for
withdrawn, and so on; a zero entry means "did not finish". It is
straightforward to translate this table into a Plackett-Luce
log-likelihood function using the hyper2 package; for simplicity we
will consider only the r top top-ranked drivers\footnote{In the
Plackett-Luce likelihood function, the performance of lower-ranked
players can be weakly informative about higher-ranked players'
strengths. For example, we see that Vettel retired twice (in
Singapore and Japan), so any player who placed in those venues will
effectively steal'' strength from Vettel, and generallygive'' it
to Hamilton or Bottas.}. Although it has many terms, the overall loglikelihood
expression is of the general form
library("hyper2",quietly=TRUE) library("magrittr",quietly=TRUE)
top <- 11 H2016 <- ordertable2supp(F1_table_2016[seq_len(top),]) H2017 <- ordertable2supp(F1_table_2017[seq_len(top),]) H2018 <- ordertable2supp(F1_table_2018[seq_len(top),]) H2019 <- ordertable2supp(F1_table_2019[seq_len(top),])
\newcommand{\pham}{p_\mathrm{Ham}} \newcommand{\pvet}{p_\mathrm{Vet}} \newcommand{\pbot}{p_\mathrm{Bot}} \newcommand{\prai}{p_\mathrm{Rai}} \newcommand{\pric}{p_\mathrm{Ric}} \newcommand{\pver}{p_\mathrm{Ver}} \newcommand{\pper}{p_\mathrm{Per}} \newcommand{\poco}{p_\mathrm{Oco}} \newcommand{\psai}{p_\mathrm{Sai}} \newcommand{\phul}{p_\mathrm{Hul}} \newcommand{\pmas}{p_\mathrm{Mas}}
\begin{equation} \frac{\pham^{20}\,\pmas^{16}\,\pbot^{19}\ldots }{\parbox{4in}{$(\pper+\poco)(\pver+\pper+\psai+\phul+\pmas)\ \rule{10mm}{0mm}(\pric+\pper+\pmas)(\psai+\pmas)\ \rule{20mm}{0mm} (\pver + \pper+ \poco + \psai+ \phul)\ldots $}} \end{equation}
Finding the maximum likelihood estimate for the players' strengths is
straightforward. The hyper2 package includes a suite of numerical
optimization routines, and because they have access to derivatives,
convergence is rapid. A graphical diagram of the strengths is given
in figure \ref{piechartstrength}
mH2016 <- maxp(H2016) mH2017 <- maxp(H2017) mH2018 <- maxp(H2018) mH2019 <- maxp(H2019)
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, seasons 2016-2019"} par(mfrow=c(2,2)) par(mar = c(2,2,2,2)/2) pie(mH2016[c(1,6,11,2,7,3,8,4,9,5,10)],cex=0.7,main="2016") pie(mH2017[c(1,6,11,2,7,3,8,4,9,5,10)],cex=0.7,main="2017") pie(mH2018[c(1,6,11,2,7,3,8,4,9,5,10)],cex=0.7,main="2018") pie(mH2019[c(1,6,11,2,7,3,8,4,9,5,10)],cex=0.7,main="2019")
```r
library("qcc")
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, 2016"} pareto.chart(mH2016[c(1,6,11,2,7,3,8,4,9,5,10)],main="2016",ylab="strength")
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, 2017"} pareto.chart(mH2017[c(1,6,11,2,7,3,8,4,9,5,10)],main="2017",ylab="strength")
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, 2018"} pareto.chart(mH2018[c(1,6,11,2,7,3,8,4,9,5,10)],main="2018",ylab="strength")
```r of the strengths of the top-ranked 11 drivers in Formula 1 motor racing, 2019"} pareto.chart(mH2019[c(1,6,11,2,7,3,8,4,9,5,10)],main="2019",ylab="strength")
We see in figure \ref{piechartstrength} that the driver with the
largest estimated strength was Rosberg in 2016 at about 30\%, and in
years 2017-8-9 was Hamilton (at about 29\%, 37\%, and 42\%
respectively). As an illustration of the value of likelihood methods
(as opposed to points-based methods), a likelihood ratio test
[samep.test(), supplied with the package] rejects the null that
Hamilton and Vettel have the same strength in 2018 ($H_0\colon \pham=
\pvet$) with a likelihood ratio of $e^{2.76}\simeq 15.8$,
corresponding (by Wilks's theorem) to an asymptotic $p$-value of about
$0.02$. We may also use the reified entity concept to test the null
hypothesis that Hamilton's strength remains constant 2016, where
Rosberg had the highest estimated strength, and 2017-9, where Hamilton
did; we fail to reject this null.
H2016a <- H2016 %>% psubs("Hamilton","Hamilton_pre") H2017a <- H2017 %>% psubs("Hamilton","Hamilton_post") H2018a <- H2018 %>% psubs("Hamilton","Hamilton_post") H2019a <- H2019 %>% psubs("Hamilton","Hamilton_post") Ha <- H2016a + H2017a + H2018a + H2019a
Do a statistical test:
ignorethis <- samep.test(Ha,c("Hamilton_pre","Hamilton_post")) ignorethis
Do another statistical test:
ignore <- samep.test(H2018,c("Hamilton","Vettel"))
Likelihood scoring vs points scoring
Applying the current points system, for example, to the 2017 results table we would rank the drivers as follows:
real <- ordertable2points(ranktable_to_ordertable(wikitable_to_ranktable(F1_table_2017)),formula1_points_systems(top)$real) real <- round(real[seq_len(top)],2) real <- sort(real,decreasing=TRUE) real[] <- seq_along(real)
\begin{multline} \mbox{Hamilton}\succ \mbox{Vettel}\succ \mbox{Bottas}\succ \mbox{Raikkonen}\succ\mbox{YES I KNOW} \mbox{Ricciardo}\succ \mbox{Verstappen}\succ\ \mbox{Perez1}\succ \mbox{Ocon}\succ \mbox{Sainz}\succ \mbox{Hulkenberg}\mbox{YES I KNOW}\succ \mbox{Massa} \end{multline}
(this happens to be identical to the ranking after the extra point was awarded for fastest lap). However, if we were award points according to Zipf's law we would have:
zipf <- ordertable2points(ranktable_to_ordertable(wikitable_to_ranktable(F1_table_2017)),formula1_points_systems(top)$zipf) zipf <- round(zipf[seq_len(top)],2) zipf <- sort(zipf,decreasing=TRUE) zipf[] <- seq_along(zipf)
\begin{multline} \mbox{Hamilton}\succ \mbox{Vettel}\succ \mbox{Bottas}\succ \mbox{Ricciardo}\succ \mbox{Verstappen}\succ \mbox{Raikkonen}\succ\mbox{YES I KNOW}\ \mbox{Perez}\succ \mbox{Ocon}\succ \mbox{Massa}\succ \mbox{Sainz}\succ \mbox{Hulkenberg} \end{multline}
Thus these two systems agree on the first three places but fourth is awarded to Raikkonen under the F1 system and Ricciardo under Zipf. Compare the likelihood ranking:
jj <- mH2017 jj <- sort(jj,decreasing=TRUE) jj[] <- seq_along(jj) jj
\begin{multline} \mbox{Hamilton}\succ \mbox{Vettel}\succ \mbox{Bottas}\succ \mbox{Raikkonen}\mbox{YES I KNOW}\succ \mbox{Ocon}\succ \mbox{Ricciardo}\succ\ \mbox{Verstappen}\succ \mbox{Perez}\succ \mbox{Massa}\succ \mbox{Sainz}\succ \mbox{Hulkenberg}\mbox{YES I KNOW} \end{multline}
Thus, for 2017 at least, we see that the current points system agrees with likelihood ranking for the top four places, while Zipf's law agrees to three.
p <- formula1_points_systems()[[1]] # real points system dp2016 <- ordertable2points(as.ordertable(F1_table_2016[seq_len(top),]),p) dp2016 <- sort(dp2016,decreasing=TRUE) dp2016[] <- seq_along(dp2016) m2016 <- maxp(H2016) m2016 <- sort(m2016,decreasing=TRUE) m2016[] <- seq_along(m2016)
dp2017 <- ordertable2points(as.ordertable(F1_table_2017[seq_len(top),]),p) dp2017 <- sort(dp2017,decreasing=TRUE) dp2017[] <- seq_along(dp2017) m2017 <- maxp(H2017) m2017 <- sort(m2017,decreasing=TRUE) m2017[] <- seq_along(m2017)
dp2018 <- ordertable2points(as.ordertable(F1_table_2018[seq_len(top),]),p) dp2018 <- sort(dp2018,decreasing=TRUE) dp2018[] <- seq_along(dp2018) m2018 <- maxp(H2018) m2018 <- sort(m2018,decreasing=TRUE) m2018[] <- seq_along(m2018)
dp2019 <- ordertable2points(as.ordertable(F1_table_2019[seq_len(top),]),p) dp2019 <- sort(dp2019,decreasing=TRUE) dp2019[] <- seq_along(dp2019) m2019 <- maxp(H2019) m2019 <- sort(m2019,decreasing=TRUE) m2019[] <- seq_along(m2019)
```r top-ranked 11 drivers in Formula 1 seasons 2016-2019. Horizontal axis gives official (points-based) order, and the vertical axis gives the likelihood horder. Thus, taking 2017 as an example, the points-based ordering would be Hamilton first, then Vettel, then R\"{a}ikk\"{o}nen; YES I KNOW IT IS WRONG: RMARKDOWN IS DEFECTIVE while the likelihood ordering is (reading vertically) Hamilton, Vettel, Bottas ALSO R PLOTTING SUCKS HERE AND I WILL HAVE TO CHANGE THE PDF BY HAND WITH INKSCAPE TO MAKE IT LOOK GOOD"}
par(mfrow=c(2,2)) par(mar = c(3,3,3,3))
par(pty='s') # square plot ox <- dp2016
print("m2016") dput(m2016) print("dp2016") dput(dp2016)
oy <- ordertrans(m2016,names(dp2016)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2016 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2017 oy <- ordertrans(m2017,names(dp2017)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2017 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2018 oy <- ordertrans(m2018,names(dp2018)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2018 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2019 oy <- ordertrans(m2019,names(dp2019)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2019 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
```r par(pty="s") ox <- dp2016 oy <- ordertrans(m2016,names(dp2016)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2016 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty="s") ox <- dp2017 oy <- ordertrans(m2017,names(dp2017)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2017 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty="s") ox <- dp2018 oy <- ordertrans(m2018,names(dp2018)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2018 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty="s") ox <- dp2019 oy <- ordertrans(m2019,names(dp2019)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2019 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
We can plot one ranking against the other\footnote{The historically correct points awarded to the drivers differs from that calculated here. That is for two reasons: firstly, ``fastest lap'' points are not included here, and also the truncation of the order table to the first 11 drivers can increase the rank of a driver if non-first-11 drivers are placed}, shown in Figure \ref{orderingbypointsandlikelihood}. We thus see a comparison between the two ordering systems. Taking 2017 as an example, drivers Hamilton and Vettel are are respectively first and second according to both ranking ranking procedures; but R\"{a}ikk\"{o}nen YES I KNOW IT IS WRONG: RMARKDOWN IS DEFECTIVE and Bottas are third and fourth, and fourth and third, according to points and likelihood respectively. We define the \emph{degree of agreement} between the two ranking systems as the maximal value of $r$ such that places $1,2,\ldots r$ all match. Thus, from figure \ref{orderingbypointsandlikelihood}, the degree of agreement between likelihood ranking and points ranking for the years 2016-2019 would be 8,2,6, and 2 respectively.
However, the points system used is essentially arbitrary. We could use, for example, Zipf's law to rank the drivers: award one point to the winner, half a point to second place, one third of a point to third, and so on; see figure \ref{orderingbypointsandlikelihood_zipf} in which we see generally poorer agreement between points-based ranks and likelihood-based ranks, with a degree of agreement of 4,2,2,2 for the years 2016-2019 respectively. This might be an indication that using a Zipf law points allocation is objectively worse than the current points system. There are a number of plausible points systems that might be used:
pz <- formula1_points_systems()$zipf dp2016z <- ordertable2points(as.ordertable(F1_table_2016[seq_len(top),]),pz) dp2016z <- sort(dp2016z,decreasing=TRUE) dp2016z[] <- seq_along(dp2016z) dp2017z <- ordertable2points(as.ordertable(F1_table_2017[seq_len(top),]),pz) dp2017z <- sort(dp2017z,decreasing=TRUE) dp2017z[] <- seq_along(dp2017z) dp2018z <- ordertable2points(as.ordertable(F1_table_2018[seq_len(top),]),pz) dp2018z <- sort(dp2018z,decreasing=TRUE) dp2018z[] <- seq_along(dp2018z) dp2019z <- ordertable2points(as.ordertable(F1_table_2019[seq_len(top),]),pz) dp2019z <- sort(dp2019z,decreasing=TRUE) dp2019z[] <- seq_along(dp2019z)
```r \label{orderingbypointsandlikelihood_zipf} but points ranking calculated by Zipf's law. REF DOES NOT WORK HERE Note the generally poorer agreement between the two ranking systems ALSO R PLOTTING SUCKS HERE AND I WILL HAVE TO CHANGE THE PDF BY HAND WITH INKSCAPE TO MAKE IT LOOK GOOD"} par(mfrow=c(2,2)) par(mar = c(3,3,3,3))
par(pty='s') # square plot ox <- dp2016z oy <- ordertrans(m2016,names(dp2016z)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2016 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2017z oy <- ordertrans(m2017,names(dp2017z)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2017 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2018z oy <- ordertrans(m2018,names(dp2018z)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2018 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2019z oy <- ordertrans(m2019,names(dp2019z)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order",ylab="likelihood order",main='2019 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
```r par(pty='s') # square plot ox <- dp2016z oy <- ordertrans(m2016,names(dp2016z)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order (Zipf)",ylab="likelihood order",main='2016 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2017z oy <- ordertrans(m2017,names(dp2017z)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order (Zipf)",ylab="likelihood order",main='2017 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2018z oy <- ordertrans(m2018,names(dp2018z)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order (Zipf)",ylab="likelihood order",main='2018 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray', cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
par(pty='s') # square plot ox <- dp2019z oy <- ordertrans(m2019,names(dp2019z)) plot(ox,oy,asp=1,pty='s',xlim=c(0,11.2),ylim=c(0,11.2),pch=16, xlab="points order (Zipf)",ylab="likelihood order",main='2019 season') par(xpd=TRUE) # allow drivers' names to appear outside plotting region for(i in seq_along(ox)){text(ox[i],oy[i],names(ox)[i],pos=4,col='gray',cex=0.6) } par(xpd=FALSE) # stop diagonal line from protruding beyond plotting region abline(0,1)
\begin{itemize} \item The current Formula 1 system $(25,18,15,12,10,8,6,4,2,1)$ \item The inaugural Formula 1 system $(8,6,4,3,2)$ \item Zipf's law $(1,\frac{1}{2},\frac{1}{3},\frac{1}{4},\ldots)$ \item Linear decrease $(n,n-1,n-2,\ldots,3,2,1,0)$ \item Halving system: $(1,\frac{1}{2},\frac{1}{4},\frac{1}{8},\ldots)$ \item A ``winner takes all'' system $(1,0,\ldots)$ \end{itemize}
We note that some of these may be generalized; we might consider a more general linear points system $(r,r-1,r-2,\ldots,3,2,1,0)$ for fixed integer $0<r<n$; the halving system can be generalized to a geometric distribution; and the winner takes all system can be replaced by giving equal points to the top $r$ competitors $(\underbrace{1,1,\ldots ,1}_{r},0,0,\ldots)$, for some integer $r<n$.
agreement2016 <- rep(NA,length(formula1_points_systems(top))) agreement2017 <- rep(NA,length(formula1_points_systems(top))) agreement2018 <- rep(NA,length(formula1_points_systems(top))) agreement2019 <- rep(NA,length(formula1_points_systems(top))) for(i in seq_along(formula1_points_systems(top))){ pa <- formula1_points_systems(top)[[i]] dp2016a <- ordertable2points(as.ordertable(F1_table_2016[seq_len(top),]),pa) dp2016a <- sort(dp2016a,decreasing=TRUE) dp2016a[] <- seq_along(dp2016a) agreement2016[i] <- sum(cumprod(names(dp2016a) == names(m2016))) dp2017a <- ordertable2points(as.ordertable(F1_table_2017[seq_len(top),]),pa) dp2017a <- sort(dp2017a,decreasing=TRUE) dp2017a[] <- seq_along(dp2017a) agreement2017[i] <- sum(cumprod(names(dp2017a) == names(m2017))) dp2018a <- ordertable2points(as.ordertable(F1_table_2018[seq_len(top),]),pa) dp2018a <- sort(dp2018a,decreasing=TRUE) dp2018a[] <- seq_along(dp2018a) agreement2018[i] <- sum(cumprod(names(dp2018a) == names(m2018))) dp2019a <- ordertable2points(as.ordertable(F1_table_2019[seq_len(top),]),pa) dp2019a <- sort(dp2019a,decreasing=TRUE) dp2019a[] <- seq_along(dp2019a) agreement2019[i] <- sum(cumprod(names(dp2019a) == names(m2019))) } # i loop closes agreement4yearsallp <- cbind( agreement2016, agreement2017, agreement2018, agreement2019 ) rownames(agreement4yearsallp) <- names(formula1_points_systems(top)) colnames(agreement4yearsallp) <- 2016:2019 wanted <- c(1,2,24,23,3,29)
It is straightforward to calculate the degree of agreement for the observed rank table for years 2016-2019, shown in the following table:
agreement4yearsallp[wanted,]
It is clear that there is no points system that is the best for all four years. However, observing that both the likelihood ranking and the points ranking are random variables in this paradigm suggests a method whereby we can objectively assess a given points system. Using sampling techniques we can repeatedly generate an order table \emph{in silico}, using estimated driver strengths from the observed tables. For each of, say, 1000 such synthetic tables, calculate drivers' maximum likelihood \citeauthor{plackett1975} strengths, and also their points awarded according to any given points system. We then compare rankings generated by the \citeauthor{plackett1975} strengths and the points awarded and note the degree of agreement between the two, as measured by the number of rankings correctly predicted. This furnishes an objective assessment of the points system used.
Numerical results
load("formula1_results_2016.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist rownames(OO)[3] <- "wta" OO2016 <- OO load("formula1_results_2017.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist rownames(OO)[3] <- "wta" OO2017 <- OO load("formula1_results_2018.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist rownames(OO)[3] <- "wta" OO2018 <- OO load("formula1_results_2019.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist rownames(OO)[3] <- "wta" OO2019 <- OO
summarytablemaker <- function(OO){ out <- data.frame( means = apply(OO,1,mean), winner_correct = rowSums(OO>1)/ncol(OO), all_correct = rowSums(OO==top)/ncol(OO) ) rownames(out) <- names(pointslist) return(out) } s2016 <- summarytablemaker(OO2016) s2017 <- summarytablemaker(OO2017) s2018 <- summarytablemaker(OO2018) s2019 <- summarytablemaker(OO2019)
```r systems when compared with a Placket-Luce strength ordering"} jj <- cbind(s2016[,1],s2017[,1],s2018[,1],s2019[,1]) # mean
typecol <- c(rep("red",2),rep("orange",10),rep("yellow",11),rep("green",1),rep("blue",7),rep("purple",7))
matplot(jj,type="n",xlab='points system',ylab='mean correct places')
points(s2016[,1],col=typecol,pch=1)
points(s2017[,1],col=typecol,pch=2)
points(s2018[,1],col=typecol,pch=3)
points(s2019[,1],col=typecol,pch=4)
legend("topright",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
matplot(jj[wanted,],type="n",xlab='points system',ylab='mean correct places',axes=FALSE) jjlabels <- names(formula1_points_systems(top)[wanted]) jjlabels[5] <- "wta" axis(1,at=1:6,labels=jjlabels) axis(2,at=1:4) points(s2016[wanted,1],pch=1,type='b',col='red') points(s2017[wanted,1],pch=2,type='b',col='orange') points(s2018[wanted,1],pch=3,type='b',col='green') points(s2019[wanted,1],pch=4,type='b',col='blue')
legend("topright",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
par(xpd=TRUE) text(x=6,y=s2016[wanted[6],1],pos=4,"2016",col='red') text(x=6,y=s2017[wanted[6],1],pos=4,"2017",col='orange') text(x=6,y=s2018[wanted[6],1],pos=4,"2018",col='green') text(x=6,y=s2019[wanted[6],1],pos=4,"2019",col='blue')
```r systems"} jj <- cbind(s2016[,2],s2017[,2],s2018[,2],s2019[,2]) # prob of winner #matplot(jj,type="n",xlab='points system',ylab='probability of winner correct') #points(s2016[,2],col=typecol,pch=1) #points(s2017[,2],col=typecol,pch=2) #points(s2018[,2],col=typecol,pch=3) #points(s2019[,2],col=typecol,pch=4) #legend("topright",lty=0,pch=1:4,legend=c("2016","2017","2018","2019")) matplot(jj[wanted,],type="n",xlab='points system',ylab='probability of winner correct',axes=FALSE) axis(1,at=1:6,labels=jjlabels) axis(2) points(s2016[wanted,2],pch=1,type='b',col='red') points(s2017[wanted,2],pch=2,type='b',col='orange') points(s2018[wanted,2],pch=3,type='b',col='green') points(s2019[wanted,2],pch=4,type='b',col='blue') #legend(x=3.4,y=0.67,lty=0,pch=1:4,legend=c("2016","2017","2018","2019")) par(xpd=TRUE) text(x=6,y=s2016[wanted[6],2],pos=4,"2016",col='red') text(x=6,y=s2017[wanted[6],2],pos=4,"2017",col='orange') text(x=6,y=s2018[wanted[6],2],pos=4,"2018",col='green') text(x=6,y=s2019[wanted[6],2],pos=4,"2019",col='blue')
```r systems"} jj <- cbind(s2016[,3],s2017[,3],s2018[,3],s2019[,3]) # total rank order correct
matplot(jj,type="n",xlab='points system',ylab='probability of total rank order correct')
points(s2016[,3],col=typecol,pch=1)
points(s2017[,3],col=typecol,pch=2)
points(s2018[,3],col=typecol,pch=3)
points(s2019[,3],col=typecol,pch=4)
legend("topright",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
matplot(jj[wanted,],type="n",xlab='points system',ylab='probability of total rank order correct',axes=FALSE) axis(1,at=1:6,labels=jjlabels) axis(2) points(s2016[wanted,3],pch=1,type='b',col='red') points(s2017[wanted,3],pch=2,type='b',col='orange') points(s2018[wanted,3],pch=3,type='b',col='green') points(s2019[wanted,3],pch=4,type='b',col='blue')
legend("topleft",lty=0,pch=1:4,legend=c("2016","2017","2018","2019"))
par(xpd=TRUE) text(x=5,y=s2016[wanted[5],3],pos=3,"2016",col='red') text(x=5,y=s2017[wanted[5],3],pos=1,"2017",col='orange') text(x=5,y=s2018[wanted[5],3],pos=3,"2018",col='green') text(x=5,y=s2019[wanted[5],3],pos=3,"2019",col='blue') ```
We now assess the six points systems using the methodology outlined above, using 1000 in silico trials. For each of the six points systems, each of the 1000 trials results in a single non-negative integer: the degree of agreement between the Plackett-Luce ranking and the rankings according to the points system considered. There are three measures that might be used to assess the distribution of degree of agreement: (1), the mean degree of agreement $\overline{d}$; (2), the probability of correctly predicting the winner, $\operatorname{Pr}(d\geqslant 1)$; and (3), the probability of correctly predicting the complete order statistic $\operatorname{Pr}(D=n)$. These summaries are shown for each of the six points systems in figures \ref{simulated_mean}, \ref{simulated_winner_correct} and \ref{simulated_complete_agreement} respectively.
We see that the linear points system $(n,n-1,n-2,\ldots,3,2,1)$ gives the highest mean number of places, and also the highest probability of correctly predicting the winner. However, the probability of predicting the complete order statistic is maximized using the winnner takes all system $(1,0,0,\ldots)$.
Discussion
Many competitive situations involve a ranking of the participants and one common method for combining repeated observations is to assign points for the winner, second place, third place, etc. The overall ranking is decided on the basis of accumulated points. However, points systems are unavoidably arbitrary. Assuming that participants attempt to optimize their placing means that the arbitrary nature of any points system.
Maximum likelihood estimation of participants' Plackett-Luce strengths furnishes a ranking system that does not suffer from the arbitrariness of a points system.
The points system used in Formula 1 motor racing is a source of lively debate from many perspectives, with changes being controversial. Here I treat total points scored as a random variable and compare different point allocation schemes against objective Plackett-Luce ranks. Of the six points systems considered here, a linear system or a winner-takes-all-system appear to be closest to objective Plackett-Luce, depending on the exact definition of ``closest''.
\bibliography{chess}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.