Influence of advanced footwear technology on sub-2 hour marathon and other top running performances

A.1 Extremes for identically distributed variables

Extreme value theory (EVT) is a branch of statistics which studies the tails of probability distributions. It was first developed for block maxima (Gumbel 1958) analysis, but the Peaks Over Threshold (POT) method (Davison and Smith 1990) is often preferred, as it uses all the most extreme data, rather than just the maxima, typically leading to more efficient inference. Let X be a random variable with distribution function F, if there exist random sequences a _n , b _n > 0 such that

as n → ∞ is a non-degenerate limiting distribution, then for a large enough threshold u we can use the approximation

where σ _u = σ + ξ(u − μ) > 0, a₊ = max(a, 0). If ξ < 0 then x must lie in the interval [0, x _H ], where x _H = u − σ _u /ξ is the upper limit of the distribution, whereas if ξ ≥ 0, x can take any positive value. The limit distribution H _u , called Generalized Pareto distribution (GPD) motivates an approximation for large u, giving a model for the distribution of the exceedances above such threshold, regardless of the distribution F.

Given a large enough sample of n independent identically distributed (IID) observations, in the POT approach a threshold u is carefully chosen, and exceedances can be used to estimate the parameters of the GPD. Threshold choice can be rather subjective and case-dependent, and is subject to a bias-variance trade-off. In this paper we base our choice on graphical diagnostics; however, other alternative methods might also be suitable; see Scarrott and MacDonald (2012) for a detailed review of these techniques.

It is remarkable that the rate of the frequency of exceedances above the threshold u can be derived in a fashion that gives way to a more complete perspective of exceedances modelling, using point process models. Let X _i be IID random variables with distribution function F, we define

where 1 ( A ) is an indicator whether the event A occurs. It follows that N _n (x) ∼ Binomial(n, 1 − F(a _n x + b _n )) with mean n 1 − F ( a n x + b n ) , and using the classical Poisson limit of the binomial distribution,

where λ = { 1 + ξ ( x − μ ) / σ } + − 1 / ξ .

Therefore we can construct a model for extreme tails with two components: a model for the number of exceedances, given by (4), which is Poisson distributed with mean λ = { 1 + ξ ( x − μ ) / σ } + − 1 / ξ , and a model for the distribution of the exceedances, which is GPD distributed, following H _u (x).

Consider the sequence of point processes on R 2 (Coles 2001)

where the scaling 1/(n + 1) in the first coordinate ensures that the time axis is continuous on (0, 1), and the sequences a _n , b _n are defined in (1). More precisely, on regions of the form [0, 1] × (u, ∞), where u is large enough such that (2) approximately holds, we have that P _n → P as n → ∞, where P is a non-homogeneous Poisson Process. Consequently, the integrated measure Λ of P on A 1 , u = [ 0,1 ] × ( u , ∞ ) is given by

and its intensity function is

with x > u and 0 < t ≤ 1. For statistical inference we assume that for large enough n, P _n ∼ P is a good approximation. The scaling coefficients a _n , b _n , can be absorbed into the intensity function, so we work directly with the series i n + 1 , X i : i = 1 , … , n . Therefore, for a region of the form A 1 , u = [ 0,1 ] × ( u , ∞ ) , containing n points x = ( t 1 , x 1 ) , … , ( t n , x n ) , the likelihood for the parameters θ = (μ, σ, ξ) is

A.2 Extremes of non-stationary sequences

The extreme value models derived so far are built on the assumption of IID variables. However, in our work, non-stationarity data arise due to the improvement of racing conditions over time, and the potential Vaporfly shoes. Therefore, we relax the identically distributed assumption by introducing a time-dependent structure, while keeping independence assumption. Indeed, the time variation for parameters θ(t) = {μ(t), σ(t), ξ(t)} will translate into a time-dependent rate of exceedances, and distribution of such exceedances. Under this covariate structure, the intensity of the non-homogeneous Poisson process P will be

Now, in the general case where we have n points x = ( t 1 , x 1 ) , … , ( t n , x n ) in the region A T , u = [ 0 , T ] × ( u , ∞ ) , the integrated intensity becomes

and the full likelihood is

The parameters θ ( t ) = ( μ ( t ) , σ ( t ) , ξ ( t ) ) are estimated by maximizing (11), and with such estimates, for a given time t, predictions about the number of exceedances can be made by integrating (9). The excess distribution at time t will be given by

where σ _u (t) = σ(t) + ξ(t){u − μ(t)}.

A.3 Model

For most disciplines (and specially for marathon-men) a linear dependence on time for the scale parameter of the GP distribution of the exceedances was best suited in AIC terms. The following parametrisation was used to incorporate such structural time dependence.

where d ∈ D is the superscript denoting discipline d, y(t) is the year corresponding to time t, ξ ( d ) , μ 0 ( d ) ∈ R , σ 0 ( d ) ∈ R + are the shape, location, and scale parameter of the Poisson process, β ∈ R controls the linear trend in σ^(d)(t) and μ^(d)(t), γ ( d ) ∈ R represents Vaporfly shoes effect, 1 is the indicator function, and 2018 is the year when the shoes started to be widely used in official races. Note that this parametrisation enforces the GPD scale parameter for exceedances above u _d to change linearly with time.

A.4 Expected running times of next new world record

As derived in Spearing et al. (2021), the expected new world record time for discipline d at year y will be

where σ r d ( d ) ( y ) = σ 0 ( d ) + ξ r d − μ 0 ( d ) + δ y ( t ) , X y * ( d ) is the random variable denoting the running-time of a new world record for discipline e, set in year y, and r _d is the current (2019) world record of discipline d, so that r _d = max(X ^(d) ), with X ^(d) the set of all observations for discipline d.

A.5 Probability of breaking a world record in a given year

Let N y ( d ) be the number of exceedances of the threshold u _d for discipline d during year y, it is Poisson distributed with mean

Therefore, let X 1 : N y ( d ) ( d ) = X i ( d ) , i = 1 , … , N y ( d ) , where X i ( d ) ∼ i i d H u ( d ) ( y ) , if we denote by Pr ( R y ( d ) ) the probability that a world record for discipline d is set in year y,

where H ̄ u ( d ) ( r d , y ) ≔ 1 − H u ( d ) ( r d , y ) .

A.6 Time until next world record is set

Let T^(d) be the random variable describing the waiting time until a new world record is set for an discipline e, if we define t _y = y − 2020, the probability F T ( d ) ( t y ) = Pr ( T ( d ) < t y ) that a world record for discipline e is set before some year y is

We can further estimate the expected waiting time until the world record is broken for any discipline e, which has the following expression

where Pr(R _y ) is the probability that the world record is broken at year y, as described in (19).

A.7 Adjusting for AFT effect

Let x > u be a running-time recorded during year y > 2018, when Vaporfly and other shoes with technological advances are widely used in official races. We denote by x _c the corrected or equivalent time of x if such shoes were not used. Its expression can be derived as in Spearing et al. (2021), obtaining

where Λ C ( d ) ( A y , u ) has the form of Λ ( d ) ( A y , u ) but with the corrected parameters

A.8 Breaking the 2 h marathon

Let Pr(B2 = y) be the probability the 2-h marathon being broken in a given year y, it follows from (19) that

where 2h ≔ −7200 and marM refers to the marathon-men discipline. Additionally, we could compute the cumulative probability of achieving a sub 2-h marathon before year y, which follows from (20)