order to explore the region of the PBH parameter space which may be probed by
such experiments. We show that current bounds on the PBH abundance in the
range $1-10^{3}M_{\odot}$, coming from microlensing and the impact of
accretion of the CMB, could be potentially improved up to two orders of
magnitude by future 21 cm power spectrum observations.

  7. 7.

Removing Astrophysics in 21 cm maps with Neural Networks

The 21 cm line stands as a promising tool to map the distribution of matter in
the universe, given that it allows tracing the location of neutral Hydrogen,
which tracks the DM clumps. This may allow eventually achieving tomographic 3D
density maps of the IGM. However, the link between the HI and DM density
distribution is not unequivocal, since the brightness temperature strongly
relies on the thermal and ionization state of the medium. Furthermore, once
reionization is advanced, the signal from HII regions is mostly suppressed. In
Ref. [7], we attempt to recover the underlying matter density field from a 21
cm map, making use of machine learning methods. Specifically, a Convolutional
Neural Network (CNN) based on the so-called U-Net architecture [576] is
trained with 1000 numerical simulations of the brightness temperature field
and its corresponding density maps, aimed to learn how to recover the latter
from the former. These simulations include different initial conditions for
each realization, varying also three relevant astrophysical parameters, in
order to become robust to distinct astrophysical scenarios.

It is shown that, once it is trained, the CNN is able to successfully extract
the density field from a given 21 cm map. To test its accuracy, several
statistics of the matter density are computed, such as the power spectrum, the
cross-correlation coefficient or the probability distribution function,
matching the corresponding ones from the true maps within a few percent down
to scales $k\sim 2$ Mpc-1. The network can be applied to maps on different
redshift, although its performance worsens at low redshifts, $z\sim 10$, when
reionization is already ongoing and the link between the HI and the DM
distribution becomes more diffuse. Furthermore, part of the trained layers of
the network are employed to predict the underlying astrophysical parameters,
showing that the CNN actually learns astrophysical information, substracting
it from the 21 cm signal in order to predict the density field.

The remaining of this part in the original edition includes the articles as
they were published in the journals. Only the first pages of each publication
are kept for this version. The complete texts can be found in the official
thesis repository of the University of València, as well as in Refs. [1, 2, 3,
4, 5, 6, 7].

See pages 1 of Chapters/Papers/WDM1.pdf

See pages 1 of Chapters/Papers/WDM2.pdf

See pages 1 of Chapters/Papers/EarlyEoR.pdf

See pages 1 of Chapters/Papers/IDM1.pdf

See pages 1 of Chapters/Papers/IDM2.pdf

See pages 1 of Chapters/Papers/PBH.pdf

See pages 1 of Chapters/Papers/CNN.pdf

##  Part III Resumen de la tesis

### Resumen de la tesis

####  1 Introducción

Durante las últimas décadas, el enorme desarrollo de la cosmología nos ha
permitido alcanzar una profunda comprensión de nuestro universo, culminando en
el establecimiento del paradigma estándar cosmológico $\Lambda$CDM. Este
modelo asume diversos componentes para explicar la evolución del universo y
las estructuras que observamos. Por un lado, la conocida como energía oscura,
que conforma cerca de un 70% de la energía actual. Esta presenta un
comportamiento en la evolución cósmica equivalente a la constante cosmológica
$\Lambda$, lo que causa la expansión acelerada del universo detectada a partir
de supernovas lejanas. Por otra parte, diversas observaciones astrofísicas a
lo largo del siglo XX, especialmente de las curvas de rotación galácticas y de
cumúlos galácticos, hicieron evidente la presencia de un cierto tipo de
materia no interactiva ni radiativa, conocida como materia oscura (Dark
Matter, DM), que comprende cerca de un 25% de la energía cosmológica. Además,
la mayor parte de observaciones sugieren que esta materia es relativamente
fría, sin dispersión de velocidades térmica apreciable, siendo pues conocida
como materia oscura fría (Cold Dark Matter, CDM). Finalmente, la materia
bariónica del Modelo Estándar de física de partículas tan solo contribuye en
un 5% a la fracción energética actual. Este paradigma cosmológico es capaz de
explicar con precisión tanto observaciones de galaxias y supernovas lejanas
como las anisotropías del Fondo Cósmico de Microondas (Cosmic Microwave
Background, CMB) y la nucleosíntesis primordial. Pese a ello, la naturaleza de
la materia oscura y de la energía oscura permanecen esencialmente
desconocidas.

Las observaciones de las anisotropías del CMB, el remanente fósil de fotones
del universo primigenio, nos ofrecen información muy precisa sobre los
diversos componentes del universo temprano, así como acerca del crecimiento de
las perturbaciones primordiales. Por otra parte, los catálogos galácticos han
permitido crear mapas cósmicos del universo cercano con gran detalle, los
cuales permiten conocer cómo se distribuye la materia e inferir los procesos
que han llevado a la formación de las estructuras cosmológicas que vemos
actualmente. Sin embargo, entre el desacoplamiento del CMB y las galaxias de
nuestro vecindario, existe un gran intervalo temporal esencialmente incierto y
desconocido. Durante la época de dominación de materia, las regiones más
densas del universo comienzan a colapsar formando las primeras estructuras
cosmológicas. Este proceso se acrecienta durante las Edades Oscuras, periodo
tras la Recombinación conocido así debido a la ausencia de fuentes de luz
(estelares y primordiales) que permitan observarlo directamente. Con el paso
del tiempo, llega un momento en que las inhomogeneidades son tan grandes que
producen halos suficientemente masivos como para que la materia de su interior
se enfríe y colapse formando estrellas. La época del nacimiento de las
primeras estrellas se conoce como el Amanecer Cósmico, e implica la emisión
por parte de las primeras galaxias de gran cantidad de radiación energética,
en forma de rayos UV y X, que puede ser absorbida por los átomos del Medio
Intergaláctico (Intergalactic Medium, IGM). Esto se espera que caliente el gas
difuso presente entre las galaxias y lo ionice. Cuando la formación estelar es
lo suficientemente alta, el medio circundante alrededor de las galaxias se
ioniza completamente, formando burbujas de hidrógeno ionizado que se expanden
durante la Época de Reionización (Epoch of Reionization, EoR), hasta que todo
el medio intergaláctico alcance un estado ionizado. Si bien esta cronología
cósmica es generalmente aceptada, todavía son desconocidos los detalles,
duración y extensión de cada periodo. Aunque la observación del espectro de
cuásares y galaxias, así como la polarización del CMB, pueden proporcionar
cierta información relativa a la extensión y la consumación de la EoR, aún
carecemos de una evaluación precisa de la historia de ionización del medio
intergaláctico.

La señal cosmológica de 21 cm se constituye como una prometedora herramienta
para entender la evolución del universo. Esta radiación emerge de la
transición hiperfina del hidrógeno neutro, y dada la gran abundancia de este
elemento en el universo, podría permitir realizar un cartografiado de la
distribución de materia en función del tiempo. Dada su sensibilidad respecto
al estado térmico del medio, su detección aportaría valiosa información sobre
la época de formación de primeras galaxias durante el Amanecer Cósmico.
Además, dado que solo está presente en átomos neutros, proporcionaría un medio
para poder conocer con precisión la historia de ionización durante la EoR. Por
otra parte, la señal de 21 cm podría servir como una herramienta excepcional
para ahondar en la naturaleza de la materia oscura. Aunque el modelo de
materia oscura fría es el aceptado generalmente, existen otros candidatos bien
motivados teóricamente que podrían proporcionar la materia oscura presente en
el universo. Estos modelos de materia oscura alternativos podrían producir un
impacto en la evolución del medio intergaláctico y tener por tanto efectos
observables. Es por ello que la precisa determinación de las distintas épocas
de la cronología cósmica resulta de gran importancia para acotar los distintos
escenarios posibles de materia oscura. En esta tesis, se examina el potencial
de la señal de 21 cm, junto con otras medidas del estado de ionización del
medio intergaláctico, con el fin de acotar la composición, naturaleza y
evolución de la materia oscura.

####  2 Candidatos alternativos de materia oscura

El caso paradigmático de materia oscura fría considera partículas masivas
débilmente interactivas (Weakly Interacting Massive Particles, WIMP), con
masas usualmente del orden del GeV o superior, que estuvieron en equilibrio
con el plasma cósmico en el universo primigenio. Estas interacciones dejarían
de ser eficientes para mantener el acoplamiento con las partículas del Modelo
Estándar en cierto momento, dejando un remanente que no interactúa con el
resto de materia bariónica y radiación. Si bien esta hipótesis está bien
motivada y es capaz de aportar la abundancia de materia oscura observada, las
partículas constituyentes no han sido observadas en experimentos dedicados a
su detección. Por otra parte, modelos de materia oscura fría presentan
diversas discordancias con observaciones cosmológicas a escalas pequeñas, que
pueden resumirse en el hecho de que simulaciones de materia oscura fría
predicen más subestructuras de las halladas en las observaciones. Pese a que
estos problemas podrían ser solucionados recurriendo a ciertos efectos
astrofísicos, como retroalimentación estelar o de supernovas, otra posibilidad
consiste en considerar modelos alternativos de materia oscura capaces de
solucionar dichas discrepancias.

Una de las propuestas más atractivas se basa en la materia oscura templada
(Warm Dark Matter, WDM), cuyas partículas tendrían una masa del orden del keV.
Ejemplos de este tipo de materia oscura podrían ser los gravitinos de teorías
supersimétricas, o bien neutrinos estériles, presentes en muchas extensiones
del Modelo Estándar de partículas. Debido a su masa, su dispersión de
velocidades y su temperatura serían notablemente mayores a las del caso frío.
Esto sería relevante en el crecimiento de las fluctuaciones de materia,
actuando como una presión efectiva, y contrarrestando el colapso gravitatorio
en escalas pequeñas. Como consecuencia, estructuras con tamaños por debajo de
una cierta escala serían fundamentalmente suprimidas. Esta escala actuaría
como una escala de Jeans efectiva, y estaría determinada principalmente por la
masa de la materia oscura templada y su abundancia. En estos modelos, la
cantidad de halos de masas bajas sería reducida, explicando así algunas de las
discrepancias presentes con la materia oscura fría. Este hecho tendría a su
vez un impacto en la evolución del medio intergaláctico, ya que produciría
menos galaxias de las esperadas, retrasando ciertos periodos de la cronología,
como el amanecer cósmico o la reionización. Por ello, es de esperar que el
análisis de observables del estado de ionización y de la señal de 21 cm sean
sensibles a dichos efectos, y puedan potencialmente constreñir modelos de
materia oscura templada.

