Prediction of chronological and biological age from laboratory data

Age is highly predictable from blood laboratory analytes

We trained a random forest model [21] to predict chronological age (in years) using data from 67,563 individuals ranging in age from 1 to 85 years (mean: 36.2, standard deviation: 23.1) from nine CDC National Health and Nutrition Examination Survey (NHANES) cohorts spanning 1999-2016 (Figure 1), a representative sample of the non-institutionalized population of the United States. The model included 356 features consisting of laboratory analytes (e.g. serum glucose, creatinine). Many of the analytes contained a large proportion of missing data (Supplementary Table 3), which was dealt with by imputing missing values using mean imputation. We evaluated model performance both using five-fold cross-validation and held-out data (Methods). Hyperparameters were selected by grid search (Methods). We define our baseline model for chronological age prediction as a linear regression model without regularization, using age as the response variable and the 356 laboratory analytes as covariates. Mean absolute error (MAE) for the baseline linear regression model was 10.53 (SE: 0.07) years in five-fold cross-validation and 10.52 years in the 20% held-out dataset. The R² for the baseline model was 0.63 (SE: 0.01) in the five-fold cross-validation and 0.62 in the held-out set. In our best random forest model, MAE was 4.80 (SE: 0.013) years in cross-validation and 4.76 years in the 20% held- out dataset. The R² from the random forest model was 0.92 (SE: 0.0005) in the five-fold cross-validation and 0.92 in the held-out set.

We also trained separate random forest models for the four main United States Census [22] age groups: [1,18), [18,45), [45,65), 65+. The predictive accuracy of the models differed substantially across age groups (pairwise R² comparisons were significant while adjusting for multiple comparisons; Methods) (Figure 2; Table 1). The model for the [1,18) age group had the Random forest models were trained on data from individuals across different age ranges within the dataset. Mean absolute error and R-squared are shown in the table as measurements of model performance on held-out data.

most accurate predictions of the four and the model trained on the 65+ age group had the least accurate predictions of the four, as measured by R² for age prediction in years. In the held-out dataset, the MAE for [1,18) was 0.87 years and the R² was 0.94. For [18,45), the MAE in the held-out dataset was 5.15 years and the R² value was 0.44; for [45,65), the MAE was 2.20 years and R² value was 0.77; for the 65+ cohort, the MAE was 4.30 years and the R² value was 0.25 (Figure 2).

Feature importance differs substantially across age groups

In order to estimate the predictive power of specific laboratory analytes in a comparable manner, we computed variable importance scores (calculated using the decrease in node impurity using the Gini impurity measure; Methods) for each of the age-specific models across the 356 laboratory analytes. We define the Top-10 set for each age bin as the 10 laboratory analytes with largest variable importance scores for the random forest model trained on that age group, denoted e.g. Top-10_[1,18) for the [1,18) age group (similarly for Top-5). We computed relative variable importance scores for each analyte for each age range, and the total relative importance of the Top-5 and Top-10 sets for each age group, denoted e.g. |Top-10_[1,18)| for the [1,18) age range.

The analytes in the Top-5 differed substantially across age ranges (Figure 3). The only analyte that appears in the Top-5 for multiple age groups is alanine aminotransferase, which appears in both the [1,18) and 65+ groups. |Top-10_[1,18)| was 0.751, |Top-10_[18,45)| was 0.294; |Top-10_[45,65)| was 0.542; |Top-10₆₅₊| was 0.225. |Top-5_[1,18)| was 0.567; |Top-5_[18,45)| was 0.169; |Top-5_[45,65)| was 0.464; |Top-5₆₅₊| was 0.144 (Supplementary Table 1). Top-5_[1,18) consisted of hepatitis B, alkaline phosphatase, lactate dehydrogenase, aspartate aminotransferase, and alanine aminotransferase; Top-5_[18,45) consisted of total cholesterol, serum vitamin E, serum cholesterol, glycohemoglobin, and hepatitis B antibody; Top-5_[45,65) consisted of herpes simplex 1, herpes simplex 2, toxoplasma, HIV 1,2 combo test, and varicella. Top-5₆₅₊ consisted of alanine aminotransferase, blood urea nitrogen, lymphocyte percentage, creatinine, and homocysteine (Figure 3).

Chronological vs. biological age in laboratory data

