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)
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")

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.

```r
H2016a <- hyper2(brackets(H2016),d=powers(H2016),pnames=sub("Hamilton","Hamilton_pre", pnames(H2016)))
H2017a <- hyper2(brackets(H2017),d=powers(H2017),pnames=sub("Hamilton","Hamilton_post",pnames(H2017)))
H2018a <- hyper2(brackets(H2018),d=powers(H2018),pnames=sub("Hamilton","Hamilton_post",pnames(H2018)))
H2019a <- hyper2(brackets(H2019),d=powers(H2019),pnames=sub("Hamilton","Hamilton_post",pnames(H2019)))
Ha <- H2016a
Ha <- combine(Ha,H2017a)
Ha <- combine(Ha,H2018a) 
Ha <- combine(Ha,H2019a)
ignorethis <- samep.test(Ha,c("Hamilton_pre","Hamilton_post"))
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.

```r
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 oy <- ordertrans(m2016,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,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,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,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:

```r
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,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,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,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,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)

\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$.

```r
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
OO2016 <- OO

load("formula1_results_2017.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist
OO2017 <- OO

load("formula1_results_2018.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist
OO2018 <- OO

load("formula1_results_2019.rda") # this created by running f1points_Omaker.R; file contains R variables OO and pointslist
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) axis(1,at=1:6,labels=names(formula1_points_systems(top)[wanted])) 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=names(formula1_points_systems(top)[wanted]))
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=names(formula1_points_systems(top)[wanted])) 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}



Try the hyper2 package in your browser

Any scripts or data that you put into this service are public.

hyper2 documentation built on March 4, 2021, 9:09 a.m.