Otra candidato alternativo podría tratarse de materia oscura interactiva
(Interacting Dark Matter, IDM) con radiación. Pese a que hay cotas muy
restrictivas respecto a las interacciones de la materia oscura con fotones, no
es posible excluirlas por completo. Colisiones elásticas de partículas de
materia oscura interactiva con fotones en el universo temprano producirían un
acoplamiento de los primeros a los segundos, reduciendo las fluctuaciones de
la materia oscura interactiva a escalas pequeñas, e imprimiéndoles en cierta
manera el espectro oscilatorio típico de las perturbaciones adiabáticas de la
radiación. En consecuencia, esto resultaría en una supresión del número de
estructuras cosmológicas en pequeñas escalas, de una manera similar al caso de
materia oscura templada, lo que ha sido propuesto como solución a algunas de
las discordancias del escenario de materia oscura fría con los datos. De la
misma manera, el número de halos que albergase formación estelar se reduciría,
y por tanto, observaciones de la absorción y emisión en 21 cm y límites sobre
la historia de ionización serían capaces de acotar el espacio de parámetros de
la materia oscura interactiva compatible con los datos.

Finalmente, cabe la posibilidad de que los constituyentes de la materia oscura
no sean partículas microscópicas, sino objetos astrofísicos macroscópicos
masivos, conocidos como MACHOs. La mayoría de estos candidatos, como podrían
ser planetas gigantes y enanas marrones, son descartados por no ser capaces de
originar las fluctuaciones primordiales patentes en el CMB y la cantidad de
bariones en la nucleosíntesis primordial. No obstante, el colapso de agujeros
negros primordiales (Primordial Black Holes, PBHs) durante el universo
temprano sería consistente con dichas restricciones, por su naturaleza previa
a la formación estelar, y podrían formar parte de al menos una fracción de la
materia oscura. Esta hipótesis recobró especial atención tras las
observaciones por parte de las colaboraciones LIGO y VIRGO de diversos eventos
de emisión de ondas gravitatorias originadas en la fusión de agujeros negros
de varias masas solares [41]. Contrariamente a los agujeros negros
astrofísicos estándares, limitados a masas iguales o superiores a $\sim
3M_{\odot}$, los agujeros negros primordiales podrían ser formados en
principio con cualquier masa. Estos objetos podrían además actuar como germen
de los agujeros negros supermasivos observados en los centros de la mayoría de
las galaxias, cuyas masas son difícilmente explicables mediante mecanismos de
acrecimiento estándares, a no ser que partiesen de una masa inicial anterior a
la época de formación estelar. Uno de los aspectos más interesantes de los
agujeros negros primordiales como candidatos de la materia oscura es la gran
cantidad de efectos observables que su existencia implicaría, permitiendo
acotar su abundancia hasta niveles potencialmente muy bajos. Ejemplos de ello
son efectos de lente gravitatoria, la emisión de partículas debido a su
evaporación, la inyección de energía en el medio circundante por el
acrecimiento de materia, o el incremento en las fluctuaciones primordiales
debido a su carácter discreto. Estos dos últimos efectos podrían afectar
potencialmente a la evolución del medio intergaláctico, y por tanto, a la
línea de 21 cm.

####  3 Metodología

La mayor parte de la tesis se focaliza en explorar diversos modelos de materia
oscura alternativos al estándar de materia oscura fría mediante simulaciones
de la evolución del universo, con el fin de contrastarlos con datos
observacionales para acotar sus propiedades. Para ello, se han llevado a cabo
simulaciones del crecimiento de las fluctuaciones del campo de materia y la
formación de halos, así como de la evolución térmica, y de ionización del
medio intergaláctico. En la mayor parte de los artículos incluidos, se ha
hecho uso del programa de acceso libre 21cmFAST [378], que emplea una serie de
aproximaciones semi-analíticas para calcular la evolución del medio
intergaláctico. Esencialmente, este código resuelve la evolución del campo de
densidad empleando teoría de perturbaciones lagrangiana de segundo orden. A
partir del mapa de densidad, el programa calcula las regiones ionizadas
alrededor de las zonas con mayor densidad, siguiendo la prescripción de
Furlanetto-Zaldarriaga-Hernquist (FZH) [509], la cual se basa en la aplicación
del formalismo de conjunto de excursiones (excursion set formalism) [175].
Estas simulaciones concuerdan con cálculos más rigurosos empleando
simulaciones hidrodinámicas o de N cuerpos para escalas relativamente grandes,
a partir de $\sim 1$ Mpc. Sin embargo, para el análisis presentado en la Ref.
[2], se han empleado simulaciones hidrodinámicas más precisas, basadas en el
programa de Reionización Cósmica en Computadoras (Cosmic Reionization on
Computers, CROC) [514], que emplean un tratamiento preciso de la evolución de
la densidad, del transporte radiativo, y por tanto, un cálculo más robusto del
campo de ionización. Ambos métodos asumen diversos modelos astrofísicos cuyos
parámetros no están bien determinados por las observaciones, y por tanto, los
resultados finales son sensibles a la elección específica de dichos
parámetros.

Para el análisis de resultados y el contraste con datos observacionales, se
emplea tanto un enfoque probabilista frecuentista, basado en el máximo de la
función de verosimilitud (o, equivalentemente, en el mínimo del $\chi^{2}$),
como métodos bayesianos, muestreando la función de probabilidad posterior
mediante Cadenas de Markov de Monte Carlo (Monte Carlo Markov Chains, MCMC).
Finalmente, en la Ref. [7], se hacen uso de métodos de aprendizaje profundo
(deep learning), empleando redes neuronales convolucionales.

####  4 Resultados de la tesis

Esta sección comprende un resumen de los trabajos de investigación originales
que conforman el cuerpo de la tesis doctoral. A continuación, se presentan la
motivación y los resultados más relevantes de cada artículo.

#####  4.1 Materia oscura cálida e historia de ionización del Universo

En la Ref. [1] se estudia el impacto de modelos de materia oscura templada en
la evolución térmica y de ionización del medio intergaláctico. Debido al
efecto de dispersión de velocidades de partículas con masas $\sim$ keV, la
formación de estructuras es suprimida en escalas pequeñas con respecto al caso
estándar de materia oscura fría. Esto implica una reducción en el número de
halos de baja masa que pueden albergar formación estelar, capaces de emitir
radiación UV y rayos X responsables de ionizar y calentar el medio
intergaláctico. La disminución del número de tales fuentes conduce a un
retraso del amanecer cósmico, de la época de calentamiento y de la EoR.
Haciendo uso de simulaciones semi-analíticas, la evolución de la temperatura
cinética y de la fracción de ionización se calcula en varios escenarios de
materia oscura templada, variando la masa de sus partículas entre $1$ y $4$
keV, un rango motivado físicamente para tratar de resolver los problemas a
pequeña escala presentes con la materia oscura fría. Dado que la historia
térmica del medio intergaláctico es muy sensible a los procesos estelares y
radiativos involucrados, también consideramos diferentes escenarios
astrofísicos. Se modifican tres parámetros astrofísicos: la temperatura virial
mínima (relacionada con la masa mínima del halo para producir formación
estelar), la eficiencia de ionización de rayos UV y la eficiencia de
calentamiento de los rayos X (relacionada con la luminosidad en rayos X de
fuentes astrofísicas, como binarias de rayos X).

Los resultados obtenidos en lo que respecta a la fracción de ionización de las
simulaciones se han confrontado con límites observacionales sobre el estado de
ionización durante la EoR. Entre los tipos de datos tenidos en cuenta, se
hallan las restricciones y los límites inferiores de la fracción de ionización
en el rango de desplazamientos al rojo (redshifts) $8\geq z\gtrsim 6$ de
espectros de cuásares y las galaxias emisoras de radiación Ly$\alpha$. Por
otro lado, también empleamos la profundidad óptica de Thomson inferida de las
anisotropías del CMB, obtenida mediante la integral de la fracción de
ionización sobre el redshift, en lugar de una restricción dependiente del
corrimiento al rojo. Se muestra que los límites que emplean solo datos de CMB
son menos restrictivos que los obtenidos usando el estado de ionización a un
tiempo dado, dado que la profundidad óptica de Thomson es una cantidad
integrada, y no una función del tiempo.

Encontramos que algunos parámetros astrofísicos están fuertemente degenerados
entre ellos, y a su vez, con la masa de la materia oscura templada. Más
específicamente, el retraso de la EoR presente en un modelo de materia oscura
templada puede ser contrarrestado aproximadamente por una alta eficiencia de
ionización (que impulsa una reionización más temprana) o por una temperatura
virial mínima más baja (lo que permite que halos de menor masa formen
estrellas y contribuyan a emitir más radiación ionizante). Por esa razón, los
límites inferiores de la masa de la materia oscura templada obtenidos por este
método son más débiles que aquellos inferidos a partir del bosque de
Ly$\alpha$, aunque son capaces de limitar masas de la materia oscura templada
menores de $\sim 1$ keV.

#####  4.2 Materia oscura cálida y reionización cósmica

Con el fin de profundizar en la fenomenología de la materia oscura templada,
en la Ref. [2] se continúa investigando el impacto de estos modelos en la
evolución cosmológica y la formación de estructuras. En este caso, se emplean
precisas simulaciones hidrodinámicas, en lugar de cálculos semi-analíticos,
proporcionando así resultados más fiables. Hacemos uso de simulaciones basadas
en el programa Cosmic Reionization on Computers (CROC) [514], que representan
la vanguardia en el modelado numérico de la EoR. En las condiciones iniciales
de cada simulación, el espectro de potencias es modificado adecuadamente para
incluir la supresión característica de escenarios de materia oscura templada.
Se asume una masa de las partículas de $m_{\rm WDM}=3$ keV, valor que se
encuentra alrededor de las restricciones obtenidas mediante el bosque de
Ly$\alpha$.

