Reanalysis of the evidence for a limit to human lifespan

Recently a paper was published in the journal Nature titled “Evidence for a limit to human lifespan” which received wide publicity, by nature itself, the New York Times, the Atlantic, the Guardian, and the BBC to name but a few. However some of the methods in this paper have been criticised:

This @nature paper should be renamed ‘Evidence for a limit to Peer-Review’

(ht @StuartBuck1)https://t.co/OsHzPEZScV

— Amitabh Chandra (@amitabhchandra2) October 7, 2016

So, Nature paper on ‘human lifespan limit’ makes this inference from those TWO points at the very end!? Total toss. https://t.co/p1L8FnH9Cw pic.twitter.com/v7P9YMrJxM

— Stuart Ritchie (@StuartJRitchie) October 7, 2016

In this analysis I look at figure 2 specifically which argues that the maximum age of death has plateaued.

I downloaded the data from the International Database on Longevity at the Max Planck Institute for Demographic Research. The terms of the data access do not permit third party sharing so the raw data is not uploaded to GitHub but you can download it yourself if you want to rerun the following analyses.

First I load the data into R, tidy up some of the columns, and subset to the same individuals used the in paper. (Not sure why they didn’t just use all 668 rather than just 534). Here is the breakdown by country:

library(data.table)
library(ggplot2)
library(extRemes) # Generalised Extreme Value distribution functions
library(scales) # to define double log transform
# load in data
data.files <- list.files(path = "raw-data")
stopifnot(length(data.files) > 0)
age.data <- data.table()
i = 0
for (file in data.files){
  i <- i + 1
  temp <- fread(paste0("raw-data/", file))
  
  age.data <- rbind(age.data, temp)
}
# correct column names
setnames(age.data, make.names(colnames(age.data)))
# extract year of death only
age.data$death_year <- as.numeric(gsub("[0-9]+/[0-9]+/", "", age.data$Date.of.death))
age.data$age_inyears <- age.data$Age.days. / 365.25
full.age.data <- age.data
# subset to only people with year of death in certian countries as used by the paper
age.data <- age.data[Country.of.death %in% c("FRA","GBR","JPN","USA")]
# check the same number of people as in the paper
age.data[, .N, by = Country.of.death]

Now let’s recreate figure 2A.

Raw Data

# define function for getting 2nd highest value, 3rd highest etc.
maxN <- function(x, N = 2){
  len <- length(x)
  if(N > len){
    warning('N greater than length(x)')
    return(as.numeric(NA))
  }
  sort(x, partial = len - N + 1)[len - N + 1]
}
# calculate mean, Xth maximum etc summary values for each year of death
summarised <- age.data[,.(mean_age = mean(age_inyears),
                    n = .N,
                    sem = sd(age_inyears) / .N,
                    maximum = max(age_inyears),
                    second_maximum = maxN(age_inyears, 2),
                    third_maximum = maxN(age_inyears, 3),
                    fourth_maximum = maxN(age_inyears, 4),
                    fifth_maximum = maxN(age_inyears, 5),
                    sixth_maximum = maxN(age_inyears, 6)),
                 by = death_year]
# Repoduce Figure 2A
MRAD <- ggplot(summarised, aes(x=death_year, y=maximum)) +
  geom_point() +
  xlab("Year of Death") +
  ylab("Maximum Reported Age at Death")
MRAD

With regression lines

MRAD +
  geom_smooth(method = "lm", data = summarised[death_year > 1994]) +
  geom_smooth(method = "lm", data = summarised[death_year <= 1994])

The authors of the paper fitted two separate regression lines to this data arguing that after 1995 there was a change in the trend (a seemingly arbitrary choice of breakpoint - the choice of a broken vs linear trend has been analysed elsewhere).

You can see from the confidence intervals on the regression lines that the gradient for the second segment is actually consistent with being the same as the first segment. In the paper the authors calculate a p-value of 0.27 for the gradient of the second segment (null hypothesis = 0) and conclude “no further increases were observed”. They apply the same reasoning in a reply to a post-publication review on publons “The latter is not significant, so we conclude that the MRAD is essentially flat”. However, you can not accept the null hypothesis based on p > 0.05, you can only reject a null hypothesis. In this case a p-value of greater than 0.05 suggests that there is not enough data to conclude that the gradient is different from 0 (perhaps the null hypothesis should really by that the gradient is the same as the first segment, although the p-value is still non significant). The 95% confidence interval for the second segment gradient is −0.83 to +0.20 which includes the point estimate of the first segment gradient of 0.15 (using non rounded age values here).

