Abstract
We use a combination of extreme value statistics, survival analysis and computer-intensive methods to analyse the mortality of Italian and French semi-supercentenarians. After accounting for the effects of the sampling frame, extreme-value modelling leads to the conclusion that constant force of mortality beyond 108 years describes the data well and there is no evidence of differences between countries and cohorts. These findings are consistent with use of a Gompertz model and with previous analysis of the International Database on Longevity and suggest that any physical upper bound for the human lifespan is so large that it is unlikely to be approached. Power calculations make it implausible that there is an upper bound below 130 years. There is no evidence of differences in survival between women and men after age 108 in the Italian data and the International Database on Longevity, but survival is lower for men in the French data.
1. Introduction
Solid empirical understanding of human mortality at extreme age is important as one basis for research aimed at finding a cure for ageing (described, e.g. in [1]), and is also an element in the hotly debated and societally important question whether the current increase in expected lifespan in developed countries, of about three months per year since at least 1840 [2], can continue. The limit to human lifespan, if any, also attracts considerable media attention (e.g. [3–5]).
Einmahl et al. [6] analyse data on mortality in The Netherlands and conclude that ‘there indeed is a finite upper limit to the lifespan’ for both men and women. Their dataset, provided by Statistics Netherlands and consisting of about 285 000 ‘Dutch residents, born in the Netherlands, who died in the years 1986–2015 at a minimum age of 92 years’, had not undergone any validation procedure. As might be expected, the vast majority died before their 100th birthdays: 99.5% lived 107 or fewer years, and 97% died at age 101 or younger. The cohorts for analysis were taken to be the calendar years of death, and truncation of lifetimes was not taken into account. Hanayama & Sibuya [7] also find an upper lifespan limit of 123 years for Japanese persons aged 100 or more, by fitting a generalized Pareto distribution to 1-year and 4-year birth cohorts, taking into account the sampling scheme. In both cases, any plateauing of mortality may be confounded with the increase in hazard between ages 100 and 105, and this would invalidate the extrapolation to yet higher ages.
The validity of conclusions about mortality at extreme age depends crucially on the quality of the data on which they are based [8], as age misrepresentation for the very old is common even in countries with otherwise reliable statistical data [9]. Motivated by this, demographic researchers from 13 countries contribute to the International Database on Longevity (
An earlier release of the
The Italian National Institute of Statistics (
Data analysis must take into account the sampling scheme underlying such data. The
In our analysis, we take the sampling frame into account, pinpoint the age, if any, at which mortality attains a plateau, and disentangle the effects of age and of birth cohort. We also compare mortality in the
We use the generalized Pareto distribution from extreme value statistics in the main analysis, supplemented by fits of the Gompertz distribution, which is standard in demography. We first outline our main results and conclusions; the appendix gives a more detailed description of our methods.
2. Results for ISTAT data
Lifetimes beyond 105 years are highly unusual and the application of extreme value models [14] is warranted. We use the generalized Pareto distribution,
The corresponding hazard function, often called the ‘force of mortality’ in demography, is the density evaluated at excess age x, conditional on survival to then, i.e.
The choice of a threshold u above which equation (2.1) models exceedances appropriately is a basic problem in extreme value statistics and is surveyed by Scarrott & MacDonald [15]. If u is high enough for equation (2.1) to provide an adequate approximation to the distribution of exceedances, then the shape parameter γ is approximately unchanged if a higher threshold u′ is used, and the scale parameters for u and u′ have a known relationship, so a simple and commonly used approach to the choice of threshold is to plot the parameters of the fitted distributions for a range of thresholds [16] and to use the lowest threshold above which parameter estimates stabilize. This choice balances the extrapolation bias arising if the threshold is too low with the increased variance incurred when taking u too high to retain an adequate number of observations.
Figure 2a shows that for age thresholds close to 105 years the estimated shape parameters for excess life lengths are negative, with 95% confidence intervals barely touching zero, but there is no systematic indication of non-zero shape above 107 years. Figure 2b displays the estimated scale parameter of the exponential model fitted to life lengths exceeding the threshold. The scale parameters decrease for ages 105–107 but show no indication of change after age 107, where the scale parameter estimate is 1.45. Parameter stability plots suggest an exponential model and hence a constant hazard after age 107 or so, where a mortality plateau seems to be attained. A more formal analysis supporting such a threshold is given in appendix A.3.
The upper part of table 1 shows results from fitting equation (2.1) and the exponential distribution to the
threshold | 105 | 106 | 107 | 108 | 109 | 110 | 111 |
---|---|---|---|---|---|---|---|
n_{u} | 3836 | 1874 | 947 | 415 | 198 | 88 | 34 |
σ | 1.67 (0.04) | 1.70 (0.06) | 1.47 (0.08) | 1.47 (0.11) | 1.33 (0.15) | 1.22 (0.23) | 1.5 (0.47) |
γ | −0.05 (0.02) | −0.07 (0.03) | −0.02 (0.04) | −0.01 (0.06) | 0.03 (0.09) | 0.12 (0.17) | 0.06 (0.30) |
σ_{e} | 1.61 (0.03) | 1.60 (0.04) | 1.45 (0.05) | 1.45 (0.08) | 1.36 (0.11) | 1.35 (0.17) | 1.58 (0.32) |
p-value | 0.04 | 0.01 | 0.70 | 0.82 | 0.74 | 0.44 | 0.84 |
p_{∞} | 0.02 | 0.01 | 0.35 | 0.41 | 0.63 | 0.78 | 0.58 |
n_{u} | 9808 | 5026 | 2471 | 1208 | 548 | 241 | 106 |
σ | 1.68 (0.02) | 1.58 (0.03) | 1.53 (0.04) | 1.43 (0.06) | 1.37 (0.08) | 1.33 (0.13) | 1.27 (0.18) |
γ | −0.06 (0.01) | −0.04 (0.01) | −0.04 (0.02) | −0.02 (0.03) | 0.01 (0.05) | 0.05 (0.08) | 0.09 (0.11) |
σ_{e} | 1.61 (0.02) | 1.53 (0.03) | 1.48 (0.03) | 1.41 (0.05) | 1.39 (0.07) | 1.38 (0.11) | 1.37 (0.16) |
p-value | 8 × 10^{−7} | 0.01 | 0.06 | 0.60 | 0.78 | 0.46 | 0.32 |
p_{∞} | 4 × 10^{−7} | 4 × 10^{−3} | 0.03 | 0.30 | 0.61 | 0.77 | 0.84 |
The estimated scale parameter obtained by fitting an exponential distribution to the
We investigated effects of year of birth, but found none; see appendix A.4.
3. Results for France 2019 data
Estimation for the
The lower part of table 1 indicates that the exponential and generalized Pareto models fit equally well above 108 years, so we prefer the more parsimonious exponential fit; see appendices A.3 and A.5 and figure 7 for more detail, including a formal check on the suitability of the chosen thresholds. For French persons older than 108, the exponential scale parameter is estimated to be 1.41 (years) with 95% confidence interval (1.32, 1.51), the exponential hazard is estimated to be 0.71 (years^{−1}) with 95% confidence interval (0.66, 0.76) and the estimated probability of surviving at least one more year is 0.49 with 95% confidence interval (0.47, 0.52).
Table 2 shows that estimates of the scale parameter for the exponential distribution for women and men for the
n | σ_{e} (95% CI) | n | σ_{e} (95% CI) | n | σ_{e} (95% CI) | |
---|---|---|---|---|---|---|
women | 375 | 1.45 (1.31, 1.60) | 1110 | 1.46 (1.36, 1.56) | 728 | 1.31 (1.21, 1.41) |
men | 40 | 1.41 (0.98, 1.85) | 94 | 0.90 (0.70, 1.10) | 72 | 1.49 (1.12, 1.86) |
all | 415 | 1.45 (1.31, 1.59) | 1204 | 1.41 (1.32, 1.51) | 800 | 1.33 (1.23, 1.42) |
4. Power
Our analysis above suggests that constant hazard adequately models the lifetimes over 108 years, and extrapolated indefinitely this would imply that there is no limit to the human lifespan. One might wonder whether an increasing hazard would be detectable, however, as the number of persons attaining such ages is relatively small. To assess this, we performed a simulation study described in appendix A.6, mimicking the sampling schemes of the
Any biological limit to their lifespan should be common for all humans, whereas differences in mortality rates certainly arise due to social and medical environments and can be accommodated by letting hazards vary by factors such as country or sex. With the overlap dropped we can treat the datasets as independent and compute the power for a combined likelihood ratio test of γ = 0 (infinite lifetime) against alternatives with γ < 0 (finite lifetime). For concreteness of interpretation, we express the results in terms of the implied upper limit of lifetime, i.e. $\iota =u-\sigma /\gamma $. Figure 3a shows the power curves for the three datasets individually and pooled. The power of the likelihood ratio test for the alternatives $\iota \in \{125,130,135\}$ years, for example, is 0.45/0.32/0.24 for the
Similar calculations give the power for testing the null hypothesis γ = 0 against alternatives γ < 0. Forcing all datasets to have the same shape parameter would allow them to have different endpoints so we reject the overall null hypothesis if we reject the exponential hypothesis, γ = 0, for any of the three datasets. The power of this procedure is also shown in figure 3. The resulting combined power exceeds 0.8 for γ < −0.065 and equals 0.97 for the alternative γ = −0.09, giving strong evidence against a sharp increase in the hazard function after 108 years.
5. Gompertz model
The hazard function of the generalized Pareto distribution cannot model situations in which the hazard increases to infinity but the upper limit to lifetimes is infinite. This possibility is encompassed by the Gompertz distribution [17], which has long been used for modelling lifetimes and often provides a good fit to data at lower ages (e.g. [18]). When the Gompertz model is expressed in the form
The Gompertz distribution has infinite upper limit to its support, so it cannot be used to assess whether there is a finite upper limit to the human lifespan. Its hazard function, σ^{−1} exp(βx/σ), is finite but increasing for all x (β > 0) or constant (β = 0). The limiting distribution for threshold exceedances of Gompertz variables is exponential, and this limit is attained rather rapidly, so a good fit of the Gompertz distribution for lower x would be compatible with good fit of the exponential distribution for threshold exceedances at higher values of x.
Computations summarized in appendix A.7 show that the exponential model, and hence also the Gompertz model with very small β, give equally good fits to the Italian and the French datasets above age 107, and that the Gompertz and generalized Pareto models fit equally well above age 105.
6. Conclusion
Table 2 shows no differences between survival after age 108 in the
There was no indication of differences in survival for women in
There is high power for detection of an upper limit to the human lifespan up to around 130 years, based on fits of the generalized Pareto model to the three databases. This does not mean such ages will be reached sometime soon, however, as the probability of surviving until 130 conditional on reaching 110 years approximately equals that of seeing heads on 20 consecutive tosses of a fair coin. This event has a probability of less than one in a million and is highly unlikely to occur in the near future, though the increasing number of supercentenarians makes it possible that the maximum reported age at death will rise to 130 years during the present century [19].
7. Discussion
The results of the analysis of the newly available
Although many fewer men than women reach high ages, no difference in survival between the sexes is discernible in the
If the
Data accessibility
The
Authors' contributions
L.R.B., A.C.D. and H.R. wrote the paper and appendix. All authors contributed to the statistical analyses. L.R.B. wrote the
Competing interests
We declare we have no competing interests.
Funding
A.C.D. is supported by the Swiss National Science Foundation.
Acknowledgements
We thank contributors to the International Database on Longevity and the Istituto Nazionale di Statistica for providing the data.
Appendix A
Here we describe the methods used to obtain our results, produce the figures and perform our inferences, and add further information on goodness of fit and hazard estimates.
A.1. Reproducibility
Alternative analyses may be performed using the
The
A.2. Truncation and censoring
Figure 1 shows that many persons in the
Consider a sampling interval $\mathcal{C}=(b,e)$ of calendar time during which individuals aged over a threshold of u years were observed. Let x = age − u denote the excess age of an individual who dies aged older than u, having reached age u at calendar time t. Assume that the excess ages x are independent with cumulative distribution function F(x; θ), probability density function f(x; θ), and survival function S(x, θ) = 1 − F(x, θ), where θ is a vector of parameters to be estimated.
The likelihood contribution for someone who died in the interval $\mathcal{C}$ is then
Inappropriate analysis can lead to misleading results: for example, fitting an exponential distribution to the
A.3. Threshold stability
For a formal assessment of threshold stability, we fit a piecewise generalized Pareto model over the K regions above the thresholds u_{1} < · · · < u_{K}, each with its own shape parameter γ_{k} but with scale parameters constrained to ensure that the density function is continuous at the thresholds [22]; this reduces to the usual generalized Pareto model if γ_{1} = · · · = γ_{K}. The p-values for testing γ_{1} = · · · = γ_{K} against the alternative of different values of γ with K = 4 thresholds at the 0, 0.2, … , 0.8 quantiles of the exceedances over 108 or 110 years are shown in table 3. They cast no doubt on the chosen thresholds for either the generalized Pareto model or with a similar test with a piecewise exponential model.
gen. Pareto | exponential | |
---|---|---|
0.47 | 0.61 | |
0.44 | 0.55 | |
0.77 | 0.77 |
A.4. ISTAT cohort effects
The local hazard estimates in figure 4 are obtained by splitting the likelihood contribution of individuals into yearly blocks, using disjoint intervals with (b, e) = (a, a + 1) years to avoid using the same individuals several times. For the highest interval we set e = ∞ and include all individuals who survived into that interval.
The parameter stability plots in figure 5 and the estimated hazard plots in figure 4 show roughly constant mortality for those born in 1886–1905 for the entire age range. Mortality is lower for ages 105 and 106 for persons born in 1906–1910, but equals that for the earlier group at ages 107 and above. This reduction in mortality for the later births implies that plateauing for the entire
A.5. Graphical diagnostics
A quantile-quantile (or QQ-) plot is a standard diagnostic of the fit of a specified distribution to data, but it must be modified to accommodate censoring or truncation. Figure 6a graphs the ordered ages at death, y_{i}, of uncensored
To assess the variability of the plot, we simulate new ages at death from the fitted model, conditioning on the birth dates, truncation time and censoring indicator; this amounts to simulating new lifetimes from a left- and right-truncated exponential distribution, since individuals whose death is observed during the sampling frame cannot exceed the age they would reach at c_{2}, i.e. 31 December 2015 for these data. Both $\hat{F}$ and $\stackrel{~}{G}$ are re-estimated using the simulated samples and evaluated at a grid of fixed values y′ ∈ {y_{1}, …, y_{m}}. The approximate 95% pointwise and simultaneous simulation envelopes shown in the panel are obtained using ${\hat{F}}_{b}^{-1}\{{\stackrel{~}{G}}_{b}({y}^{\prime})\}$ ([27], §4.2.4).
For left-truncated and right-censored data, we can also compare the non-parametric, Nelson–Aalen, conditional cumulative hazard function estimate with its parametric counterparts; uncertainty assessment for the former is discussed in ([23], p. 210). QQ-plotting for left- and right-truncated observations is awkward, but the cumulative hazard can again be estimated fairly readily. Figure 6b and figure 7 show conditional cumulative hazard plots for the
A.6. Power
To assess the statistical power of our procedures, we computed the maximum profile likelihood estimates ${\hat{\sigma}}_{\gamma}$ for the original data with γ fixed at the values −0.25, … , 0 and simulated new excess lifetimes for each γ, conditioning on the sampling frame. In the
We proceeded similarly for the endpoint $\iota $. In order to simulate from generalized Pareto distributions with $\iota $ fixed at a given value, we reparametrized the log-likelihood function in terms of $\iota $ and σ and estimated the scale parameters ${\hat{\sigma}}_{\iota}$ for each of the original three datasets with $\iota $ fixed, then used the relation $\iota =u-\sigma /\gamma $ to obtain the three implied shape parameters ${\hat{\gamma}}_{\iota}$. We then simulated new datasets with the sampling scheme described above, but using the three sets of parameter pairs $({\hat{\gamma}}_{\iota},{\hat{\sigma}}_{\iota})$. For the joint test of the null hypothesis of an exponential distribution, $\iota =\mathrm{\infty}$, against the alternative of a common but finite $\iota $, we reparametrized the likelihood in terms of $\iota $ and allowed the three datasets to have different values of σ.
Figure 3a shows how the empirical proportion of rejections for a test of nominal size 5% based on the directed likelihood root statistic $r=-\{2{\ell}_{\mathsf{g}\mathsf{p}}(\hat{\iota},\hat{\mathit{\sigma}})-2{\ell}_{\mathsf{e}\mathsf{x}\mathsf{p}}({\stackrel{~}{\mathit{\sigma}}}_{e}){\}}^{1/2}$ varies with $\iota =u-\sigma /\gamma $ for γ < 0, for the
A.7. Gompertz model
The reciprocal hazard function of the Gompertz distribution r(x) = σ exp(− βx/σ) encodes the speed of convergence to the limiting extreme value distribution [30]; even if β > 0, r′(x) = β exp(− βx/σ) → 0 exponentially fast as x → ∞. Any improvement in the fit of the Gompertz model for exceedances over some threshold would be shown by evidence that β > 0. This distribution places a point mass of exp(1/β) at x = ∞ when β < 0, so we allow only β ≥ 0.
Table 4 summarizes the fit of this model for various thresholds and the
threshold | 105 | 106 | 107 | 108 | 109 | 110 |
---|---|---|---|---|---|---|
ISTAT | ||||||
n_{u} | 3836 | 1874 | 946 | 415 | 198 | 88 |
σ | 1.67 (0.05) | 1.71 (0.07) | 1.48 (0.08) | 1.47 (0.12) | 1.36 (0.1) | 1.35 (0.2) |
β | 0.05 (0.02) | 0.09 (0.04) | 0.02 (0.05) | 0.02 (0.07) | 0 | 0 |
p-value | 0.02 | 0.01 | 0.37 | 0.45 | 1 | 1 |
France 2019 | ||||||
n_{u} | 9808 | 5026 | 2471 | 1208 | 548 | 241 |
σ | 1.7 (0.03) | 1.59 (0.03) | 1.54 (0.05) | 1.43 (0.06) | 1.39 (0.1) | 1.38 (0.1) |
β | 0.07 (0.01) | 0.05 (0.02) | 0.05 (0.03) | 0.02 (0.04) | 0 | 0 |
p-value | 0 | 0 | 0.02 | 0.31 | 1 | 1 |
In general, p-values obtained using a parametric bootstrap are more reliable than asymptotic approximations such as $\frac{1}{2}{\chi}_{0}^{2}+\frac{1}{2}{\chi}_{1}^{2}$ and are preferable for such comparisons. In the present case, this entails simulating independent datasets like the original data from the boundary exponential distribution (the null hypothesis, β = 0), and estimating the p-value by the empirical proportion of these simulated datasets in which the likelihood ratio statistics w* are no smaller than the original value w, i.e. ${\hat{Pr}}_{0}^{\ast}({w}^{\ast}\ge w)$, where the asterisk indicates a parametric bootstrap simulation and the subscript indicates that the simulation is under the null hypothesis; see ch. 4 of [27]. This approach was used to obtain the p-values in table 4, which show that the exponential model is statistically indistinguishable from the Gompertz model from 108 years onwards, though the Gompertz model fits better at 106 years and below for the
Table 5 compares the fits of the Gompertz, generalized Pareto and exponential models to the
threshold | 105 | 106 | 107 | 108 | 109 | 110 |
---|---|---|---|---|---|---|
Gompertz | 0.01 | 1.65 | 2.21 | 0.22 | 0.68 | 1.65 |
GP | 0.00 | 2.17 | 2.21 | 0.23 | 0.56 | 1.05 |
exponential | 4.16 | 8.21 | 2.36 | 0.28 | 0.68 | 1.65 |
A.8. Differences between men and women
The imbalance in the number of men and women limits our ability to detect any effects of gender on mortality at great age. To illustrate this, we conducted a simulation study based on the
For the
A.9. Hazard estimates
To construct a local hazard estimate based on the limiting generalized Pareto distribution, we note that this distribution has reciprocal hazard function r(x) = (σ + γx)_{+}, where a_{+} = max (a, 0) for real a. A more flexible functional form is r(x) = {σ + γx + g(x)}_{+}, where g(x) is a smooth function of x that tends to zero as x increases. For exploratory purposes, we take $g(x)=\sum _{k=1}^{K}{\beta}_{k}{b}_{k}(x)$, where
Lifetime data are recorded to the nearest day and often there are ties for small x, so we use a discrete version of the model, with x ∈ {1, … , x_{max}} days, where x_{max} corresponds to 16 years after age 105. We let h(x) = 1/r(x) denote the hazard function, set $H(x)=\sum _{z=1}^{x}h(z)/365$ for compatibility with the continuous case, and obtain survivor and probability mass functions $Pr(X>x)=\mathrm{exp}\{-H(x)\}$ and $Pr(X=x)=h(x)\mathrm{exp}\{-H(x)\}$ for x ∈ {1, …, x_{max}}.
Let X denote a survival time (days) beyond 105 years. For each individual, the available data are of the form (s, d, x), where s = 0 if observation of X began at 105 years and, if not, s > 0 is the age in days above 105 at which observation of X began, d = 1 indicates death, X = x, and d = 0 indicates right-censoring, X > x. The likelihood is then a product of terms of the form
Figure 8a shows a local estimate of the hazard function for the
References
- 1.
Vijg J, Campisi J . 2008 Puzzles, promises, and a cure for ageing. Nature 544, 1065-1071. (doi:10.1038/nature07216) Crossref, ISI, Google Scholar - 2.
Oeppen J, Vaupel JW . 2002 Broken limits to life expectancy. Science 296, 1029-1031. (doi:10.1126/science.1069675) Crossref, PubMed, ISI, Google Scholar - 3.
Guarino B . 2018 Good news for human life spans—at age 105, death rates suddenly stop going up. Washington Post. 28 June. Google Scholar - 4.
Saplakoglu Y . 2018 Have humans reached their limit on life span? These researchers say no. Live Science. 28 June. Google Scholar - 5.
Hoad P . 2019 ‘People are caught up in magical thinking’: was the oldest woman in the world a fraud? The Guardian. 30 November. Google Scholar - 6.
Einmahl JJ, Einmahl JH, de Haan L . 2019 Limits to human life span through extreme value theory. J. Am. Stat. Assoc. 114, 1075-1080. (doi:10.1080/01621459.2018.1537912) Crossref, ISI, Google Scholar - 7.
Hanayama N, Sibuya M . 2016 Estimating the upper limit of lifetime probability distribution, based on data of Japanese centenarians. J. Gerontol. A 71, 1014-1021. (doi:10.1093/gerona/glv113) Crossref, Google Scholar - 8.
Gavrilov LA, Gavrilova NS . 2019 Late-life mortality is underestimated because of data errors. PLOS Biol. 17, 1-7. (doi:10.1371/journal.pbio.3000148) Crossref, ISI, Google Scholar - 9.
Poulain M . 2010 On the age validation of supercentenarians. In Supercentenarians (eds H Maier, J Gampe, B Jeune, JM Robine, JW Vaupel), pp. 4–30. Heidelberg, Germany: Springer. Google Scholar - 10.
Gampe J . 2010 Human mortality beyond age 110. In Supercentenarians (eds H Maier, J Gampe, B Jeune, JM Robine, JW Vaupel), pp. 219–230. Heidelberg, Germany: Springer. Google Scholar - 11.
Rootzén H, Zholud D . 2017 Human life is unlimited—but short (with discussion). Extremes 20, 713-728. (doi:10.1007/s10687-017-0305-5) Crossref, ISI, Google Scholar - 12.
Rootzén H, Zholud D . 2018 Rejoinder to discussion of the paper ‘Human life is unlimited—but short’. Extremes 21, 415-424. (doi:10.1007/s10687-018-0325-9) Crossref, ISI, Google Scholar - 13.
Barbi E, Lagona F, Marsili M, Vaupel JW, Wachter KW . 2018 The plateau of human mortality: demography of longevity pioneers. Science 360, 1459-1461. (doi:10.1126/science.aat3119) Crossref, PubMed, ISI, Google Scholar - 14.
de Haan L, Ferreira A . 2006 Extreme value theory: an introduction. Berlin, Germany: Springer. Crossref, Google Scholar - 15.
Scarrott C, MacDonald A . 2012 A review of extreme-value threshold estimation and uncertainty quantification. REVSTAT—Stat. J. 10, 33-60. ISI, Google Scholar - 16.
Davison AC, Smith RL . 1990 Models for exceedances over high thresholds (with discussion). J. R. Stat. Soc. B 52, 393-442. (doi:10.1111/j.2517-6161.1990.tb01796.x) Google Scholar - 17.
Gompertz B . 1825 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. Phil. Trans. R. Soc. Lond. 115, 513-585. Link, Google Scholar - 18.
Thatcher AR . 1999 The long-term pattern of adult mortality and the highest attained age (with discussion). J. R. Stat. Soc. A (Stat. Soc.) 162, 5-43. (doi:10.1111/rssa.1999.162.issue-1) Crossref, PubMed, ISI, Google Scholar - 19.
Pearce M, Raftery AE . 2021 Probabilistic forecasting of human maximum lifespan by 2100 using Bayesian population projections. Demogr. Res. 44, 1271-1294. (doi:10.4054/DemRes.2021.44.52) Crossref, ISI, Google Scholar - 20. R Core Team. 2021 R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Google Scholar
- 21.
Davison AC . 2018 ‘The life of man, solitary, poore, nasty, brutish, and short’: discussion of the paper by Rootzén and Zholud. Extremes 21, 365-372. (doi:10.1007/s10687-018-0329-5) Crossref, ISI, Google Scholar - 22.
Northrop PJ, Coleman CL . 2014 Improved threshold diagnostic plots for extreme value analyses. Extremes 17, 289-303. (doi:10.1007/s10687-014-0183-z) Crossref, ISI, Google Scholar - 23.
Andersen P, Borgan O, Gill R, Keiding N . 1993 Statistical models based on counting processes. New York, NY: Springer. Crossref, Google Scholar - 24.
Waller LA, Turnbull BW . 1992 Probability plotting with censored data. Am. Stat. 46, 5-12. ISI, Google Scholar - 25.
Tsai WY, Jewell NP, Wang MC . 1987 A note on the product-limit estimator under right censoring and left truncation. Biometrika 74, 883-886. (doi:10.1093/biomet/74.4.883) Crossref, ISI, Google Scholar - 26.
Shen PS . 2010 Nonparametric analysis of doubly truncated data. Ann. Inst. Stat. Math. 62, 835-853. (doi:10.1007/s10463-008-0192-2) Crossref, ISI, Google Scholar - 27.
Davison AC, Hinkley DV . 1997 Bootstrap methods and their application. Cambridge, UK: Cambridge University Press. Crossref, Google Scholar - 28.
Turnbull BW . 1976 The empirical distribution function with arbitrarily grouped, censored and truncated data. J. R. Stat. Soc. B 38, 290-295. (doi:10.1111/j.2517-6161.1976.tb01597.x) Google Scholar - 29.
Nair VN . 1984 Confidence bands for survival functions with censored data: a comparative study. Technometrics 26, 265-275. (doi:10.1080/00401706.1984.10487964) Crossref, ISI, Google Scholar - 30.
Smith RL . 1987 Approximations in extreme value theory. Technical Report 205, Center for Stochastic Processes, University of North Carolina at Chapel Hill. Google Scholar - 31.
Self SG, Liang KY . 1987 Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. J. Am. Stat. Assoc. 82, 605-610. (doi:10.1080/01621459.1987.10478472) Crossref, ISI, Google Scholar - 32.
Efron B, Tibshirani RJ . 1993 An introduction to the bootstrap. New York, NY: Chapman & Hall. Crossref, Google Scholar