Se calculan diferentes observables de reionización para explorar el impacto de
la escala de supresión de estructuras, como la función de masa de los halos y
las funciones de luminosidad UV de la galaxia, que muestran una supresión (en
masas bajas y magnitudes altas, respectivamente) con respecto al caso
estándar. Los telescopios actuales no parecen ser capaces de discriminar entre
estos casos de materia oscura fría y templada. Sin embargo, el futuro
telescopio James Webb Space Telescope (JWST) sería capaz de alcanzar el
extremo de la función de luminosidad de galaxias de magnitudes débiles, siendo
apto para restringir tales masas de la materia oscura templada. También se
discute el efecto en la función de distribución de probabilidad cumulativa
para la profundidad óptica Ly$\alpha$ efectiva $\tau_{\rm GP}$ en el medio
intergaláctico posterior a la reionización. En comparación con el escenario de
materia oscura fría en un tiempo dado, el caso de materia oscura templada
presenta distribuciones con colas más extendidas hacia valores de $\tau_{\rm
GP}$ altos, que corresponden a regiones de alta densidad, ya que se ionizan
posteriormente. Comparando con los datos observacionales del bosque de
Ly$\alpha$, el caso de materia oscura fría resulta preferido respecto al caso
de materia oscura templada, aunque no está claro si el ajuste de los
parámetros astrofísicos en el escenario de materia oscura templada podría
mejorar su concordancia con los datos. Por último, también se calcula el
espectro de potencias de la línea de 21 cm, que en el caso de materia oscura
templada muestra un balance entre una supresión debida a la dispersión térmica
de velocidades, y un realzamiento dado por la reionización retardada. Esto
conduce a un comportamiento no trivial, que lo hace difícilmente detectable
por la primera generación de experimentos de 21 cm.

#####  4.3 ¿Hubo un componente de reionización temprano en nuestro universo?

La historia de ionización del universo durante la EoR está lejos de
encontrarse bien determinada. Aunque la mayoría de los límites son
consistentes con un proceso de reionización relativamente rápido que termina a
redshifts $z\sim 6$, aún puede ser posible una contribución temprana, por
ejemplo, de escenarios de materia oscura no canónicos que inyecten energía en
el medio, como materia oscura desintegrándose o agujeros negros primordiales.
En la Ref. [3], esta cuestión se examina a la luz de los datos de Planck del
CMB. La historia de ionización se modela para que sea coherente con el
espectro de polarización EE del CMB, que es sensible a la reionización en
multipolos bajos. Se emplean varios métodos para emular la evolución de la
fracción de ionización $x_{e}$. En primer lugar, se consideran dos
parametrizaciones funcionales frecuentes en la literatura, como la
parametrización simétrica y la asimétrica en $z$. Por otro lado, se aplica una
reconstrucción de $x_{e}$ basada en el Análisis de Componentes Principales
(Principal Component Analysis, PCA), extrayendo los modos propios de la matriz
de Fisher de los componentes $C_{l}$ del espectro EE. Se espera que aquellos
con las varianzas más pequeñas contengan la información cosmológica más
valiosa. Un subconjunto de estas funciones propias se emplea como base para
escribir como su combinación el historial de ionización alrededor de un modelo
por defecto dado, ya sea asumiendo un valor constante o una de las
parametrizaciones mencionadas anteriormente. Finalmente, se emplea un tercer
método basado en el polinomio de interpolación de Hermite cúbico por partes
(Piecewise Cubic Hermite Interpolating Polynomial, PCHIP), considerando como
parámetros libres la fracción ionizada en varios $z$ diferentes, que se
interpolan suavemente a través del PCHIP.

Comparamos los diferentes enfoques para recuperar la evolución de la
ionización con un análisis empleando Cadenas de Markov de Monte Carlo (MCMC),
que muestra que algunos de ellos prefieren un componente de reionización
temprano, incluido algún patrón oscilatorio, como el uso de la combinación de
PCA sumada a un modelo por defecto. Sin embargo, cuando se incluye la
combinación PCA como argumento de la parametrización, explotando su forma
funcional, la evolución de $x_{e}$ resultante es más suave y no muestra un
inicio temprano. La aplicación de los criterios frecuentista de información de
Akaike y de información bayesiano confirman que ninguna de las formulaciones
en particular es preferida por los datos sobre el resto, y por lo tanto, que
los datos de CMB no indican una reionización temprana, contrariamente a
afirmaciones anteriores en la literatura [554].

#####  4.4 Una nueva mirada al escenario de la materia oscura interactiva

En la Ref. [4], los modelos fenomenológicos de materia oscura que involucran
la dispersión elástica con fotones, conocidos como materia oscura interactiva,
se discuten en el contexto de la evolución térmica del medio intergaláctico y
del número de galaxias satélites de la Vía Láctea. En el caso de la materia
oscura interactiva, el acoplamiento entre partículas de materia oscura y
fotones produce un efecto de amortiguación debido a las colisiones, induciendo
una supresión y un patrón oscilatorio en las escalas más pequeñas en el
espectro de potencias, de manera similar a los inducidos por el acoplamiento
barión-fotón en el caso estándar. Este efecto puede reducir el número de
estructuras a pequeña escala y halos de baja masa, de manera similar a lo que
ocurre con la materia oscura templada. La magnitud de este amortiguamiento
está determinada por la sección eficaz de dispersión elástica sobre la masa,
$\sigma_{\gamma\rm DM}/m_{\rm DM}$.

Se han realizado simulaciones semi-analíticas del estado térmico y de
ionización del medio intergaláctico para investigar el impacto de escenarios
con materia oscura interactiva. La supresión del crecimiento de estructuras
debido a la amortiguación de colisiones retrasa la formación de las primeras
estrellas y galaxias, retardando todos los procesos astrofísicos como el
calentamiento por rayos X y la reionización. Análogamente a lo visto en la
sección 4.1 (Ref. [1]) con respecto a la masa de la materia oscura templada,
la sección eficaz de la materia oscura interactiva está fuertemente degenerada
con varios parámetros astrofísicos, principalmente con la eficiencia de
ionización y la temperatura virial mínima. Por otro lado, se calcula el número
de subhalos esperados, que en un escenario de materia oscura interactiva se
reduce con respecto al escenario de materia oscura fría. Este hecho permite
acotar la sección eficaz elástica al comparar las predicciones teóricas con
las galaxias satélite de la Vía Láctea observadas. Sin embargo, la restricción
a partir de las galaxias satélite depende en gran medida de la masa total
asumida de la Vía Láctea, que no está completamente determinada con certeza.
Dentro de su rango de masa esperado, los límites de la sección eficaz pueden
variar hasta un orden de magnitud. Combinando este límite superior con los
límites de la historia de ionización del medio intergaláctico, se pueden
romper algunas degeneraciones astrofísicas.

#####  4.5 Microfísica de materia oscura y observaciones de 21 cm

La colaboración de EDGES presentó la supuesta primera detección de la señal
global de 21 cm, presentando un espectro de absorción con una amplitud de
aproximadamente $500$ mK, que duplica las expectativas en el caso estándar
$\Lambda$CDM. Esto resulta ser difícilmente explicable solo a partir de
distintos modelos astrofísicos, posiblemente requiriendo escenarios no
estándar con un mecanismo de enfriamiento adicional para explicar esta
absorción, como materia oscura interaccionando con bariones. Sin embargo, la
localización en el tiempo de la señal de absorción medida por EDGES,
perfectamente consistente con las predicciones estándar en escenarios
$\Lambda$CDM, se puede emplear para restringir otros escenarios no canónicos
de materia oscura que pueden afectar la evolución del medio intergaláctico. En
la Ref. [5], se explotan dos aspectos de la ubicación temporal de la señal de
EDGES para restringir modelos de materia oscura interactiva y templada, que
son capaces de retrasar el inicio de la formación estelar durante el amanecer
cósmico.

Por un lado, se exige que el mínimo de la temperatura global de brillo de 21
cm calculada a partir de simulaciones aparezca no más tarde del pico de
absorción de EDGES, situado en $z\simeq 17$. Por otro lado, al examinar el
comienzo del régimen de absorción en los datos de EDGES, se puede inferir
cuándo el acoplamiento Ly$\alpha$ debe ser lo suficientemente fuerte como para
acoplar la temperatura de espín con la temperatura cinética del medio
intergaláctico. Siguiendo estas ideas, para ser coherente con los datos, se
puede suponer que el coeficiente de acoplamiento Ly$\alpha$ debe alcanzar el
orden uno a $z\simeq 20$. Por tanto, es posible descartar los modelos materia
oscura interactiva y templada que presentan una época de absorción demasiado
retrasada y no pueden cumplir las condiciones anteriores. Ambos criterios
mencionados permiten extraer restricciones competitivas en la masa de la
materia oscura templada y en la sección eficaz de la interactiva. Estos
límites son fuertemente dependientes del escenario astrofísico específico,
debilitando los límites en escenarios con bajas temperaturas viriales. No
obstante, asumiendo el escenario predeterminado de estrellas de población II
(de baja metalicidad) con una temperatura virial de $10^{4}$ K como fuente de
radiación Ly$\alpha$, los límites obtenidos mediante este método pueden
mejorar en gran medida aquellos obtenidos a partir de otros análisis.

#####  4.6 Restringiendo la abundancia de agujeros negros primordiales con
cosmología de 21 cm

En la Ref. [6] se estudia en detalle los efectos cosmológicos de la existencia
de agujeros negros primordiales como parte de la materia oscura en la señal
cosmológica de 21 cm. Hay dos consecuencias principales de agujeros negros de
origen primitivo de masa solar o superior para la evolución del medio
intergaláctico. Por un lado, la presencia de agujeros negros primordiales en
el universo temprano implicaría una contribución adicional al espectro de
potencias de la materia dado por un espectro de ruido blanco de Poisson,
independiente de la escala, debido a su naturaleza discreta, evolucionando
como un modo de isocurvatura. Esto magnificaría la formación de estructuras a
pequeña escala, donde este término se vuelve dominante con respecto a las
fluctuaciones adiabáticas estándares. Por lo tanto, este incremento afectaría
principalmente a los minihalos, aquellos halos de baja masa que no son lo
suficientemente masivos para iniciar la formación de estrellas. La temperatura
de brillo promedio de 21 cm esperada proveniente de estos minihalos es
calculada, que difiere de la señal global estándar del medio intergaláctico.
En el artículo se demuestra que, incluso con el incremento en el espectro de
potencias a pequeñas escalas debido a la contribución de Poisson, no sería
suficiente para producir una señal detectable, contrariamente a lo que se ha
afirmado anteriormente en la literatura [316, 317].