# calculate 95% confidence intervals on the second segment gradient
confint(lm(maximum ~ death_year, summarised[death_year > 1994]), "death_year")

                2.5 %    97.5 %
death_year -0.8276784 0.1956178

cat("First segment point estimate: \n"); coef(lm(maximum ~ death_year, summarised[death_year <= 1994]))['death_year']

First segment point estimate: 
death_year 
 0.1533292 

cat("P-value when H0 = 0.1533: \n"); summary(lm(maximum ~ death_year, summarised[death_year > 1994], offset = 0.1533292 * death_year))$coefficients["death_year", "Pr(>|t|)", drop = FALSE]

P-value when H0 = 0.1533: 
             Pr(>|t|)
death_year 0.06819191

However this analysis is quite sensitive to the choice of breakpoint. In the above mentioned review response the authors re-analysed the data and found that a breakpoint of 1999 was a better fit. Although the package used (“segmented”) fits a continuous piecewise regression, I will continue using the method above to illustrate a different choice of date anyway. Replotting the regression lines using this breakpoint shows the confidence intervals more clearly supporting the downward trend and the p-value (with H0 = 1st segment gradient) is now significant at the 0.05 threshold. The upper 95% confidence interval is now −0.2 suggesting a downward trend (rather than a plateau).

MRAD +
  geom_smooth(method = "lm", data = summarised[death_year <= 1998]) +
  geom_smooth(method = "lm", data = summarised[death_year > 1998])

  
first.segment <- lm(maximum ~ death_year, summarised[death_year <= 1998])
second.segment <- lm(maximum ~ death_year, summarised[death_year > 1998])
# calculate 95% confidence intervals on the second segment gradient
cat("First segment gradient point estimate & confidence intervals: \n"); coef(first.segment)['death_year']; confint(first.segment, "death_year")

First segment gradient point estimate & confidence intervals: 
death_year 
 0.1926284 
                2.5 %   97.5 %
death_year 0.09503971 0.290217

cat("Second segment gradient point estimate & confidence intervals: \n"); coef(second.segment)['death_year']; confint(second.segment, "death_year")

Second segment gradient point estimate & confidence intervals: 
death_year 
-0.6874939 
               2.5 %     97.5 %
death_year -1.164361 -0.2106265

cat("P-value for second segment when H0 = 0.193: \n"); summary(lm(maximum ~ death_year, summarised[death_year > 1998], offset= 0.1926284 * death_year))$coefficients["death_year", "Pr(>|t|)", drop = FALSE]

P-value for second segment when H0 = 0.193: 
              Pr(>|t|)
death_year 0.004034243

Higher order maximums (2nd, 3rd etc)

The authors note that due to the fact that each of these data points is just a single individual the apparent plateau they observe could be due to random fluctuation. To strengthen their argument they looked at the 2nd highest reported age at death, 3rd highest etc and claimed that these series showed the same pattern. However the data points were only plotted for the 1st MRAD and only cubic smoothing splines for the remaining. Fitting a cubic spline could be misleading / overfitting and each series should probably be processed in the same manner as figure 2A if one is to conclude that they show the same pattern. Below I plot each series individually so the actual data is visible. The cubic splines show downward trends towards the end although with increasing uncertainty and linearity. Similarly with the linear regressions the gradient of the second segments are lower than the first segments although with increasing consistency between the two (note variable y-axis).

Raw Data

# Figure 2B
summary_melt <- melt(summarised, id.vars = c("death_year","mean_age","n","sem"))
summary_melt <- summary_melt[!is.na(value)]
# Recreate figure 2A for each maximum type
higher_orders <- ggplot(summary_melt, aes(death_year, value)) +
  facet_wrap(~ variable, ncol = 2, scales = "free") +
  geom_point() +
  xlab("Year of Death") +
  ylab("Age at Death")
