Modelling the age distribution of longevity leaders
Abstract
Human longevity leaders with remarkably long lifespan play a crucial role in the advancement of longevity research. In this paper, we propose a stochastic model to describe the evolution of the age of the oldest person in the world by a Markov process, in which we assume that the births of the individuals follow a Poisson process with increasing intensity, lifespans of individuals are independent and can be characterized by a gamma–Gompertz distribution with time-dependent parameters. We utilize a dataset of the world’s oldest person title holders since 1955, and we compute the maximum likelihood estimate for the parameters iteratively by numerical integration. Based on our preliminary estimates, the model provides a good fit to the data and shows that the age of the oldest person alive increases over time in the future. The estimated parameters enable us to describe the distribution of the age of the record holder process at a future time point.
keywords:
human longevity, world’s oldest person, stochastic model, gamma–Gompertz lifespan distributionIntroduction
Exceptionally long human lifespans are one of the cornerstones of demography and mortality research. Studying the group of record holders may reveal not only the underlying mortality mechanism of a population but also potentially shed some light on the future developments of human longevity. As life expectancy increases [1] with deaths shifting to older ages, the distribution of deaths at the oldest-old ages [2, 3] gains more interest of demographers, actuaries and decision makers of numerous disciplines.
The pattern of adult human mortality has already been described [4] but there is still a debate about the exact distribution of deaths at adult ages. More details on the necessity of a correctly specified model for the underlying mortality process and its impact on further research are discussed in [5, 6]. Numerous publications review the ongoing discussion on the existence of a mortality plateau [7, 8, 9, 10, 11, 12], and the levelling-off of adult human death rates at the oldest ages is supported by the findings in [13, 14, 15] while others cast some doubt on this observation [16, 17, 18]. If there is a mortality plateau then the distribution of deaths at the oldest-old ages must be gamma-Gompertz and human lifespan can increase further without any maximum [19]. Further studies discuss the existence of a limit to human lifespan with more focus on the extreme value distribution aspect of the deaths at the oldest-old ages for various populations and different lifetime distributions [20, 21, 22, 23, 24]. These models can be helpful in determining the plausibility of longevity leaders as well.
We contribute to this discussion by proposing a stochastic model to describe the evolution of the age of the world’s oldest person. Based on our estimates the model provides a good fit to the titleholder data since 1955, collected by the Gerontology Research Group [25]. With the model results, it is possible to predict the age of the oldest person in the world in the future. When should we expect to see the next Jeanne Calment, the supercentenarian with the longest human lifespan ever documented? Will her record ever be surpassed? Our results provide a prediction for the age distribution of the record holder in the coming decades to answer these questions.
Results
Our model describes the evolution of the age of the oldest living person under the following assumptions. We assume that the births of individuals follow a Poisson process with time-dependent intensity [26]. The lifespans of individuals in the population are independent and their distribution may depend on the date of birth. Then the age of the record holder in the population evolves in time as a Markov process with explicit transition probabilities. As the first main result of this paper, we explicitly compute the distribution of the age of the record holder for any given birth rate parameter and lifespan distribution. The detailed mathematical description and the properties of the general model are described in Section 1.
We apply our general result to the case which approximates the human birth rate over the world and the human lifespan. We specify the intensity function of the Poisson process of births to have an exponential growth in time. The underlying force of mortality is chosen so that it follows an extension of the Gompertz mortality model [4] and the lifespan distribution of individuals is given by the gamma–Gompertz () distribution with time-dependent parameters. This distribution adequately captures the slowing down of senescence mortality at the oldest old ages. Given the growth parameters of the birth rate, we fit the model parameters to the statistics of the oldest person titleholder data using maximum likelihood method. The optimal parameters of the model fit well to the data. It shows in particular that the age of the oldest person alive increases over time, and it will most likely increase further in the future. We compute the expected value and a confidence interval for the age of the world’s oldest person using the fitted model parameters for each year between 1955 and 2019 shown by the green curves on Figure 1. The detailed discussion of the model specification, likelihood calculations as well as the parameter fitting are given in Section 2. Section 3 contains calculations related to the gamma–Gompertz–Makeham generalization of the gamma–Gompertz distribution.
Our results enable us to predict the age distribution of the world’s oldest person at future time points. We compute the probability density of the age of the world’s oldest person in different years not only in the past but also in the future. These densities are shown in Figure 2. When comparing the age distribution of the oldest person in the world in different years to the age of Jeanne Calment at her death, we find that on Jan 1st 2060 we can expect that the age of the world’s oldest person will exceed her age with probability around . This also means that with high probability her age record will already be broken by that time.
In Figure 1 two extreme outliers with unexpectedly long lifetimes can be observed. Jeanne Calment died at the age of 122.45 years in 1997, and Sarah Knauss died at the age 119.27 in 2000. In our model, the probability of observing an age greater than or equal to their actual age at the time of their death is for Calment and for Knauss. See details in Subsection 2.4.
The fact that Calment and Knauss are outliers among the oldest old in the world became even more evident when we performed a backtesting of our model. We estimated the parameters based on the data on the world’s oldest person between 1955 and 1988 where the ending date is the time when Calment became the world’s oldest person. The model-estimated mean and confidence interval of the world’s oldest person using the full data and the partial dataset (before Calment) are shown in Figure 1 by green and red, respectively. The estimate using the data until 1988 is less reliable after 2000 which is shown by the fact that the observed data is out of the confidence interval in the majority of the time after 2000. When we compare the two confidence intervals, we can conclude that, based on the data before 1988, Calment and Knauss already had extremely high ages at their deaths. Adding the remaining data set, the estimated mean age of the world’s oldest person becomes lower. Hence we can conclude that now we can consider the ages of the two outliers between 1988 and 2000 at their death to be more extreme than based on the information available in 1988.
The other important observable is the reign length of a record holder. The numerical value of the expected reign length with our estimated model parameters is in 1955 and it is in 2019. The empirical value of the reign length is which is not much less than the model-based estimate.
Our approach to studying the age of the oldest old is completely novel because it takes into account jointly the age and the time of birth of individuals. Although the age and the reign length of the world’s oldest person depend in a complex non-linear way on the total lifespan and the time of birth of supercentenarians, we compute explicitly the probability distribution of the age of the oldest person. Hence, the performance of our predictions cannot be directly compared to previous results in the literature because in the usual approach, the oldest person in each cohort is considered separately, and it is not relevant whether this person was ever the oldest in the population, see e.g. the extreme value method in [20]. Our model contributes to the mathematical understanding of the evolution of the oldest individual, which is the extra benefit compared to a prediction using the trend in the data, e.g. in linear regression. In this way, we not only observe but also prove mathematically that the dynamics of the birth process and that of the lifespan distribution which we consider in this paper necessarily imply the increase of the expected age of the world’s oldest person.
1 Mathematical model for the age of the oldest person
In this section, we provide the mathematical definition of a general model for the age of the world’s oldest person, where the births of individuals follow a Poisson process and their lifespans are independent. Under the assumptions of the general model, the age of the record holder in the population evolves in time as a Markov process with explicit transition probabilities. In Subsections 1.1–1.2, we describe the exact distribution of the age of the record holder in this generality for any given birth rate parameter and lifespan distribution using the two-dimensional representation of the age process of the oldest person. In the time-homogeneous case with constant birth rate and identical lifespan distributions the reign length distribution of a record holder is computed in Subsections 1.3–1.5. We explain the role of the entry age parameter in Subsection 1.6.
1.1 Model description and two-dimensional representation
The model is formally defined as follows. Let be the birth rate parameter which depends on time and let and be a family of cumulative distribution functions and density functions corresponding to non-negative random variables which are also time-dependent. We assume that individuals are born according to a Poisson point process at rate and that the lifespan of an individual born at time is given by so that lifespans are independent for different individuals.
Let denote the age of the oldest person in the population at time . The process is Markovian. The Markov property holds because at any time the history of the process provides information about the lifetime of individuals born before the current record holder while any transition of depends only on the lifetime of the current record holder and of those born after them.
The evolution of the Markov process is the following. It has a deterministic linear growth with slope due to the ageing of the current record holder. This happens until the death of the record holder. Additionally, given that for some with , the process has a downward jump at time at rate which is the hazard rate of the distribution at . This corresponds to the possibility that the record holder dies at time which happens at rate . The conditional distribution of the jump is given by
(1) |
for all . The jump distribution in (1) has an absolutely continuous part supported on with density
(2) |
and a point mass at with probability
(3) |
As we shall see in the relevant parameter regime the probability of the point mass at is negligible. The transition formula in (1) can be proven using the description below.
We introduce a two-dimensional representation of the process as follows. Let be a marked Poisson process in where forms a Poisson point process on with intensity and is sampled independently for each according to the distribution . The point represents an individual born at time with lifespan for all , that is, the individual is alive in the time interval and their age at time is if . Hence the marked Poisson process contains all relevant information about the age statistics of the population at any time. In particular the age of the oldest person can be expressed in terms of as
(4) |
where the indicator is exactly if the th person is alive at time .
The transition distribution formula (1) can be seen using the two-dimensional representation as follows. Given that the current record holder dies at time at age the event means that nobody with age between and can be alive at time . This event can be equivalently characterized in terms of the Poisson process of birth at rate thinned by the probability that the person is still alive at time . Indeed the event can be expressed as a Poisson process of intensity at time given by for not having any point in the time interval . This probability appears exactly on the right-hand side of (1). In other words, for any people are born at time at rate . On the other hand, the probability for a person born at time to be alive at time (that is, at age ) is . See Figure 3 for illustration.
1.2 Exact distribution of the oldest person’s age process
We assume that all birth events are already sampled on together with the corresponding lifespans. Then the distribution of can be computed explicitly for all using the two-dimensional representation. For all the density
(5) |
for all and the point mass at
(6) |
characterize the distribution of which can be seen as follows. We mention that the point mass at is negligible in the application.
Similarly to the proof of the transition formula in (1) the event for any is the same as the event that nobody with age at least is alive at time . We express this event in terms of the Poisson process of birth at rate thinned by the probability that the person is still alive at time . The event means that a Poisson process of intensity at time given by for does not have any point in yielding
(7) |
In other words for any individuals are born at time at rate . A person born at time is alive at time at age with probability . (5)–(6) follow by differentiation in (7) and by taking the limit. See Figure 4 for illustration.
1.3 Homogeneous model
The exact computation of the reign length distribution (see Subsection 1.5) can only be performed in a special case of our general model described in Subsection 1.1. We introduce this special case as the homogeneous model where individuals are born at the times of a Poisson process of constant rate for all . The lifespan of individuals are independent and identically distributed with a fixed density and cumulative distribution function for all which does not depend on time.
In the homogeneous model the jump distribution given in (2)–(3) simplifies to
(8) | ||||
The distribution of does not depend on time in this case hence it is a stationary distribution as well. The formulas for the density of and point mass at reduce in the homogeneous case to
(9) | ||||
where the integral is equal to the expected lifespan.
The equilibrium condition for the homogeneous density can be written as
(10) |
After differentiation and using the fact that
(11) |
one can derive from (10) the second order differential equation
(12) |
The point mass at satisfies
(13) |
1.4 The peaks process
In the homogeneous model the sequence of peaks in forms a discrete time Markov chain. By peak we mean a local maximum of with value being equal to the lifespan of the last record holder. Each time the oldest person dies the process has a peak with a downward jump following it. Let denote the age of record holders at which they die which are the values of the peaks of the process . The sequence forms a discrete time Markov chain. The Markov property follows by the fact that ages at death of previous record holders only give information on people born before the current record holder but transitions depend on the lifespan of the current record holder and that of people born after them.
The stationary density of is given by
(14) |
The formula can be seen as follows. To have a record holder who dies at age there has to be a person who has lifespan which gives the factor in the numerator on the right-hand side of (14). The exponential factor is by the two-dimensional representation equal to the probability that no people born before the record holder who just died can be alive at the time the record holder dies. The denominator on the right-hand side of (14) makes a probability density function.
The density of can also be characterized by the following description. It satisfies the integral equation
(15) |
which comes from the possible transitions of the peak process as follows. If the previous record holder had a total lifetime then at the death the process jumps down to some value at rate or to with probability . The density of the age at which a person dies who becomes a record holder at age is . From the integral equation in (15) one can derive the second order differential equation for the function given by
(16) |
which is satisfied by in accordance with (14).
1.5 Reign length distribution
In the homogeneous model, let denote the reign length of the th record holder, that is, the time length for which this person is the oldest person of the population. The density of the random reign length is given by
(17) |
The density formula in (17) can be derived based on the stationary density of the peaks process given by (14) as follows.
It holds for the density of the reign length that
(18) |
based on the decomposition with respect to the previous value of the peaks process . The integral is the density of the convolution of the density with an independent exponential distribution of parameter . On the right-hand side of (18), one can use the definitions of the density given by (14) and the homogeneous jump distribution and given by (8). Then in the numerator after the exchange of the order of integrations in the first term and by using the formula for the stationary distribution given by (9) one gets the numerator of (17). In the denominator one can use the equality
(19) |
which follows by integration by parts.
Note also that the density is not equal to the remaining reign length of under the stationary distribution because it would involve the integral in place of the first term in the numerator on the right-hand side of (17).
1.6 Entry age parameter
Next we introduce another parameter which we call the entry age and we denote it by . As opposed to our original model we consider individuals as being born at age at a modified birth rate and with a modified lifespan distribution. As a result we obtain a model to the age of the world’s oldest person all values of the entry age parameter and we can fit the model parameters with different values of the entry age.
For any value of the entry age we denote by the rate at which people reach the age at time , that is,
(20) |
because the new birth process of rate is obtained by an inhomogeneous thinning of the original Poisson process of the birth events. The lifespan distribution of those born at time with age becomes the remaining lifetime distribution at age . The modified cumulative distribution function and density are given by
(21) |
2 Model specification and parameter fitting
In Subsection 2.1 we specify the general model introduced and discussed in Section 1, that is, we assume that the birth rate parameter increases exponentially in time and that the lifespan distribution is given by the gamma–Gompertz–Makeham distribution with time-dependent parameters. We provide the details of the computation of the likelihood as a function of the model parameters in Subsection 2.2. We show the way to maximize the likelihood and how the optimal parameters can be found using the Nelder–Mead method in Subsection 2.3. With these values of the parameters, the age of the world’s oldest person and the reign length of the record holder can be computed as described in Subsection 2.4.
2.1 Model specification: birth rate parameter and lifespan distribution
For the rest of the paper we specifiy our general model described in Section 1 to the following choice of the birth rate parameter and of the lifespan distribution. We choose the value of the entry age to be and we fit the model parameters for all three values of separately.
First we specify the intensity function of the Poisson process of births with an exponential growth in time. For any of the three values of the entry age we assume that the birth rate at age is given by
(22) |
where the numerical values of the parameters and are obtained by linear regression of the logarithmic data of newborns, people at age and published by the United Nations since 1950. We extrapolate the linear regression backwards in time and we use the numerical values shown in Table 1.
We assume that the underlying force of mortality is chosen so that the lifespan distribution of individuals follows the gamma–Gompertz () distribution with cumulative distribution function and density
(23) | ||||
for where are positive parameters. We mention that the gamma–Gompertz–Makeham () distribution differs from the distribution by the presence of a non-negative extrinsic mortality parameter which appears as an additive term in the force of mortality. See (29) for the definition of the distribution. In our model, we exclude the extrinsic mortality for the following two reasons. Since the extrinsic mortality becomes irrelevant at high ages and we aim to model the front-end of the death distribution at the oldest-old ages, we do not expect to obtain a reliable estimate on the extrinsic mortality using the data about the world’s oldest person. On the other hand, as explained later in Subsection 2.3, the likelihood maximization provides unrealistic lifespan distributions even for the model if one tries to optimize in all the parameters at the same time. In order for the algorithm to result in a distribution close to the actual human lifespan distribution, the number of model parameters had to be decreased.
For our model we suppose that in the lifespan distribution, parameters , the rate of aging and , the magnitude of heterogeneity are constants over time and that they only depend on the value of the entry age parameter . The parameter , the initial level of mortality at the entry age for individuals born at time , depends on time given by the exponentially decreasing function
(24) |
where the exponent and the constant only depends on the entry age . The reason for subtracting in (24) is only technical, the numerical values of the parameters do not become tiny with this definition.
In the model with entry age , we assume that the birth rate is given by (22) and we fit the gamma–Gompetz distribution with parameters and given by (24) for the modified distribution function and density in (21). This means that we search for the best fitting values of the parameters which results in an approximation of the remaining lifetime distribution at the age .
2.2 Likelihood calculations
The aim of the maximum likelihood method is to give an estimate to the parameters , , and for by finding those values for which the likelihood of the full sample is the largest. The sample is obtained from the the historical data on the world’s oldest person available in [25]. We transform this information into a list of triples of the form for where is the th time in the sample when the oldest person dies at age and the new record holder has age at time . Then the data has to satisfy the consistency relation since the two sides express the date of birth of the same person.
In the model with entry age , the likelihood of the th data point given the previous data point is equal to
(25) |
for all except for in which case the factor in the denominator is missing. In (25) above we use the transition probabilities of the model with entry age given by
(26) |
as a generalization of (2). The explanation of the left-hand side of (25) is that the person died at time at age had age at time . The previous data point ensures that this person has already reached age hence we condition their lifetime distribution on this fact. The transition probabilities in (26) are obtained similarly to (2) with the difference that a person at age with at time had age at time .
Note that when computing the likelihood of the full data by multiplying the right-hand side of (25) for different values of the consistency relation of the data implies that the factor of the th term cancels with the factor coming from the st term. Hence the log-likelihood of the full sample is given by
(27) | ||||
where we suppress the dependence of the parameters on the entry age. Note that the last two terms do not depend on the parameters hence we can omit these terms in the maximization of the log-likelihood.
2.3 Likelihood maximization
We implemented the calculation of the log-likelihood function given by (27) in Python. We used numerical integration to obtain the integrals on the right-hand side of (27). We mention that the general integral formula in (30) could not be used because the parameter of the gamma–Gompertz distribution in the integrand depends on the integration variable on the right-hand side of (27).
In order to maximize the value of the log-likelihood function we applied the Nelder–Mead method [27] which is already implemented in Python. We mention that initially we used the gamma–Gompertz–Makeham distribution as lifespan distribution, see (29) for the definition, which contains the extra parameter to be fitted but it turned out that the number of model parameters has to be reduced. The behaviour of the optimization algorithm in the five parameters using the gamma–Gompretz–Makeham model was very similar the case of four parameters in the gamma–Gompertz model. Running the optimization in the full set of parameters ( in the gamma–Gompertz–Makeham model or in the gamma–Gompertz model), it turned out that after a few rounds the parameter started to decrease dramatically and reached values below . The resulting lifespan distribution seemed very unrealistic with almost no mortality before the age of 100. This happened for all values of the entry age .
We explain this phenomenon by the fact that historical data about the oldest person in the world only gives information about the behaviour of the lifespan distribution between the ages 107 and 123. The simple optimization in the four parameters simultaneously yields an excellent fit for the tail decay of the lifespan distribution with the historical data but the result may be very far from the actual human lifespan. This would limit the practical relevance of our results.
The mathematical reason for the fact that the four-parameter optimization does not result in a satisfactory approximation to the human lifespan distribution is the following. In these cases, the optimization procedure diverges to those regimes of the parameter space where the corresponding gamma–Gompertz distribution is degenerate. One can prevent reaching these unrealistic combinations of parameters by reducing the amount of freedom in the optimization. Hence we specify some of the parameters a priori and we perform the optimization in the remaining ones so that it provides a good fit to the data on the age of the oldest old as well as a realistic lifespan distribution.
We believe that the most robust of the four parameters of the model is which is the exponent in the time dependence of the mortality rate. By setting the rate of aging the algorithm gives the optimal triple with the best likelihood which is very stable under changing the initial values of these parameters. The running time is also very short.
The Nelder–Mead algorithm, being a numerical maximization method, heavily relies on the tolerance parameter, which determines the minimal improvement required for the algorithm to continue running. If this parameter is set too high, the algorithm might stop before reaching the optimum. Conversely, if set too low, the algorithm might take excessively long to converge. To address this, we drew inspiration from dynamic learning rate algorithms used in neural network training and developed the following meta-algorithm.
First, we run the Nelder–Mead optimization. Based on the improvement from the starting point, we dynamically adjust the tolerance factor, similar to how learning rates are modified during neural network training. We then run the optimization again, recalibrating the tolerance factor based on the observed improvement, and repeat the process. This iterative adjustment allows us to get closer to the optimum, a hypothesis supported by our practical experience with this meta-algorithm. Following this meta-algorithm, only a few calls of the Nelder–Mead method is enough to reach the optimum. The Python codes for the likelihood calculations as well as the Nelder–Mead optimization implemented to this problem are available in [28].
The numerical values of the resulting parameters for the three choices of the entry age are shown in Table 2. The survival probability functions with the parameters given in Table 2 for individuals born in 2000 corresponding to the entry age are shown on Figure 5 as a function of the age. We also computed the optimal values of the parameters for other values of the rate of aging as a sensitivity analysis. The resulting parameter values for the choices , and are shown in Table 3.
We mention as an alternative approach that scaling the parameters could enhance the optimization process, but this requires prior knowledge of the range within which the parameters vary. This range could be determined through our iterative application of the Nelder–Mead algorithm.
2.4 Computation of the oldest person’s age and of the reign length
We observe that the model with the parameters in Table 2 fits well to the titleholder data. We focus on two statistics of the process in order to support this observation about the comparison: the age of the world’s oldest person and the reign length of the record holder. In the case of both statistics exact formulas are only available for the homogeneous model introduced in Subsection 1.3 where the birth rate is constant as well as the lifespan distribution does not depend on time. Hence we apply an approximation where the error is negligible compared to the difference from the statistics computed using the data.
In the general model the distribution of the age of the world’s oldest person at time is given by the density in (5) and by the point mass at in (6). For the numerical computations, we ignore the point mass which is below the round-off error in the numerical results. The difficulty in computing the mean age of the oldest person at time is that parameter of the gamma–Gompertz–Makeham distribution function in the exponent of (5) also depends on the integration variable .
In our approximation we fix the value of the parameter of the distribution in in (5) to a value which is equal to in (24) with some delay . The delay is chosen so that the mean age of the oldest person computed using as parameter for all times in the distribution function in (5) is equal to the same value . For a given , this value of can be obtained as the fixed point of the contraction map
(28) |
which provides a reasonable approximation for the mean age of the world’s oldest person. For the comparison with the data and for the prediction, we use the model with entry age . Hence in (28), the function is given in (22) with and is the distribution function with parameters given by the values in Table 2 and with in (24). This approximation is reasonable because the distribution of the age of the oldest person is highly concentrated.
The fixed point of the map in (28) as the expected age of the world’s oldest person can be found in a few steps of iterations. We show the result on Figure 1. We applied 10 iterations using the fitted model parameters for each year between 1955 and 2019. By computing the standard deviation of the age of the oldest person as well we obtain the mean and a confidence interval for the age.
The predictions for the age distribution of the world’s oldest person in the future shown on Figure 2. We obtained them by using the exact age distribution formula in Subsection 1.2 along with the numerical values of the parameters given in Tables 1 and 2 for entry age .
In our model, the distribution function of the age of the world’s oldest person is given in (7) where the distribution function can be substituted with the estimated parameter values at any time. In this way, the probability of observing an age greater than or equal to Calment’s or Knauss’ actual age at the time of their death can be computed exactly. The numerical values are for Calment and it is for Knauss.
The backtesting mentioned in the Results section is performed as follows. We estimated the best parameter values with entry age based on the reduced data on the world’s oldest person between 1955 and 1988 where the ending date is the time when Calment became the world’s oldest person. The resulting parameters are numerically not very far from the optimal parameters in Table 2 but the difference is more visible on Figure 1. The figure shows the model-based mean age and confidence interval for the age of the world’s oldest person computed using the full data as well as the data until 1988.
For the reign length of record holders, we again used the expected age at a given time obtained as the fixed point of the iteration in (28). The numerical value of the expected reign length obtained from the iteration is in 1955 and it is in 2019. The empirical value of the reign length is computed from the data by dividing the total length of the time interval between 1955 and 2019 by the number of record holders.
3 Methods
In this section, we provide supplementary information related to the main result of this paper. We perform explicit computations with the gamma–Gompertz–Makeham model and we express the integral of the survival function in terms of a hypergeometric function.
The cumulative distribution function and the density of the gamma–Gompertz–Makeham () distribution are given by
(29) | ||||
for where are positive parameters. The positivity of parameters implies the finiteness of all moments and, in particular, the convergence of the integral of the survival function . In the homogeneous model, the integral of the survival function appears in the density of the distribution of in (9) and in the stationary density of the peaks process in (14). We show below that in the gamma–Gompertz–Makeham model the integral of the survival function can be computed explicitly and it is given by
(30) |
where is the hypergeometric function. See 15.1.1 in [29] for the definition and properties.
We prove (30) based on the following integral representation 15.3.1 in [29] of the hypergeometric function
(31) |
which holds whenever . First we prove an identity for complex parameters which satisfy and we compute
(32) | ||||
where we applied a change of variables in the second equality above and we applied the hypergeometric identity (31) in the last equality with , , , together with the observation that with these values of the parameters the prefactor of the integral on the right-hand side of (31) simplifies to . Note that the condition for (31) to hold is satisfied by our assumption which also makes the integrals in (32) convergent.
References
- [1] Oeppen, J. & Vaupel, J. W. Broken limits to life expectancy. \JournalTitleScience 296, 1029–1031 (2002).
- [2] Canudas-Romo, V. Three measures of longevity: Time trends and record values. \JournalTitleDemography 47, 299–312 (2010).
- [3] Vaupel, J. W., Villavicencio, F. & Bergeron-Boucher, M.-P. Demographic perspectives on the rise of longevity. \JournalTitleProceedings of the National Academy of Sciences 118, e2019536118 (2021).
- [4] Gompertz, B. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. \JournalTitlePhilosophical Transactions of the Royal Society of London 115, 513–583 (1825).
- [5] Németh, L. & Missov, T. I. Adequate life-expectancy reconstruction for adult human mortality data. \JournalTitlePlos one 13, e0198485 (2018).
- [6] Missov, T. I., Németh, L. & Dańko, M. J. How much can we trust life tables? sensitivity of mortality measures to right-censoring treatment. \JournalTitlePalgrave Communications 2, 1–10 (2016).
- [7] Maier, H., Jeune, B. & Vaupel, J. W. Exceptional lifespans (Springer Nature, 2021).
- [8] Dang, L. H. K. et al. The question of the human mortality plateau. \JournalTitleDemographic Research 48, 321–338 (2023).
- [9] Vijg, J. & Le Bourg, E. Aging and the Inevitable Limit to Human Life Span. \JournalTitleGerontology 63, 432–434 (2017).
- [10] Alvarez, J.-A., Villavicencio, F., Strozza, C. & Camarda, C. G. Regularities in human mortality after age 105. \JournalTitlePloS one 16, e0253940 (2021).
- [11] Belzile, L. R., Davison, A. C., Gampe, J., Rootzén, H. & Zholud, D. Is there a cap on longevity? A statistical review. \JournalTitleAnnual Review of Statistics and Its Application 9, 21–45 (2022).
- [12] Rootzén, H. & Zholud, D. Human life is unlimited–but short. \JournalTitleExtremes 20, 713–728 (2017).
- [13] Barbi, E., Lagona, F., Marsili, M., Vaupel, J. W. & Wachter, K. W. The plateau of human mortality: Demography of longevity pioneers. \JournalTitleScience 360, 1459–1461 (2018).
- [14] Modig, K., Andersson, T., Vaupel, J., Rau, R. & Ahlbom, A. How long do centenarians survive? Life expectancy and maximum lifespan. \JournalTitleJournal of internal medicine 282, 156–163 (2017).
- [15] Wilmoth, J. R. & Robine, J.-M. The world trend in maximum life span. \JournalTitlePopulation and Development Review 29, 239–257 (2003).
- [16] Gavrilov, L. & Gavrilova, N. Mortality measurement at advanced ages: A study of the social security administration death master file. \JournalTitleNorth American Actuarial Journal 15, 442–447 (2011).
- [17] Newman, S. J. Errors as a primary cause of late-life mortality deceleration and plateaus. \JournalTitlePLoS biology 16, e2006776 (2018).
- [18] Camarda, C. G. The curse of the plateau. measuring confidence in human mortality estimates at extreme ages. \JournalTitleTheoretical Population Biology 144, 24–36 (2022).
- [19] Missov, T. I. & Vaupel, J. W. Mortality implications of mortality plateaus. \JournalTitlesiam REVIEW 57, 61–70 (2015).
- [20] Gbari, S., Poulain, M., Dal, L. & Denuit, M. Extreme value analysis of mortality at the oldest ages: a case study based on individual ages at death. \JournalTitleNorth American Actuarial Journal 21, 397–416 (2017).
- [21] Hanayama, N. & Sibuya, M. Estimating the upper limit of lifetime probability distribution, based on data of japanese centenarians. \JournalTitleJournals of Gerontology Series A: Biomedical Sciences and Medical Sciences 71, 1014–1021 (2016).
- [22] Einmahl, J. J., Einmahl, J. H. & de Haan, L. Limits to human life span through extreme value theory. \JournalTitleJournal of the American Statistical Association 114, 1075–1080 (2019).
- [23] Li, J. & Liu, J. A modified extreme value perspective on best-performance life expectancy. \JournalTitleJournal of Population Research 37, 345–375 (2020).
- [24] Milholland, B. Jeanne Calment, actuarial paradoxography and the limit to human lifespan. \JournalTitleRejuvenation Research 23, 17–18 (2020).
- [25] Group, G. R. World’s oldest person titleholders since 1955. https://grg.org/Adams/C.HTM (2018). [Online; accessed 14-06-2023].
- [26] Brillinger, D. R. A biometrics invited paper with discussion: the natural variability of vital rates and associated statistics. \JournalTitleBiometrics 42, 693–734 (1986).
- [27] Nelder, J. A. & Mead, R. A Simplex Method for Function Minimization. \JournalTitleThe Computer Journal 7, 308–313 (1965).
- [28] Kiss, C. Modelling the age of the oldest person in the world. https://github.com/csabi0312/modelling-the-age-of-the-oldest-person-in-the-world (2024).
- [29] Abramowitz, M. & Stegun, I. A. Pocketbook of Mathematical Functions (Verlag Harri Deutsch, Thun-Frankfurt am Main, 1984).
Acknowledgements
We thank Katalin Kovács for some useful advice at an early stage of the project which led to this collaboration. The work of Cs. Kiss and B. Vető was supported by the NKFI (National Research, Development and Innovation Office) grant FK142124. B. Vető is also grateful for the support of the NKFI grant KKP144059 “Fractal geometry and applications” and for the Bolyai Research Scholarship of the Hungarian Academy of Sciences. L. Németh was supported by MaRDI, funded by the Deutsche Forschungsgemeinschaft (DFG), project number 460135501, NFDI 29/1 “MaRDI – Mathematische Forschungsdateninitiative.
Author contributions
B. V. initiated the research. Cs. K. and B. V. derived the model, prepared the scripts and figures, and carried out the estimations. Cs. K., L. N. and B. V. analyzed the results and wrote the manuscript.
Data availability
The titleholder data are freely available at https://grg.org/Adams/C.HTM
Additional information
The authors declare no competing interests.