Por otro lado, al igual que los agujeros negros astrofísicos estándares, los
agujeros negros primordiales acrecerían materia circundante debido a su
atracción gravitatoria. Durante la caída hacia el PBH, las partículas podrían
acelerarse enormemente, liberando radiación UV y de rayos X al ambiente
circundante. Estos fotones inyectarían energía en el medio cuando se
absorbieran, calentándolo e ionizándolo aún más. En el artículo, calculamos la
señal de 21 cm que surge del entorno local de los PBH, mostrando que estaría
muy por debajo de las sensibilidades observables. Sin embargo, la radiación
energética puede escapar de los alrededores cercanos de los agujeros negros
primordiales y depositarse en el medio intergaláctico, afectando notablemente
la señal global de 21 cm y el espectro de potencias. Debido a la ionización
adicional y al medio más caliente, se suprimiría la amplitud de absorción de
la temperatura de brillo y sus fluctuaciones, implicando una transición a un
régimen de emisión en los casos más extremos. Realizamos un análisis MCMC con
simulaciones semi-analíticas de la señal del medio intergaláctico de 21 cm,
variando la fracción de agujeros negros primordiales como materia oscura,
$f_{\rm PBH}$, y su masa, $M_{\rm PBH}$, junto con cuatro parámetros
astrofísicos. Las sensibilidades esperadas de los próximos interferómetros
HERA y SKA se emplean para explorar la región del espacio de parámetros de
agujeros negros primordiales que puede ser testada por tales experimentos.
Finalmente, mostramos que los límites actuales en la abundancia de agujeros
negros primordiales en el rango de $1-10^{3}M_{\odot}$, provenientes del
efecto de lente gravitatoria y el impacto de la acrecimiento en el CMB,
podrían mejorarse potencialmente hasta en dos órdenes de magnitud con
observaciones del espectro de potencias de 21 cm.

#####  4.7 Eliminando la astrofísica en mapas de 21 cm con redes neuronales

La línea de 21 cm se erige como una prometedora herramienta para mapear la
distribución de materia en el universo, dado que permite localizar la
ubicación del hidrógeno neutro en los cúmulos de materia oscura. Esto podría
permitir lograr mapas tomográficos de densidad en 3D del medio intergaláctico.
Sin embargo, la relación entre la distribución de hidrógeno neutro y densidad
de materia oscura no es unívoca, ya que la temperatura de brillo depende en
gran medida del estado térmico y de ionización del medio. Además, una vez que
la reionización es muy avanzada, la señal de las regiones de hidrógeno
ionizado se suprime en su mayor parte. En la Ref. [7], se trata de recuperar
el campo de densidad de materia subyacente a partir de un mapa de 21 cm,
utilizando métodos de aprendizaje automático. En concreto, se entrena una red
neural convolucional basada en la arquitectura U-Net [576] con 1000
simulaciones numéricas del campo de temperatura de brillo y sus
correspondientes mapas de densidad, con el objetivo de aprender a recuperar
este último a partir del mapa de 21 cm. Estas simulaciones incluyen diferentes
condiciones iniciales para cada realización, variando también tres parámetros
astrofísicos relevantes, con el fin de volverse robustos a distintos
escenarios astrofísicos.

Se demuestra que, una vez entrenada, la red neuronal es capaz de extraer con
éxito el campo de densidad de un mapa de 21 cm dado. Para probar su precisión,
se calculan varios estadísticos del campo de densidad de materia, como el
espectro de potencias, el coeficiente de correlación cruzada o la función de
distribución de probabilidad, coincidiendo las predicciones con las
correspondientes de los mapas reales con alta precisión hasta escalas tan
pequeñas como $k\sim 2$ Mpc -1. La red se puede aplicar a mapas en diferentes
momentos de la evolución cósmica, aunque su rendimiento empeora con instantes
más recientes, $z\sim 10$, cuando la reionización ya está en curso y el
vínculo entre la distribución de hidrógeno neutro y materia oscura se vuelve
más difuso. Además, parte de las capas entrenadas de la red se emplean para
predecir los parámetros astrofísicos subyacentes, lo que muestra que la red
neuronal realmente aprende información astrofísica, con el fin de sustraerla
de la señal de 21 cm para predecir el campo de densidad.

####  5 Conclusiones

Los trabajos de investigación incluidos en la presente tesis doctoral
pretenden ejemplificar diversos enfoques para estudiar la naturaleza y
propiedades de la materia oscura a partir del estudio de la evolución del
medio intergaláctico. Diferentes modelos alternativos de materia oscura son
considerados, como materia oscura templada, interactiva y agujeros negros
primordiales, los cuales presentan diferentes efectos en la formación de
estructuras y en la historia térmica. Los casos de materia oscura templada e
interactiva presentan una supresión en las fluctuaciones a pequeñas escalas,
dada por la dispersión de velocidades térmica y por el amortiguamiento
colisional, respectivamente, retrasando la formación de estructuras
cosmológicas y con ello todos los procesos astrofísicos. Por otra parte,
agujeros negros primordiales de masa solar o mayor pueden acrecer materia,
radiando e inyectando posteriormente energía en el medio, calentándolo e
ionizándolo. Todos estos efectos dejarían un gran impacto en la evolución del
medio intergaláctico, afectando a la época de formación de primeras estrellas,
y a la EoR. Como se ha mostrado, tests observacionales del progreso de la
reionización, como las anisotropías del CMB y los espectros de absorción de
quásares lejanos, podrían por tanto ser capaces de conferir cotas sobre los
diferentes modelos de materia oscura. Las restricciones sobre el espacio de
parámetros de la materia oscura obtenidas de este modo son, sin embargo,
altamente sensibles al modelado astrofísico del medio intergaláctico, y una
precisa determinación de los distintos procesos y fuentes astrofísicas es
necesaria para obtener límites robustos. Datos de la historia de ionización
como el CMB también pueden aplicarse a casos donde un modelo específico no es
asumido, estudiando por ejemplo procesos de reionización temprana, que podrían
ser inducidos por partículas de materia oscura que inyectasen energía en el
medio.

Especial atención merece la señal cosmológica de 21 cm, dada su gran
sensibilidad al estado térmico y de ionización del medio, siendo un observable
de gran interés. Empleando la localización temporal de la señal global de
absorción detectada por EDGES, es posible acotar diversos modelos de materia
oscura como templada e interactiva, mejorando algunas de las restricciones
actuales mediante otros métodos. Sin embargo, aún más prometedor resulta el
espectro de potencias de 21 cm, el cual sería capaz de arrojar información
sobre las inhomogeneidades del medio intergaláctico, además de ser más
fácilmente separable de los contaminantes astrofísicos, que pueden exceder en
varios órdenes de magnitud la señal cosmológica. Como ejemplo de su potencial,
se han calculado las posibles futuras cotas sobre la abundancia de agujeros
negros primordiales como materia oscura asumiendo la sensibilidad de detección
de la siguiente generación de radiointerferómetros, como SKA, mostrando que
dichas restricciones podrían mejorar en varios órdenes de magnitud las
presentes actualmente en el mismo rango de masas mediante otros métodos.

Además, se presentan aplicaciones de aprendizaje automático en cosmología, con
el fin de obtener el campo de densidad de materia subyacente y los parámetros
astrofísicos a partir de un mapa de 21 cm dado. Esto muestra, por un lado, el
gran potencial de la señal hiperfina del hidrógeno como herramienta para
estudiar el crecimiento de las inhomogeneidades y la formación y distribución
de estructuras. Y por otra parte, la gran capacidad y aplicaciones que pueden
tener los métodos basados en redes neuronales para extraer información
cosmológica, la cual puede ser inalcanzable de otro modo.


# Bayesian Analysis of Formula One Race Results: Disentangling Driver Skill
and Constructor Advantage

Erik-Jan van Kesteren & Tom Bergkamp

###### Abstract

Successful performance in Formula One is determined by combination of both the
driver’s skill and race-car constructor advantage. This makes key performance
questions in the sport difficult to answer. For example, who is the best
Formula One driver, which is the best constructor, and what is their relative
contribution to success? In this paper, we answer these questions based on
data from the hybrid era in Formula One (2014 - 2021 seasons). We present a
novel Bayesian multilevel rank-ordered logit regression method to model
individual race finishing positions. We show that our modelling approach
describes our data well, which allows for precise inferences about driver
skill and constructor advantage. We conclude that Hamilton and Verstappen are
the best drivers in the hybrid era, the top-three teams (Mercedes, Ferrari,
and Red Bull) clearly outperform other constructors, and approximately 88% of
the variance in race results is explained by the constructor. We argue that
this modelling approach may prove useful for sports beyond Formula One, as it
creates performance ratings for independent components contributing to
success.

##  1 Introduction

In most competitive sports with a large individual component (e.g., chess,
athletics, swimming, or tennis) success is determined primarily by the
relative skill (i.e., ability or talent) of the contestants. Official
competitions in such sports naturally result in official rankings representing
this skill level on an individual basis. Unlike these skill-based sports,
competitive motor racing has an additional key factor contributing to success:
the material, i.e., the race car or motor bike. In Formula One in particular,
the influence of the car on the results is considered to be substantial
(Budzinski and Feddersen,, 2020). In contrast to “spec” series, where cars
have the same specifications and are built by the same constructor, Formula
One cars are each built from the ground up by different constructors with
differing levels of technological and financial resources. These resource gaps
can lead to large differences in performance between cars, despite rules
imposed to counteract such performance differences. Arguably, the presence of
these large differences in constructor advantage have led to a single
constructor (Mercedes) dominating the sport in the “hybrid era”, from 2014 to
2020.

The dependence on materials reduces the relative influence of driver skill on
success — a race win is an entangled combination of both driver and
constructor performance. Therefore, ranking drivers in terms of skill level by
simply using competition race results is complex. Because of this problem,
perennial questions such as “who is the best Formula One driver?”, “which
constructor is the best?” and “is the driver or the constructor-team more
important to success?” are difficult to answer scientifically. These questions
are persistent in the sport; the 2016 world champion Niko Rosberg famously
posed that 80% of success in Formula One can be attributed to the car and 20%
to the driver (Bol,, 2020). In this article, we argue that it is now possible
to answer these questions due the widespread availability of race result data
and the accessibility of advanced statistical methodology. We propose a novel
Bayesian multilevel Beta regression model to answer three interrelated
questions for the hybrid era in Formula One: (a) what is the relative
influence of the driver and constructor on race results, (b) how do drivers
rank in terms of skill level, and (c) how do constructors rank in terms of
race car advantage?