higher_orders

With Cubic Spline

higher_orders +
  geom_smooth(method = "lm", formula = y ~ splines::bs(x, degree = 3), size = 0.5)

With Regression Lines

higher_orders +
  geom_smooth(method = "lm", data = summary_melt[death_year > 1994], size = 0.5) +
  geom_smooth(method = "lm", data = summary_melt[death_year <= 1994], size = 0.5)

Mean age of death

In another alternate approach the authors looked at all individuals in the dataset to calculate mean age of death and concluded that the annual average age of supercentenarians had not increased since 1968 (the start of the dataset). I recreate their plot below but with the addition of error bars representing the standard error of the mean for each point in order to visualise the uncertainty in the values.

# Figure 2C
ggplot(summarised, aes(x = death_year, y = mean_age)) +
  geom_errorbar(aes(ymin = mean_age-sem, ymax = mean_age+sem), width = 0.1) +
  geom_point() +
  xlab("Year of Death") +
  ylab("Mean Age at Death") +
  geom_smooth(method = "lm", formula = y ~ splines::bs(x, degree = 3), se=FALSE)

You can see that for the earlier points there are no error bars, this is because there is only a single data point for those years. It is therefore quite misleading to give each mean equal weighting by fitting a cubic spline to point estimates of the means alone.

A perhaps fairer approach is to recreate the graphs but using the whole dataset (note the dataset does not include anyone who died younger than 110). In this form the uncertainty in the first and last few years is much clearer, and the dip pattern fitted above is much less convincing. I would argue that a linear regression fits the data just as well and this gives an increase of ~ 0.04 years per year.

Cubic Spline

all_data <- ggplot(age.data, aes(x = death_year, y = age_inyears)) +
  geom_point(alpha = 0.3) +
  xlab("Year of Death") +
  ylab("Age at Death")
#plot(age.data$death_year, age.data$age_inyears)
#lines(smooth.spline(age.data$death_year, age.data$age_inyears), col="red")
all_data +
  geom_smooth(method = "lm", formula = y ~ splines::bs(x, degree = 3))

Linear Regression

all_data + 
    geom_smooth(method = "lm", colour="red") + 
    annotate("text", x = 1975, y = 121,
             label = paste("gradient =",format(coef(lm(age_inyears ~ death_year, age.data))[2], digits = 2)))

Sample Sizes

In the study the authors analysed maximum reported age of death (MRAD) over different years but the data for each year was from a different combination of countries and hence the sample size varies. One therefore might expect that the MRAD could change solely due to variation in the sample size (we are more likely to see high maximums when there is more data). Here I investigate the effect of using different sample sizes on the MRAD.

To get an equation for the distribution of age at death we can fit a generalised extreme value distribution to data from the UK Office of National Statistics (which fits much better than a normal distribution).

# Take data from UK Office of National Statistics life tables to fit distribution
# http://www.ons.gov.uk/generator?uri=/peoplepopulationandcommunity/birthsdeathsandmarriages/deaths/articles/mortalityinenglandandwales/2012-12-17/74c4dc50&format=csv
# And http://www.ons.gov.uk/generator?uri=/peoplepopulationandcommunity/birthsdeathsandmarriages/deaths/articles/mortalityinenglandandwales/2012-12-17/4b665795&format=csv
# load data
life_expectancy <- fread("chart7628369219653552718.csv") # females
life_expectancy_male <- fread("chart182252582428823712.csv") # males
# convert to numeric
life_expectancy$`2010` <- as.numeric(life_expectancy$`2010`) + as.numeric(life_expectancy_male$`2010`)
life_expectancy$Age <- as.numeric(life_expectancy$Age)
# convert histogram values to sample data
deaths <- numeric()
for (Age2 in life_expectancy$Age){
  deaths <- c(deaths, rep(as.numeric(Age2), life_expectancy[Age == Age2]$`2010`))
}
# fit to a generalised extreme value distribution
fitted_distribution <- fevd(-deaths)
fitted_parameters <- fitted_distribution$results$par
plot(fitted_distribution, type = "density", xlab = "Negative Age at death", main = "")

plot(fitted_distribution, type = "qq2")