In the held-out datasets (total n = 13,513), we defined ‘outliers’ as individuals with residual errors in the top 5% or bottom 5% of their age group, representing approximately 37.1 million individuals in the US (based on NHANES sample weights). For the bottom 5% of each age group, the model underestimated their age by an average of 2.20 years (sd = 0.50) for [1,18); 10.9 years (sd = 1.64) for [18,45); 6.38 years (sd = 1.80) for [45,65); and 8.64 years (sd = 1.09) for 65+. For the top 5% of each age group, the model overestimated the age of outliers by an average of 2.47 years (sd = 0.89) for [1,18); 12.23 years (sd = 1.91) for [18,45); 5.52 years (sd = 0.72) for [45,65); and 9.07 years (sd = 1.04) for 65+. We compared the outliers from each age bin with the remaining individuals in the held-out dataset to assess differences in demographic features between these groups (Figure 2). After correcting for multiple comparisons, we found no significant differences between the outlier populations from each age bin and the rest of the individuals from that age bin in gender, race, and income to poverty ratio distributions.

In order to investigate aging in males and females, we stratified 356 analytes individually by age and gender (Figure 2D, Supplementary Figure 2). We found that several analytes, including major blood labs such as red blood cell count, hemoglobin, and hematocrit, and other labs such as alkaline phosphatase, lactate dehydrogenase, and calcium, showed age-related changes that differed starkly between males and females. For hematocrit, hemoglobin, and red blood cell count, male and female values are homogeneously mixed in younger children, separate in the teenage years until male and female values exhibit different ranges, and then cross over again in the senior years. For alkaline phosphatase, this trend is reversed, and sexual dimorphism is apparent in children below 15 years old, and then gradually reduces with age.

Feature scores are highly correlated for males and females of different race/ethnicity groups

We trained separate chronological age predictors for different subpopulations spanning combinations of gender (male, female) and race (Mexican American, Other Hispanic, Non-Hispanic White, Non-Hispanic Black, Other) groups in the NHANES population. We computed correlation coefficients to compare the variable importance scores for each pair of subpopulations (e.g. Mexican American Females and Black Males). Figure 4 shows all pairwise correlations between the feature importance scores across the 10 groups. Pairs of race/ethnicity groups of the same gender all had correlation coefficients above 0.85, with the majority above 0.9. Correlations between feature importance scores of males and females in the same race groups and across race groups were considerably lower, ranging from 0.53-0.80. The strongest correlation between male and female feature importance scores across race groups occurred between white males and white females (0.80).

Models accurate for one age group fail to generalize to other age groups

We tested whether Top-5 and Top-10 for each age group could predict chronological age accurately for other age groups. Table 2 contains the MAE values for models containing only the Top-5 and Top-10 variables when trained and tested on other age groups. The best performing Top-5 model for each age group was the model trained using the Top-5 from that age group. The same was true for the Top-10 models. For the [1,18) age group, the best Top-5 model gave a MAE of 1.35 and the best Top-10 model gave a MAE of 1.16. For the [18,45) age group, the best Top-5 model gave a MAE of 6.41 and the best Top-10 model gave a MAE of 5.51. For the [45,65) age group, the best Top-5 model gave a MAE of 3.28 and the best Top-10 model gave a MAE of 2.91. And for the 65+ age group, the best Top-5 model gave a MAE of 4.63 and the best Top-10 model gave a MAE of 4.49. These results demonstrate that an analyte’s predictive power is not necessarily the same for different age bins and that the set of most predictive analytes is not consistent across age bins.

Widespread non-linearity in analyte aging trajectories

Having observed strong predictive value in the pediatric cohort, we sought to identify significant transitions in biomarker trajectories occurring between the ages of 11 and 30 for comparison against traditional age groupings. To do this, we examined 342 laboratory analytes for piecewise linearity by estimating ‘breakpoints’ in the aging curves (analyte level by age) for ages [11, 30] of each analyte using piecewise regression models [23] (Figure 5). Analytes that did not have data for children younger than 18 years were not included in the analysis. We tested the slopes of the regression lines on either side of the breakpoints for differences. Of the 342 analytes tested, 97 were significant (28.4%) for differences in slope at a Bonferroni-adjusted p-value threshold of 1.46 x 10^-4 (Methods). The median of the 97 breakpoints was 16.4 years with 50% of the breakpoints falling in the range of 15.0-17.7 years. The mode (rounded in years) was 16 years, and the maximum breakpoint was 28.9 years.

In addition to piecewise linearity for the adolescent to adult transition (11 to 30), several different categories of laboratory analyte “aging curves” were observed across the full lifespan, including the following, in order of increasing complexity: (1) linear (e.g. uric acid, iron); (2) piecewise linear (e.g. hematocrit, hemoglobin); (3) power (e.g. alkaline phosphatase, phosphorus); (4) U-shaped (e.g. measles antibody) (Figure 5).