Several attempts have been made to disentangle driver and constructor
performance. Eichenberger and Stadelmann, (2009) used linear regression with
dummy variables and several covariates to estimate driver and constructor-year
effects on race finishing position in the 1950-to-2006 Formula One seasons.
The driver-specific effects on race outcomes were then used to compute a
ranking of drivers’ skill level. Additionally, the authors found that
predictors for weather (wet vs. dry) and circuit type (street circuit vs.
permanent circuit) were relevant additions to the model. One shortcoming of
this study is the choice of outcome variable: because the number of
contestants per race changed over seasons (and sometimes even within seasons),
the interpretation of “finishing position” changed as well. This shortcoming
was addressed by Phillips, (2014), who analysed Formula One race data from
1950-to-2013. Here, the effects of driver performance, constructor performance
and season difficulty were estimated using an adjusted “points scored” outcome
variable, based on the official points scoring system used in the Formula One
between 1991 and 2002: 10 points for first place, 6 for second, down to 1
point for sixth place. With this approach, the authors ranked drivers in terms
of skill and concluded that Juan Manual Fangio was the best driver of all
time.

As Phillips, (2014) noted, there are two reasons that such models can
differentiate driver from constructor effects: (a) throughout the history of
Formula One, constructors have generally had two cars enter a race. Barring
minor differences in individual races, these cars have the same performance,
which allows for direct comparisons of a driver’s skill level against their
teammates. Furthermore, (b) drivers generally move to different constructor-
teams throughout their career, meaning their teammates also change. This
allows for simultaneous estimation of driver and constructor effects based on
race results.

A disadvantage of the aforementioned studies is that they used models with
many dummy variables (fixed effects) instead a multilevel (random effects)
model. A multilevel approach is beneficial, as it makes the models tractable
with fewer observations per driver and improves predictions for nested data
(Gelman,, 2006). A similar argument was made by Bell et al., (2016), who used
the same “points scored” outcome variable as Phillips, (2014), but used a
multilevel (random-coefficients) linear model to determine the driver and
constructor-year effects. Using this approach, the authors were the first to
directly estimate a parameter for comparing driver skill and constructor
advantage: they concluded that the constructor accounts for 86% of the
variance in points scored, and the driver for 14%.

While the approach presented by Bell et al., (2016) provides an answer to the
research questions posed above (pre-2014), we here present several
improvements which enable a more accurate insight into the current state of
Formula One. First and foremost, we argue that the “points scored” variable
leads to information loss: any result below sixth place leads to zero points,
making it impossible to differentiate constructors and drivers who
consistently finish below this threshold. Our approach, explained in section
3, is to create a Rank-Ordered Logit model for the rank ordering of finish
positions for each race. Second, we focus on the current hybrid era in Formula
One, from 2014 – 2021, using data from a publicly available resource (Newell,,
2021). This focus enables us to closely inspect changes across these most
recent seasons with comparable regulations. Third, we apply a Bayesian
workflow for model development (Gelman et al.,, 2013), visualisation (Gabry et
al.,, 2019), and comparison (Vehtari et al.,, 2017), which makes the
parameters interpretable and clarifies the model’s implications. Using our
approach, the parameters for driver skill and constructor advantage in our
model are directly interpretable as log-odds ratios of beating competitors,
similar to Elo ratings in chess (e.g., Van Der Maas and Wagenmakers,, 2005).
These parameters (and their uncertainty intervals), can then be used to create
rankings which provide clear insight into the relative performance of drivers
and constructors in Formula One.

The paper is structured as follows. First, in Section 2 we describe the data
source and processing steps we performed to obtain the predictors and the
outcome of interest. Then, in Section 3 we introduce the proposed Bayesian
multilevel regression model and its interpretation. In the same section, we
also perform model selection and model checking to validate the estimation
procedure. In Section 4 we perform inference for the 2014-to-2021 Formula One
seasons to answer the research questions surrounding driver skill and
constructor advantage. This includes a driver ranking for the 2021 season and
investigation of a counterfactual statement based on the estimated model:
would Hamilton in an Alfa Romeo beat Räikkönen in a Mercedes? Last, in Section
5 we place our contributions in context, discuss its implications, and provide
suggestions for future work. All analysis scripts and pre-processed data are
openly available in the supplementary material at
https://doi.org/10.5281/zenodo.7632045 (van Kesteren and Bergkamp,, 2023).

##  2 Data processing

We collected race results (driver id, constructor id, season (year), race
number, finishing position and status) for the 2014-to-2021 Formula One
seasons from the dataset behind the Ergast motorsports API (Newell,, 2021).
This resulted in a dataset of 160 races, with 51 unique drivers and 19
distinct constructors. In addition, we performed a data enrichment step by
scraping and parsing further race information from Wikipedia. In this step,
two predictors were added to the dataset which were previously found to be
relevant in the work by Eichenberger and Stadelmann, (2009) and Bell et al.,
(2016). We constructed a variable indicating whether the race was wet or dry,
and we also collected information about the circuit type (street circuit or
permanent circuit). As the information was not complete, especially for
weather type in the 2018 season, we manually completed this data using
publicly available race summaries. In total, in the period of interest, there
were 143 dry races and 17 wet races.

For the main analysis, non-finishers were removed from the analysis, removing
590 rows (from the original 3267 rows) in the dataset. Thus, we only examined
finished races for each driver: accidents and other reasons for non-finishing
or non-starting can be due to any number of reasons, which adds noise and
complexity to the outcome of interest. This removal of non-finishes has
implications for the interpretation of the model, which we discuss in section
3. In Appendix A we perform a sensitivity analysis for the effects of dealing
with finishing in different ways. Note that instead of the approach we take,
it is also possible to include non-finishing in the analysis by creating a
larger, joint model of non-finishes and race results conditional on finishing
(e.g., Ingram,, 2021).

A visual display of the finish positions for three drivers (Räikkönen,
Hamilton, and Giovinazzi) in the 2015-2020 seasons is shown in Figure 1. This
forms the basis of the outcome variable we construct, which is the per-race
ranking of competitors. By considering only race ranking as our outcome, we
disregard the starting position of the competitor, which is determined in a
qualifying session before the race. Thus, any conclusions about relative
performance of competitors necessarily includes both race performance and
qualifying performance: for all of our models, a race-winning performance is
no different if the competitor starts the race in first or last position. In
the next section, we explain in detail how we model this outcome variable.

Figure 1: Finish positions for Räikkönen, Hamilton, and Giovinazzi for the
seasons 2015 – 2020. Räikkönen’s move from Ferrari to Alfa Romeo in 2019 is
clearly visible in the finish positions. At Alfa Romeo, Räikkönen became
teammates with Giovinazzi, who he on average outperforms slightly.

##  3 Bayesian multilevel rank-ordered logit model

In this section, we describe in detail the process used to model the finish
position as a function of driver skill, yearly driver form, constructor
advantage, and yearly constructor form. All models in this paper were
estimated using the software package Stan (Carpenter et al.,, 2017) with
reasonable default priors for all parameter types (see section 2.1 of Bürkner,
(2017) for more details).

###  3.1 Basic model specification and parameter interpretation

For each race $r$, assume we have a set of competitors $\mathcal{C}_{r}$ which
compete for the win. The outcome variable in our model is a vector
representing a ranking of these competitors, $\boldsymbol{y}_{r}$. In our
generative model, we assume that this ranking follows a Rank-Ordered Logit
(ROL) distribution, based on a vector of latent abilities of the competitors
in that race $\boldsymbol{\vartheta}_{r}$:

$\displaystyle\boldsymbol{y}_{r}$ $\displaystyle\sim{\rm
RankOrderedLogit}(\boldsymbol{\vartheta}_{r})$ (1)

The precise definition of the ROL distribution is explained in detail in
Glickman and Hennessy, (2015), and in our extended model definition in
Appendix B. Here, note that each element of $\boldsymbol{\vartheta}_{r}$
represents the latent ability of each competitor in the race:

$\displaystyle\boldsymbol{\vartheta}_{r}$
$\displaystyle=\\{\vartheta_{c}\,|\,c\in\mathcal{C}_{r}\\}$ (2)

A competitor $c\in\mathcal{C}_{r}$ is defined as a pairing of a driver $d$ and
constructor (or team) $t$, in a particular season (or year) $s$. Thus, we
represent the skill of each competitor as the skill of the driver-constructor
pairing in a season: $\vartheta_{c}=\theta_{dts}$. In our model, this skill is
a sum of the average driver skill $\theta_{d}$, the driver’s seasonal form
$\theta_{ds}$, the constructor’s average advantage $\theta_{t}$ and the
constructor’s seasonal form $\theta_{ts}$. This results in the following
cross-classified multilevel model for each competitor’s latent ability:

$\displaystyle\vartheta_{c}=\theta_{dts}$
$\displaystyle=\theta_{d}+\theta_{ds}+\theta_{t}+\theta_{ts}$
$\displaystyle\theta_{d}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{d})$
$\displaystyle\theta_{ds}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{ds})$
$\displaystyle\theta_{t}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{t})$
$\displaystyle\theta_{ts}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{ts})$
(3)

This model form leads to a specific, natural interpretation for the parameters
$\theta$. The logit link function in the ROL likelihood, in combination with
the omission of an overall intercept, ensures that the (hypothetical) average
driver at an average constructor with an average seasonal form will on average
have $\theta_{dts}=0$, which translates into a probability of 0.5 of beating a
randomly selected other competitor. Then, $\theta_{d}$ represents the mean
driver skill as a log-odds ratio; e.g., if $\theta_{d}=0.3$, this means
(ceteris paribus) that the probability of finishing in front of the other
driver is $1/(1+e^{-0.3})\approx 0.57$. This parameter represents a deviation
from the average driver, so negative values mean worse than average skill, and
positive values mean better than average skill. We also include the seasonal
driver form parameter $\theta_{ds}$, which represents yearly deviations from
this long-term average driver skill.