params <- fitted_distribution$results$par
print(fitted_distribution$results$par)

   location       scale       shape 
-86.6589101   9.8271019   0.0372343 

We will also need to estimate the sample size (number of deaths) for each year in each country. For this I multiplied the world bank crude death rate by population size. We can then see how the total sample size varies over time in the original papers analysis.

# Calculate sample size as number of deaths for each year-country, by multiplying crude death rate by populatoin size
# from http://api.worldbank.org/v2/en/indicator/SP.DYN.CDRT.IN?downloadformat=csv
death.rate <- fread("API_SP.DYN.CDRT.IN_DS2_en_csv_v2.csv")
death.rate <- death.rate[`Country Code` %in% c("GBR","FRA","JPN","USA"),]
death.rate$`Indicator Name` <- NULL
death.rate$`Indicator Code` <- NULL
death.rate$`Country Name` <- NULL
death.rate <- melt(death.rate, id.vars = c("Country Code"), variable = "year", value = "death_rate")
setnames(death.rate,"Country Code","country")
death.rate$year <- as.numeric(as.character(death.rate$year))
# from http://api.worldbank.org/v2/en/indicator/SP.POP.TOTL?downloadformat=csv
population <- fread("API_SP.POP.TOTL_DS2_en_csv_v2.csv")
population <- population[`Country Code` %in% c("GBR","FRA","JPN","USA"),]
population$`Indicator Name` <- NULL
population$`Indicator Code` <- NULL
population$`Country Name` <- NULL
population <- melt(population, id.vars = c("Country Code"), variable = "year", value = "population")
setnames(population, "Country Code", "country")
population$year <- as.numeric(as.character(population$year))
setkey(population, "year", "country")
setkey(death.rate, "year", "country")
deaths <- population[death.rate]
deaths <- deaths[(country=="GBR" & year >= 1968 & year <= 2006) | 
             (country=="FRA" & year >= 1987 & year <= 2003) |
             (country=="USA" & year >= 1980 & year <= 2003) |
             (country=="JPN" & year >= 1996 & year <= 2007)]
deaths$population <- as.numeric(deaths$population)
deaths$death_rate <- as.numeric(deaths$death_rate)
deaths$deaths <- round((deaths$death_rate / 1000) * deaths$population)
deaths.summary <- deaths[, .(deaths = sum(deaths)), by = year]
qplot(deaths.summary$year, deaths.summary$deaths,
      geom = "point",
      xlab = "Year",
      ylab = "Total deaths")

The trend is similar to the regression lines they fit and so any bias from sample size would result in an overestimate in their favour for the gradient of both of the regression lines. However the effect of sample size on MRAD is probably not linear - maybe the population sizes used are large enough that the MRAD is effectively independent. With the sample size and an equation for the distribution of age at death we can now calculate the probability distribution of MRAD (more formally the nth order statistic) for different sample sizes. First let us look at the distributions of MRAD for the estimated minimum and maximum sample size used in the study.

# define k=1 (minimum) order statistic of minimum extream value distribution
order_pdf_k1_evd <- function(age, sample.size){
  sample.size *
    (1-pevd(age, loc = params['location'], scale = params['scale'], shape = params['shape']))^(sample.size - 1) * 
    devd(age, loc = params['location'], scale = params['scale'], shape = params['shape'])
}
max.sample.size = max(deaths.summary$deaths)
min.sample.size = min(deaths.summary$deaths)
# compare the probability density distributions for the minimum and maximum sample sizes used
age = seq(-109, -115, -0.02)
max.min.comparison <- data.table(
  density = c(order_pdf_k1_evd(age, sample.size = max.sample.size),
              order_pdf_k1_evd(age, sample.size = min.sample.size)),
  maximum.age = -c(age, age),
  sample.size = as.factor(rep(c(max.sample.size, min.sample.size), each = length(age))))
# find modes of both distributions
mode.min = -optimize(order_pdf_k1_evd,
                     interval = c(-100, -120),
                     sample.size = min.sample.size,
                     maximum = TRUE)$maximum
mode.max = -optimize(order_pdf_k1_evd,
                     interval = c(-100, -120),
                     sample.size = max.sample.size,
                     maximum = TRUE)$maximum