Data

We collected blood laboratory analyte measurements and demographic data from nine waves of the Centers for Disease Control and Prevention (CDC) National Health and Nutrition Examination Survey (NHANES) including the following cohorts: 1999-2000, 2001-2002, 2003-2004, 2005-2006, 2007-2008, 2009-2010, 2011-2012, 2013-2014, 2015-2016. Observations with zero or missing two-year survey weights were removed and all variable names with prefix ‘LBX’ were retained (Supplementary Table 2). Laboratory analytes with greater than 95% missingness were removed (i.e. analytes measured in < 3,901 of individuals), and individuals with fewer than 20 measured labs were also removed (Supplementary Table 3 shows the number of missing values for each laboratory variable). These criteria allowed for analytes to be included in the model with a large proportion of missing values (many contained over 50% missing values). We imputed all missing values using mean imputation, where each missing value was replaced by the mean value for that analyte over all individuals. The final dataset included 67,563 individual and 356 laboratory analytes. Code is available for download here: https://github.com/manrai/Age-Prediction.

Supervised learning for predicting age from laboratory data

We used random forests [21] to train the main age prediction models in this study. Random forests are machine learning models composed of an ensemble of many decision/regression trees. Each of the individual trees in the “forest” is trained on a bootstrap sample of the training data, while the features used for splitting at each node are selected from a random subset of all possible features. Random forests are robust to outliers and perform well on data with linear and non-linear features.

The random forest model used in this analysis was implemented with the scikit-learn library in Python [45]. The data was partitioned randomly using an 80%-20% train-test split, missing values were imputed using mean imputation, and all variables were normalized. We selected hyperparameters using a grid search method in which the maximum number of trees in the random forest and the maximum number of features selected for evaluation at each tree were iteratively evaluated over 50 combinations and scored using five fold cross-validation while taking into account computational time (grid combinations included: max number of trees = 25, 50, 100, 200, 400, 500; max number of features selected = 1, 5, 10, 25, 50, 100; and bootstrap = True, False). We evaluated model accuracy using five-fold cross-validation (scored with the mean absolute error criterion) on the training data and then tested on the 20% held-out dataset. Ordinary least squares linear regression was used for baseline predictions.

Statistical analysis

Variable importance was calculated using Gini impurity, or Mean Decrease Gini, as implemented in the feature_importances_ function in scikit-learn. The importance of a variable X_m is calculated as suggested by Breiman [21, 46, 47]:

${\text{Imp(X}}_{\text{m}})=\frac{1}{N_T}{\displaystyle \sum _T}{\displaystyle \sum _{t\in T:v(s_t)=X_m}p(t)\Delta i(s_t,t)}$

where N_T is the number of trees, v(s_t) is the variable in split s_t, and p(t)Δi(s_t, t) is the weighted impurity decrease (using Gini impurity) for all nodes t where X_m is used [46]. The sum of variable importance scores across all variables is 1:

${\displaystyle \sum _m^M\text{Imp(}{X}_{{}_m})}=1$

For each age range, we defined the total relative importance of the Top-5 and Top-10 variables as:

$$\left|{\text{Top}}_5\right|={\displaystyle \sum _{k\in \text{Top-5}}\text{Imp(}{X}_m{)}_k}$$

$$\left|{\text{Top}}_{10}\right|={\displaystyle \sum _{k\in \text{Top-10}}\text{Imp(}{X}_m{)}_k}$$

Piecewise regression analysis was carried out using the segmented package [48] in R. Piecewise regression models were used to estimate a breakpoint, in this case an age that marks a change in trajectory, for each of the 342 analytes used in the analysis. Breakpoint estimates and corresponding test statistics were computed using the Davies’ test (via the davies.test R function), which tests for non-zero differences in the slope parameter of a segmented relationship between the regression lines on either side of the estimated breakpoint [49]. We used a Bonferroni adjustment to set a threshold for statistical significance at p < 0.05/342, correcting for the 342 analytes.

A vector of feature importance scores was computed for models trained on subgroups consisting of gender and race combinations (10 total subgroups). Pearson product-moment correlations were then computed using pairwise complete importance scores for each subgroup against every other. R² values were computed using predictions in the 20% held-out dataset compared. Chi-squared tests were performed in R using the chisq.test function and the Kolmogorov-Smirnov test was performed using the ks.test function in R 50.