A similar interpretation holds for $\theta_{t}$, which indicates the long-term
average constructor advantage. Constructors with positive values on this
parameter tend to produce cars which are better than average, and negative
values indicate cars which are worse on average. We also include the seasonal
constructor form parameter $\theta_{ts}$, which represents yearly deviations
from this long-term average team advantage. In Appendix D we compare this
parameterization to two other ways of dealing with time-varying $\theta$
parameters. For more detailed parameter interpretations and conclusions, see
Section 4.

Note that the latent skill are assumed to be stable within seasons. This means
that for races within the same season $s$ with the same competitors
$\mathcal{C}_{r}$ we assume independent and identically distributed rankings
$\boldsymbol{y}_{r}$. Note also that with this model formulation we implicitly
assume no correlation between the random intercepts for driver and car
constructor; there are no interactions at all between driver skill and team
advantage. This means that a driver’s skill is independent of the team
advantage, i.e., the driver skill does not change when the driver moves to a
different constructor.

###  3.2 Extending the basic model

Previous work has shown that several predictors may change the race results
(Bell et al.,, 2016). The first extension we make reflects the knowledge that
wet races are different from dry races. Wet races require a specific set of
skills, which rely less on the car and more on the driver. Like Bell et al.,
(2016), we represent this knowledge by splitting the driver average skill
parameter into a random intercept parameter $\gamma_{0d}$ and a random slope
parameter $\gamma_{1d}$ as in Equation 4:

$\theta_{d}=\gamma_{0d}+\gamma_{1d}\cdot\text{wet\\_race}$ (4)

where wet_race is an indicator (dummy) variable with a 1 if the race was wet
and 0 if the race was dry (see Section 2 for details). The driver average
skill in dry races is then $\gamma_{0d}$, and in wet races it is
$\gamma_{0d}+\gamma_{1d}$.

The second extension we make reflects the knowledge that different
constructors have different car philosophies. Theoretically, high-downforce
concept cars (e.g., Red Bull cars in the hybrid era) are relatively better
suited to narrow, curvy street circuits such as the famous Monaco circuit.
However, this advantage disappears on fast, permanent circuits such as Monza
(Italy). Therefore, we add a random slope to the constructor advantage
parameter, splitting it up as in Equation 5:

$\theta_{t}=\gamma_{0t}+\gamma_{1t}\cdot\text{permanent\\_circuit}$ (5)

where permanent_circuit is an indicator (dummy) variable with a 1 if the race
was on a permanent circuit and 0 if the race was on a street circuit (see
Section 2 for details).

With these two extensions, four models are possible: (a) the basic model, (b)
a weather model, (c) a circuit type model, and (d) a weather and circuit type
model. In the next subsection, we compare these models to select our final
model, on which we perform inference.

###  3.3 Model selection

We used efficient leave-one-out cross-validation (LOO, Vehtari et al.,, 2017)
to compare the four possible models. In short, LOO uses the samples from the
posterior to compute the expected log posterior density (ELPD) for each model,
which is an alternative to the standard information criteria in Bayesian model
comparison such as the Bayes Factor (marginal density) or DIC. For the four
tested models, the results are shown in Table 1.

| $ELPD$ | $SE_{ELPD}$ | $\Delta$ | $SE_{\Delta}$
---|---|---|---|---
Circuit | -3992.29 | 84.11 |  |
Basic | -3993.36 | 83.72 | -1.07 | 6.35
Circuit + Weather | -3993.66 | 84.31 | -1.37 | 1.94
Weather | -3995.55 | 84.09 | -3.25 | 6.22
Table 1: Model comparison results showing the expected log posterior density
(ELPD), its standard error, the difference between each model and the best
model, and the standard error of this difference. Note: a higher ELPD
indicates a better fit.

The circuit model is shown to be the best model in terms of ELPD, but the
differences are small compared to the standard errors of these differences.
This indicates that in terms of out-of-sample predictive performance, the
models perform similarly. To maintain parsimony in both modelling and
interpretation, we decide use the basic model for assessment, inference, and
prediction in the following sections.

###  3.4 Model assessment

Posterior samples were obtained via Hamiltonian Monte Carlo sampling with 8
chains of 1250 samples each after 1000 burn-in iterations. The effective
sample size for all parameters was higher than 2500, and their R-hat value was
smaller than 1.01, indicating adequate convergence. Trace plots for the main
parameters of the model are shown in Appendix C.

Posterior predictive checks (PPCs) are a vital part of the Bayesian workflow
(Gabry et al.,, 2019). In a visual PPC, simulated data $\tilde{y}$ from the
posterior predictive distribution is compared to the observed data $y$. If
$\tilde{y}$ approximates $y$ well, then the model captures the outcome well.
For our model, we performed two PPCs: one for the 2015 season (early hybrid
era) and one for the 2019 season (late hybrid era). The results are shown in
Figures 2 and 3.

Figure 2: Posterior predictive check for driver finishing positions in the
2015 season. Figure 3: Posterior predictive check for driver finishing
positions in the 2019 season.

The plots show satisfactory recovery of individual performances per driver-
constructor-season combination. It is noteworthy that for some drivers, the
observed data show bimodality (e.g., Vettel in 2019 or Bottas in 2015) while
the model cannot capture these patterns. In these cases, the averages in the
observed and simulated data still seem to align closely. In general, we
conclude that the model fits the observed data well. In the next section, we
use the model to perform parameter inference.

##  4 Results

In this section, we perform inference with the model that resulted from the
modelling procedure described in Section 3. To narrow down our inference
efforts, we focus on a subset of the drivers and teams competing in the 2021
season. The results for all drivers, teams, and seasons are available in the
supplementary material (van Kesteren and Bergkamp,, 2023).

###  4.1 Driver skill

After estimation of the preferred model, it is possible answer the question of
which driver is the most skilled, while taking into account constructor
advantage, constructor form, as well as all parameter uncertainties. In order
to produce a ranking, we obtained the posterior means and 89% credible
intervals (see McElreath,, 2018) of $\theta_{d}+\theta_{ds}$ for the season
2021. These summaries are shown in Figure 4.

Figure 4: Driver skill in 2021 on the log odds ratio scale, accounting for
yearly driver skill and constructor form. Error bars indicate 89% credible
interval.

In Figure 4 it is apparent that of the 2021 drivers, Hamilton and Verstappen
are ranked as the most skilled driver. Note that this ranking comes directly
from the model, which is estimated only on the hybrid-era data. Earlier
performances by drivers such as Vettel (four-time world champion in the period
2010–2013) and Räikkönen (world champion in 2007) have not been taken into
account, explaining their lower position on the ranking.

Because the model contains a yearly form parameter, we can also visualize the
latent skill trajectories of several drivers throughout the hybrid era, with
their credible intervals. The result of this visualization for 12 drivers of
the 2021 season is shown in Figure 5.

Figure 5: Driver skill trajectories for 12 drivers in the hybrid era on the
log odds ratio scale. Ribbons indicate 89% credible interval.

The figure shows that there are slight changes in skill across seasons.
Drivers such as Verstappen, Norris, and Gasly tend to improve over years, on
average, whereas Bottas displays a slight decline. Notably, Hamilton is
consistently at the top, whereas Alonso and Sainz are consistently slightly
above average in terms of their latent skill.

###  4.2 Constructor advantage

For the constructors competing in 2021, we here investigate how much of an
advantage their car yields. We do this by computing the posterior means and
89% credible intervals of $\theta_{t}$ from the model. The result is shown in
Figure 6.

Figure 6: Long-run average constructor advantage on the log odds scale. Error
bars indicate 89% credible interval.

One feature that becomes clear from this constructor advantage plot is the
relative advantage of the _big three_ teams: Mercedes, Ferrari, and Red Bull
have the largest budget and the most resources to spend on developing their
car, which has resulted in these three teams excelling in the hybrid era.
Another interesting feature in this plot is the large uncertainty around the
teams that have competed in only a few seasons. For example, Alpine and Aston
Martin were new teams in 2021, and therefore have not had a chance compete in
many races, resulting in uncertainty around where it is placed in this
constructor ranking. Again, these parameters need to be interpreted with care:
they represent the average constructor advantage over the entire hybrid era.

The last random intercept component is the constructor-year effects
$\theta_{ts}$. These represent yearly constructor form, as a deviation from
their long-term average advantage. In Figure 7, the yearly constructor
advantage trajectories for a selection of teams is shown from 2014 to 2021.

Figure 7: Yearly constructor advantage (summing constructor and constructor-
year random effects) for Ferrari, Mercedes, Red Bull, McLaren, Haas, and
Williams in the period 2014 – 2020. Ribbons indicate 89% credible interval.

One of the most striking features from this graph is Ferrari’s drop in form
from 2019 to 2020. There is a good explanation for this: after the 2019
season, there were allegations that Ferrari’s engine in the previous year did
not match league regulations. Ferrari subsequently reached a settlement with
the Formula One Management (Formula1.com,, 2020), which left the team with a
relatively weak engine in the 2020 season. This change is aptly reflected in
the constructor form parameter for Ferrari, but also to a lesser extent in
that of McLaren, which profited from Ferrari’s drop and ended up third in the
constructor’s championship that year.

###  4.3 Season performance

By combining driver skill and constructor advantage for each entrant in the
2021 season, we can create a posterior predictive simulation for the entire
season. Here, we simulate finish positions (for the races where the entrant
actually finished), which we then translate to points using a slightly
simplified version of the 2021 points system (i.e., excluding fastest lap
points and sprint race points). From this, we compute the expected average
points per race over the whole 2021 season, as well as its 89% credible
interval. Then, we compare this to the realized average points for each
entrant using the same simplified points system. The result of this is shown
in Table 2.