ggplot(max.min.comparison, aes(x = maximum.age, y = density, group = sample.size, colour = sample.size) ) +
  geom_line() +
  geom_vline(xintercept = mode.max, colour="turquoise", linetype = "dashed") +
  geom_vline(xintercept = mode.min, colour="red", linetype = "dashed") +
  annotate("text", label = paste(signif(mode.max - mode.min, digits = 3), "years"), x = 114, y = 0.6) +
  xlab("Maximum Reported Age at Death") + ylab("Probability Density")

This shows we might expect a difference of over a year in the MRAD due to the change in sample size alone (dashed lines indicate mode). We can also look at how the modal MRAD changes over many different sample sizes.

# plot most likely maximum observed value as a function of sample size
# Linear on log scale sampling
sample.sizes = matrix(unique(round(exp(runif(300, 0, 16)))))
modal_MRAD <- apply(sample.sizes, 1,function(x) optimize(order_pdf_k1_evd, interval = c(-90, -120), sample.size = x, maximum = TRUE)$maximum)
MRAD_samplesize <- data.table(sample.size = sample.sizes[, 1], modal.maximum = - modal_MRAD)[order(-sample.size)]
ggplot(MRAD_samplesize, aes(x = sample.size, y = modal.maximum)) +
  geom_line() +
  ylim(108, 113) + 
  geom_vline(xintercept = min.sample.size, linetype = "dashed") +
  geom_vline(xintercept = max.sample.size, linetype = "dashed") +
  xlab("Sample Size") + ylab("Modal Maximum Age")

The modal MRAD increases sharply at first and then starts to plateau once the sample size increases to millions of deaths. The dashed lines indicate the estimated minimum and maximum sample sizes used in the study. A double log distribution fits this curve well for reasonable sample sizes (>20).

# plot on double log x scale (straight line)
# define transformation
double_log_trans <- function(base = exp(1)){
trans <- function(x) log((log(x, base)), base)
inv <- function(x) (base^(base^x))
trans_new("double_log", trans, inv, log_breaks(base = base), domain = c(1e-100, Inf))
}
# plot
ggplot(MRAD_samplesize[sample.size > 20], aes(x = sample.size, y = modal.maximum)) +
  geom_point() +
  scale_x_continuous(trans = "double_log", breaks = 10^(2:7)) +
  xlab("Sample Size (double log scale)") +
  ylab("Modal maximum age") +
  geom_smooth(method = "lm", size = 0.5)

We can also plot the difference from the mean MRAD for each year in the study based on changing sample size alone.

sample.sizes.actual = matrix(deaths.summary$deaths)
MRAD_expected <- apply(sample.sizes.actual, 1, function(x) optimize(order_pdf_k1_evd, interval = c(-90, -120), sample.size = x, maximum = TRUE)$maximum)
MRAD_expected_by_samplesize <- data.table(sample.size = sample.sizes.actual[, 1], modal.maximum = - MRAD_expected)
MRAD_expected_by_samplesize$death_year <- deaths.summary$year
# to adjust points?
# MRAD_expected_by_samplesize$demeaned <- as.numeric(scale(MRAD_expected_by_samplesize$modal.maximum, scale = FALSE))
# setkey(MRAD_expected_by_samplesize, death_year)
# setkey(summarised, death_year)
# summarised <- MRAD_expected_by_samplesize[summarised]
# summarised$adjusted_maximum <- summarised$maximum - summarised$demeaned
qplot(MRAD_expected_by_samplesize$death_year,
      scale(MRAD_expected_by_samplesize$modal.maximum, scale = FALSE),
      xlab = "Year",
      ylab = "Mean-Centered MRAD")

Hence the sample sizes used would probably have a noticeable although small effect on the MRAD and a correction would slightly weaken the authors conclusions by reducing the gradient of both regression lines. Even though the effect is moderate it would have been nice to see an analysis of this type reported in the study.

Whether or not there is a genuine plateau rather than a temporary fluctuation would be clearer if there was more than 7-10 years of data beyond the breakpoint, given it is now a decade on perhaps there is new data available, for example from the USA Death Master File.