Driver | Constructor | Expected points [89% CI] | Observed points
---|---|---|---
Verstappen | Red Bull | 15.30 [12.32, 18.00] | 18.00
Hamilton | Mercedes | 14.58 [11.45, 17.50] | 17.59
Bottas | Mercedes | 8.85 [5.95, 11.73] | 9.82
Sainz | Ferrari | 8.79 [5.73, 11.91] | 7.45
Norris | McLaren | 7.07 [4.23, 10.18] | 7.23
Leclerc | Ferrari | 7.00 [4.18, 9.95] | 7.32
Perez | Red Bull | 6.75 [4.00, 9.73] | 8.59
Ricciardo | McLaren | 5.71 [3.09, 8.64] | 5.41
Gasly | Alpha Tauri | 5.40 [2.86, 8.14] | 5.14
Alonso | Alpine | 4.13 [1.95, 6.73] | 3.68
Ocon | Alpine | 3.95 [1.82, 6.50] | 3.50
Stroll | Aston Martin | 2.61 [0.82, 4.86] | 1.55
Tsunoda | Alpha Tauri | 2.30 [0.68, 4.32] | 1.45
Vettel | Aston Martin | 2.08 [0.55, 4.09] | 2.18
Raikkonen | Alfa Romeo | 1.93 [0.45, 3.82] | 0.45
Russell | Williams | 1.26 [0.14, 2.86] | 1.14
Giovinazzi | Alfa Romeo | 1.95 [0.45, 3.86] | 0.14
Latifi | Williams | 0.80 [0.00, 2.14] | 0.36
Schumacher | Haas | 0.28 [0.00, 1.14] | 0.00
Mazepin | Haas | 0.12 [0.00, 0.73] | 0.00
Table 2: For entrants in the 2021 season, the model’s posterior expectations
and 89% credible interval for expected points per race (averaged over the 2021
season), and the realized 2021 average points per race (excluding points for
fastest laps). Cases where the observed average points are outside of the 89%
CI are underlined in the last column.

Generally, the expected and observed columns line up well in terms of credible
interval coverage, although the model tends to underestimate points at the top
and overestimate points at the bottom. This may be due to the regularizing
effects of the multilevel model implementation, both on average performance
and seasonal form.

###  4.4 Relative contributions of drivers and constructors

In order to investigate the contributions of drivers and constructors to the
race results, we investigate the standard deviations of the random intercepts
in the model. By investigating these standard deviations, we can make
conclusions about which matters more in terms of race results: the driver or
the car. The posteriors for the variation coefficients are shown in Table 3.

Component | Symbol | Estimate | Est.Error | Lower | Upper
---|---|---|---|---|---
Constructor advantage | $\sigma_{c}$ | 1.63 | 0.34 | 1.14 | 2.27
Constructor form | $\sigma_{cs}$ | 0.73 | 0.10 | 0.57 | 0.91
Driver skill | $\sigma_{d}$ | 0.54 | 0.12 | 0.35 | 0.76
Driver form | $\sigma_{ds}$ | 0.35 | 0.06 | 0.24 | 0.46
Table 3: Standard deviations ($\sigma$) for the random effects in the model.
Lower and upper represent bounds of the 89% credible intervals.

The standard deviation of the constructor is larger than that for the driver.
This means that on average, the constructor has a larger impact on race
results than the driver. Rephrasing this, on average the correlation in the
outcome is stronger for two different drivers driving for the same team than
for the same driver driving for different teams. This interpretation is also
exemplified by the race results shown in Figure 1: Räikkönen driving for Alfa
Romeo (2020) looks more like his teammate (Giovinazzi) than like Räikkönen
driving for Ferrari (2018).

Quantifying the relative importance of long-term constructor advantage
compared to driver skill is also possible directly from the numerical
summaries. The posterior estimates for the variances are as follows:
$\sigma^{2}_{c}\approx 2.65$, $\sigma^{2}_{cs}\approx 0.54$,
$\sigma^{2}_{d}\approx 0.29$, and $\sigma^{2}_{ds}\approx 0.12$. Following the
methodology of Bell et al., (2016, §4.2), this means that constructor effects
account for around 88% of the variance in the model (89% CI [0.775, 0.945]),
which is very similar to the aforementioned authors who reported 86%.

###  4.5 Counterfactual inference

Using samples from the posterior distributions of the parameters, we can
answer some counterfactual questions about the drivers in the model. An
example question would be: “According to the model, would Hamilton be expected
to beat Räikkönen in a race in 2021 if Hamilton drove for Alfa Romeo and
Räikkonen for Mercedes?”.

We can answer this question by computing the posterior probability of Hamilton
beating Räikkönen $\pi_{\rm ham>rai}$:

$\displaystyle\theta_{\rm ham:alfa:2021}$ $\displaystyle=\theta_{\rm
ham}+\theta_{\rm ham:2021}+\theta_{\rm alfa}+\theta_{\rm alfa:2021}$
$\displaystyle\theta_{\rm rai:merc:2021}$ $\displaystyle=\theta_{\rm
rai}+\theta_{\rm rai:2021}+\theta_{\rm merc}+\theta_{\rm merc:2021}$
$\displaystyle\pi_{\rm ham>rai}$ $\displaystyle=\frac{\exp(\theta_{\rm
ham:alfa:2021})}{\exp(\theta_{\rm ham:alfa:2021})+\exp(\theta_{\rm
rai:merc:2021})}$ (6)

The posterior of $\pi_{\rm ham>rai}$ is shown in Figure 8. It is shown that
Räikkönen is expected to beat Hamilton in this scenario ($E[\pi_{\rm
ham>rai}]<0.5$), with a clear degree of uncertainty. Note that this
counterfactual prediction is a way to summarise the model, meaning the same
assumptions that accompany the model also accompany the counterfactual
predictions. For example, for these predictions the data from before 2014 is
irrelevant, and driver talent is independent of the constructor advantage.

Figure 8: Posterior distribution of the win probability of Hamilton in an Alfa
Romeo over Räikkönen in a Mercedes in a race during the 2021 season. The
expected value of the distribution (vertical line) is 0.36, meaning that
Räikkönen is expected to beat Hamilton in this scenario.

##  5 Discussion

In this paper, we used a Bayesian multilevel rank-ordered logit model for the
rank results in Formula 1 races of the hybrid era (2014 – 2021). After model
development, comparison with leave-one-out cross-validation, and validation
with posterior predictive checks, we have made inferences about the drivers
and constructors competing in the 2021 season. In terms of the drivers,
Hamilton and Verstappen are the best drivers. In terms of the constructors,
the top three teams (Mercedes, Ferrari, and Red Bull) clearly outperform the
rest on average in the hybrid era. Additionally, the model accurately
represents changes in constructors’ seasonal form, for example reproducing a
drop in performance by Ferrari in 2020. Comparing driver contributions to team
contributions, we have concluded that the car is more important than the
driver when it comes to race results. Using our model and posterior sampling,
it is possible to obtain answers of counterfactual questions as posterior
distributions, which we have shown in the last part of Section 4.

In terms of parameter interpretation, there is an interesting parallel between
this model and Elo ratings (Elo,, 1978) in chess. Both the Elo rating and our
driver talent parameters can be transformed using an inverse logit to compare
the relative strength of the competitors, and thus how likely it is that one
competitor wins. The larger the difference between the ratings, the more
certain it is that the competitor with the higher rating wins. We have shown
such a comparison in a counterfactual situation in Section 4. Note that this
Bayesian hierarchical approach has been applied before to different sports
(e.g., in tennis; Ingram,, 2019). However, in our model not only the athletes
get ratings, but also the constructors, so comparisons can be made at this
level as well. This approach could be used in other sports where multiple
independent components contribute to competition results.

While this model does well at describing past data, for example closely
reproducing the ranking of the 2021 season, it is probably not suitable for
prediction. It uses very limited information (only the driver, constructor,
and year) and for each year a specific effect needs to be estimated. Before a
season starts, there is no data on that season, meaning that these year-
effects are unavailable (even though they are important components of this
model). For forecasting, approaches such as that of Henderson et al., (2018)
may be more suitable relative to the baseline of bookmaker odds.

There are several areas where this model may be improved. One area is in team
continuity: teams can officially change their name, when behind the scenes it
is the same team, with the same long-run performance. For example, Alpha Tauri
is a re-branding of the Italian Toro Rosso team, but it enters our dataset as
a completely new team, with understandably large credible intervals around
performance in 2020. By not accounting for team continuity across different
team names, we had in total 17 different constructors in the dataset. On the
other hand, team name changes often do go hand-in-hand with some structural
changes, and it is hard to determine the extent to which this happens: where
to draw the line between a “re-brand” and a new team?

While our model accurately represents driver and constructor performances in
the seasons 2014 – 2020, the data range could be expanded in order to provide
results that better reflect the careers of certain drivers. For example,
Räikkönen – who is at the end of his career – is in the bottom half in Figure
4 and Table 2. Because the biggest part of Räikkönen’s career is missing from
the data, his 2007 world championship has not been appropriately taken into
account. These errors propagate, perhaps overestimating the performance of
Alfa Romeo (Räikkönen’s team for that year) and underestimating the skill of
Giovinazzi (whose only teammate has been Räikkonen).

## References

  * Bell et al., (2016) Bell, A., Smith, J., Sabel, C. E., and Jones, K. (2016).  Formula for success: multilevel modelling of formula one driver and constructor performance, 1950–2014.  Journal of Quantitative Analysis in Sports, 12(2):99–112.
  * Bol, (2020) Bol, R. (2020).  How to win in formula one: is it the driver or the car?  The Correspondent.
  * Budzinski and Feddersen, (2020) Budzinski, O. and Feddersen, A. (2020).  Measuring competitive balance in formula one racing.  In Outcome Uncertainty in Sporting Events. Edward Elgar Publishing.
  * Bürkner, (2017) Bürkner, P.-C. (2017).  brms: An R package for bayesian multilevel models using stan.  Journal of statistical software, 80(1):1–28.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of statistical software, 76(1).
  * Eichenberger and Stadelmann, (2009) Eichenberger, R. and Stadelmann, D. (2009).  Who is the best formula 1 driver? an economic approach to evaluating talent.  Economic Analysis & Policy, 39(3).
  * Elo, (1978) Elo, A. (1978).  The rating of chess players, past and present.  Arco, New York.
  * Formula1.com, (2020) Formula1.com (2020).  Fia reaches ‘settlement’ with Ferrari following 2019 engine investigation.  https://www.formula1.com/en/latest/article.fia-reaches-settlement-with-ferrari-following-2019-engine-investigation.6beur1atKeTLvJHPEuHUJW.html (accessed May 2021).
  * Gabry et al., (2019) Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., and Gelman, A. (2019).  Visualization in bayesian workflow.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 182(2):389–402.
  * Gelman, (2006) Gelman, A. (2006).  Multilevel (hierarchical) modeling: what it can and cannot do.  Technometrics, 48(3):432–435.
  * Gelman et al., (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013).  Bayesian data analysis.  CRC press.
  * Glickman and Hennessy, (2015) Glickman, M. E. and Hennessy, J. (2015).  A stochastic rank ordered logit model for rating multi-competitor games and sports.  Journal of Quantitative Analysis in Sports, 11(3):131–144.
  * Hamaker and Grasman, (2015) Hamaker, E. L. and Grasman, R. P. (2015).  To center or not to center? investigating inertia with a multilevel autoregressive model.  Frontiers in psychology, 5:1492.
  * Henderson et al., (2018) Henderson, D. A., Kirrane, L. J., et al. (2018).  A comparison of truncated and time-weighted plackett–luce models for probabilistic forecasting of formula one results.  Bayesian Analysis, 13(2):335–358.
  * Ingram, (2019) Ingram, M. (2019).  A point-based bayesian hierarchical model to predict the outcome of tennis matches.  Journal of Quantitative Analysis in Sports, 15(4):313–325.
  * Ingram, (2021) Ingram, M. (2021).  A first model to rate formula 1 drivers.  https://martiningram.github.io/f1-model/ (accessed March 2022).
  * McElreath, (2018) McElreath, R. (2018).  Statistical rethinking: A Bayesian course with examples in R and Stan.  Chapman and Hall/CRC.
  * Newell, (2021) Newell, C. (2021).  Ergast developer API.  http://ergast.com/ (accessed February 2022).
  * Phillips, (2014) Phillips, A. J. (2014).  Uncovering formula one driver performances from 1950 to 2013 by adjusting for team and competition effects.  Journal of Quantitative Analysis in Sports, 10(2):261–278.
  * Van Der Maas and Wagenmakers, (2005) Van Der Maas, H. L. and Wagenmakers, E.-J. (2005).  A psychometric analysis of chess expertise.  The American journal of psychology, pages 29–60.
  * van Kesteren and Bergkamp, (2023) van Kesteren, E.-J. and Bergkamp, T. (2023).  vankesteren/f1model: Rank-ordered logit model.
  * Vehtari et al., (2017) Vehtari, A., Gelman, A., and Gabry, J. (2017).  Practical bayesian model evaluation using leave-one-out cross-validation and waic.  Statistics and computing, 27(5):1413–1432.

##  Appendix A Sensitivity analysis for race finishing

In the main paper, we remove non-finishers from the data, estimating the model
parameters data only from the races in which a certain competitor finished the
race. This choice influences the interpretation of the parameters: the skill
parameters exclude any measure of ”reliability” because only finished races
count. In this section, we perform sensitivity analysis for this choice,
investigating the effects on the results of making a different choice here:

  1. 1.

Not removing any results, i.e., ranking unfinished competitors based on race
distance.

  2. 2.

Removing only car-related non-finishes and ranking the driver-related non-
finishes as normal

  3. 3.

Removing only driver-related non-finishes and ranking the car-related non-
finishes as normal

  4. 4.

Only ranking finished competitors (i.e., original analysis from main paper)

As mentioned in the main text in section 2, yet another option is to include
non-finishes in the model directly, resulting in a joint binomial model for
finishing and rank model for ranking conditional on finishing. This approach
is outside the scope of this sensitivity analysis.

In the data, it is difficult to accurately ascertain which results belong to
car-related and driver-related errors: there are 51 different status types
other than finishing the race; some are ambiguous with respect to who caused
them, e.g., “Excluded”, “Technical”, “Damage”. For the purpose of this
sensitivity analysis, we have made the following selection:

  * Driver-related:

Collision, Disqualified, Withdrew, Retired, Accident, Collision damage, Spun
off, Excluded, Illness

  * Car-related:

ERS, Oil pressure, Engine, Technical, Gearbox, Electrical, Power Unit, Brakes,
Clutch, Exhaust, Mechanical, Turbo, Rear wing, Drivetrain, Suspension, Oil
leak, Water leak, Water pressure, Electronics, Transmission, Wheel, Power
loss, Fuel system, Front wing, Tyre, Throttle, Brake duct, Hydraulics,
Battery, Puncture, Overheating, Wheel nut, Vibrations, Driveshaft, Fuel
pressure, Seat, Spark plugs, Steering, Damage, Out of fuel, Debris, Radiator

The results of the sensitivity analysis are shown graphically in Figures 9, 10
(driver skills) and 11 (constructor advantage). A general pattern in comparing
the main paper model (”only finishers retained”) to the other three models is
the smaller variance of these random effect components; when we include non-
finishes in the ranking data, this ranking becomes more noisy so the latent
skill and advantage parameters will be closer together (and closer to 0). This
is especially visible in Figure 11, which was faceted by the model variable
for this purpose.

For the drivers, generally the skill patterns over time (Figure 9) are similar
across different models, meaning that skill increases and decreases are
captured somewhat similarly for the different models. Looking at the skill
level in 2021, we see that a notable outlier is Bottas, who drops dramatically
in skill level relative to his surrounding drivers once non-finishes are
included in the model. It is unclear what the mechanism is behind this drop,
as Bottas has not had much more non-finishes than other drivers (4 per season,
whereas the average is 3).

In summary, the results are somewhat sensitive to the choice of what to
include in the model. The conclusions in the main paper should be interpreted
carefully; we once again make clear that the interpretation of the skill
parameters only includes _finished_ races, thus do not contain any component
of reliability of the car or carefulness of the driver. For example, Pastor
Maldonado — who was notoriously crash-prone — is ranked $6^{th}$ worst if we
include all data in the model, but a much better $19^{th}$ of 38 drivers if we
only look at finished races.

Figure 9: Sensitivity analysis for driver skills inferences. Figure 10:
Sensitivity analysis for driver skills inferences for the 2021 season. Figure
11: Sensitivity analysis for constructor advantage inferences.

##  Appendix B Extended model definition

In this section, we extend the model definition in Section 3 in the main paper
into a fully specified generative model. Following Glickman and Hennessy,
(2015), the Rank-Ordered Logit (ROL) distribution implies that the skills for
each competitor in the $r^{th}$ race ($\boldsymbol{\vartheta}_{r}$) generate a
latent _performance_ vector $\boldsymbol{z}_{r}$ following independent
standard extreme value (Gumbel) distributions. The resulting rank-ordering
$\boldsymbol{y}_{r}$ is the rank-ordering of these latent performance values.
By replacing RankOrderedLogit by this generative process, the full model can
be specified as follows:

$\displaystyle\boldsymbol{y}_{r}$ $\displaystyle={\rm
rank}(\boldsymbol{z}_{r})$ $\displaystyle\boldsymbol{z}_{r}$
$\displaystyle=\\{z_{c}\,|\,c\in\mathcal{C}_{r}\\}$ $\displaystyle z_{c}$
$\displaystyle\sim{\rm Gumbel}(\vartheta_{c})$
$\displaystyle\vartheta_{c}=\theta_{dts}$
$\displaystyle=\theta_{d}+\theta_{ds}+\theta_{t}+\theta_{ts}$
$\displaystyle\theta_{d}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{d})$
$\displaystyle\theta_{ds}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{ds})$
$\displaystyle\theta_{t}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{t})$
$\displaystyle\theta_{ts}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{ts})$
(7)

where we make use of the standard Gumbel distribution with only a location
parameter $\mu$ and no scale:

$f(x;\mu)=\exp\left(x-\mu-\exp(x-\mu)\right)$ (8)

The likelihood for the model in equation 7 can be simplified for
implementation (Glickman and Hennessy,, 2015, eq. 4). Assuming, without loss
of generality, that the set of competitors for a specific race
$\mathcal{C}_{r}$ is already ordered according to the race result
$\boldsymbol{y}_{r}$, we obtain the following likelihood:

$p(\boldsymbol{y}_{r}|\boldsymbol{\vartheta}_{r})=p(z_{1}>z_{2}>...>z_{m_{r}}|\boldsymbol{\vartheta}_{r})=\prod_{i=1}^{m_{r}-1}\frac{\exp{(\vartheta_{i})}}{\sum_{j=i}^{m_{r}}\exp{(\vartheta_{j})}}$
(9)

where $m_{r}$ is the number of competitors $|\mathcal{C}_{r}|$.

##  Appendix C Traceplot

Figure 12 shows the traceplot for the main standard deviation parameters of
the rank-ordered logit model. From top to bottom, the standard deviations for
the four random effects $\beta_{d}$, $\beta_{ds}$, $\beta_{c}$, and
$\beta_{cs}$ are shown. All traceplots for the four different chains overlap
as expected from an appropriately converged model.

Figure 12: Monte carlo markov chain visualisation for the main model
parameters, indicating satisfactory convergence for all shown parameters.

##  Appendix D Comparing different dynamic skill parameter implementations

In the paper, the model does not take time into account explicitly, but it
estimates i.i.d. seasonal form parameters per driver and constructor. Here,
this approach is briefly compared to explicitly using a latent multilevel
AR(1) implementation (auto), and to a more parsimonious multilevel intercept +
slope model (slope) for the latent skill parameters. For the auto-regressive
model, we used a cluster mean centered parameterization as in Hamaker and
Grasman, (2015).

Because these models have the same likelihood form (the rank-ordered logit
likelihood), we can compare them using efficient leave-one-out cross-
validation:

| $ELPD$ | $SE_{ELPD}$ | $\Delta$ | $SE_{\Delta}$
---|---|---|---|---
Auto | -3992.7 | 83.9 |  |
Rank | -3993.8 | 83.9 | -1.1 | 1.7
Slope | -4042.9 | 81.6 | -50.2 | 18.4

The simple slope model is quite clearly worse than the other two
implementations; it is not able to capture the skills properly. For example,
looking at Giovinazzi and Räikkönen in the plot below shows that they are
estimated to be top drivers in the 2021 season, which is unlikely given their
results.

Figure 13: Posterior distribution of driver skill parameters over seasons for
the slope model. Ribbons indicate 89% credible interval.

The AR(1) and default rank model are very similar in terms of LOOCV
performance. Looking at the estimated skill trajectories, it becomes clear why
this is the case:

Figure 14: Posterior distribution of driver skill parameters over seasons for
the multilevel and AR(1) implementations. Ribbons indicate 89% credible
interval.

The driver skills have very similar estimates. The same holds for the
constructor advantage over time:

Figure 15: Posterior distribution of team advantage parameters over seasons
for the multilevel and AR(1) implementations. Ribbons indicate 89% credible
interval.

Based on this comparison, the decision was made to work with the more
interpretable basic multilevel model in the main paper.