Dataset-Agnostic Biological Age Calculator and Predictions of Mortality in NHANES (1999-2010), UK Biobank (2006-2010), and Dalian Heart Failure (2015-2023) Cohorts: Development and Validation Study

Introduction

As the global population continues to age, it is essential to develop optimized approaches to manage chronic diseases, strengthen health care infrastructure, and promote healthy aging [1], where comprehensive, accurate, and patient-specific estimation of aging status plays a key role. Biological age (BA) refers to the composite physiological state of the organism and/or its various systems [2]. A positive difference between BA and chronological age (CA) indicates accelerated aging. Its estimation is important as BA and (BA – CA) can be used for predicting longevity [3], adverse outcomes [4], and for prognostication [5].

Different models have been proposed to estimate BA based on biomarkers of different mechanisms, including DNA methylation, telomere length changes, proteomics, metabolomics, and the combination of heterogeneous biomarkers [6]. Regarding the algorithm, multiple linear regression and principal component analysis (PCA) methods have been used in processing multiple factors, while the Klemera and Doubal’s method (KDM) provides a more accurate estimation by adjusting the weights of different variables according to their correlations with age [7]. With the recent development of AI, machine learning and deep learning models including deep neural networks, convolutional neural networks, and deep convolutional long-term memory have been introduced to process multidimensional data in BA studies [7,8].

Despite these advances, there are three major limitations that preclude the widespread use of BA models: (1) the high technical demands needed to develop the models [9], (2) lack of alignment of BA models across diverse populations [10], and (3) difficulties in updating the models when new data are received or cohorts are expanded [11]. Thus, the aim of this project is to develop a tool that can address the issues above by simplifying and automating BA model development and updates.

Methods

Program Structure

The present toolkit, named BACalc, facilitates the ingestion of raw clinical datasets stored in CSV format. Multiple CSV files can be loaded and merged row-wise to form a unified dataset. Within the graphical interface, users can select (1) the CA column, (2) optional binary variables (eg, sex, smoking, drinking, etc), and (3) continuous variables (for biomarkers) used for BA estimation. When a binary split variable is defined, the full dataset is partitioned into 2 groups, and the full pipeline is run independently within each group, provided that at least a minimum number of participants are available. In this study, we developed BACalc, a BA calculator that can process raw clinical datasets stored in CSV format. Within the graphical interface, users can select (1) the CA column, (2) binary variables (eg, sex, smoking, and alcohol use), and (3) continuous variables (eg, creatinine levels), used for BA estimation. All numerical operations are implemented via the NumPy and pandas libraries, and linear modeling components are provided by scikit-learn and SciPy. This vectorized implementation allows for rapid processing of typical clinical datasets, while the graphical user interface layers provide an accessible interface for users without the need of scripting or coding. The present toolkit implements 2 methods, PCA method with Dubina correction and the KDM. Both methods operate on the same preprocessed biomarker matrix and CA vector. The program generates a total of three main output files for every user-defined group: (1) the row-wise predictions, (2) the regression statistics, and (3) the explicit linear equations and PCA loadings. Example outputs were generated using the publicly available National Health and Nutrition Examination Survey (NHANES).

The program includes an optional feature to compute mean arterial pressure (MAP) when columns containing systolic blood pressure (SBP) and diastolic blood pressure (DBP) are present. In those cases, MAP is computed according to the standard approximation:

${\displaystyle {\displaystyle \mathrm{M}\mathrm{A}\mathrm{P}=\frac{1}{3}\mathrm{S}\mathrm{B}\mathrm{P}+\frac{2}{3}\mathrm{D}\mathrm{B}\mathrm{P}}}$(1)

It can be included as an additional biomarker in the subsequent modeling steps.

Calculation of BA

The present toolkit implements 2 methods for BA calculation: the PCA method with Dubina correction and the KDM. Both methods operate on the same preprocessed biomarker matrix and CA vector.

Data Preprocessing and Standardization

Overview

Prior to model fitting, all selected variables are converted to numeric form. Any nonnumeric entries are set to missing values. For each biomarker j within a given analysis group, missing entries are imputed using the median of that biomarker in that group. This simple and robust imputation strategy avoids the introduction of bias toward extreme values while retaining all available participants for modeling.

Let x_ij denote the observed value of biomarker j for participant i in a group with N participants and m biomarkers. For each biomarker, the mean μ_j and SD σ_j are computed as follows:

${\displaystyle {\displaystyle {\displaystyle {\mu }_j=\frac{1}{N}\sum _{i=1}^N{x}_{ij},{\sigma }_j=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(x_{ij}-{\mu }_j\right)}^2}}}}$(2)

To avoid division by 0, biomarkers with zero variance are assigned σ_j=1. The standardized (z-scored) biomarker matrix z_ij is then defined as follows:

$z_{ij}=\frac{x_{ij}-{\mu }_j}{{\sigma }_j}$(3)

This standardization step ensures that all biomarkers contribute on a comparable scale, regardless of their original units or ranges.

Method 1: PCA-Dubina Algorithm

The initial part of the computation involves the generation of a biological age score (BAS). As raw biomarker data varied significantly in units and magnitude, the program first performed standardization. PCA is then performed on the standardized biomarker matrix Z, and the first principal component (PC1) is extracted.

Let l = (l₁,...,l_m) denote the loading vector corresponding to PC1. For each participant i, the PC1 score is:

${\displaystyle {\displaystyle {\displaystyle \mathrm{P}\mathrm{C}{1}_i=\sum _{j=1}^m{l}_j{z}_{ij}}}}$(4)

To ensure that higher scores correspond to higher CA, the present implementation automatically corrects the directionality of PC1 by inspecting its Pearson correlation with CA. If corr(PC1,CA) <0, the loading vector I is multiplied by −1 so that PC1 is positively associated with age.

The BAS for participant i was then calculated directly in the original biomarker space as a linear combination of the input biomarkers x_ij weighted by the component loadings normalized by the biomarker SD. Specifically, we defined weights

$w_j=\frac{l_j}{{\sigma }_j}$(5)

and an intercept term w₀ chosen so that the projection is centered:

$w_0=-{\sum }_{j=1}^m{w}_j{\mu }_j$(6)

The governing equation for BAS in the present model is therefore:

${\displaystyle {\displaystyle {\displaystyle {\mathrm{B}\mathrm{A}\mathrm{S}}_i=w_0+\sum _{j=1}^m{w}_j{x}_{ij}}}}$(7)

By construction, BAS is dimensionless. Let μ_BAS and σ_BAS denote the mean and SD of BAS in the training population, and let μ_CA and σ_CA denote the mean and SD of CA in the same group. To convert the dimensionless BAS into a timescale (y), a T-scale transformation is applied:

${\displaystyle {\displaystyle {\mathrm{B}\mathrm{A}}_i=\left({\mathrm{B}\mathrm{A}\mathrm{S}}_i-{\mu }_{\mathrm{B}\mathrm{A}\mathrm{S}}\right)\cdot \left(\frac{{\sigma }_{\mathrm{C}\mathrm{A}}}{{\sigma }_{\mathrm{B}\mathrm{A}\mathrm{S}}}\right)+{\mu }_{\mathrm{C}\mathrm{A}}}}$(8)

This transformation maps the distribution of BA to match the central tendency and spread of CA in the training sample.

However, BA estimates derived from regression or PCA often exhibit a regression to the mean effect, where older individuals appear younger, whereas younger individuals appear older. To rectify this, the Dubina correction was applied [12]. The slope b of the regression of the initial BA against CA is calculated as follows:

${\displaystyle {\displaystyle {\displaystyle b=\frac{\sum _{i=1}^N\left({\mathrm{B}\mathrm{A}}_i-{\mu }_{\mathrm{B}\mathrm{A}}\right)\left({\mathrm{C}\mathrm{A}}_i-{\mu }_{\mathrm{C}\mathrm{A}}\right)}{\sum _{i=1}^N{\left({\mathrm{C}\mathrm{A}}_i-{\mu }_{\mathrm{C}\mathrm{A}}\right)}^2}}}}$(9)

where μ_BA is the mean of the uncorrected BA values in the training group. The corrected biological age, denoted BA_c, was then given by:

${\displaystyle {\displaystyle {\mathrm{B}\mathrm{A}}_{\mathrm{c},i}={\mathrm{B}\mathrm{A}}_i+\left({\mathrm{C}\mathrm{A}}_i-{\mu }_{\mathrm{C}\mathrm{A}}\right)\cdot \left(1-b\right)}}$(10)

This step forced the slope of the predicted age against CA to unity, ensuring that deviations from the diagonal represented genuine biological acceleration or deceleration instead of artifacts of regression shrinkage. For interpretability, the program also reported the BA acceleration derived from the PCA-Dubina pipeline:

${\displaystyle {\displaystyle {\mathrm{A}\mathrm{c}\mathrm{c}\mathrm{e}\mathrm{l}}_{{\mathrm{B}\mathrm{A}}_{\mathrm{c},i}}={\mathrm{B}\mathrm{A}}_{\mathrm{c},i}-{\mathrm{C}\mathrm{A}}_i}}$(11)

This quantifies how many years older (positive values) or younger (negative values) a participant appears biologically relative to their CA.

Method 2: KDM

The second module uses the KDM, which assumes that BA is a latent variable driving the linear deterioration of biomarkers. For each biomarker j, a simple linear model is fitted with CA as the predictor:

${\displaystyle {\displaystyle x_{ij}=k_j\cdot {\mathrm{C}\mathrm{A}}_i+q_j+{ϵ}_{ij}}}$(12)

where k_j is the slope, q_j is the intercept, and ∈_ij are residuals. These parameters are estimated via ordinary least squares regression. The root mean square error (RMSE) for biomarker j is given by

${\displaystyle {\displaystyle {\displaystyle s_j=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(x_{ij}-{\hat{x}}_{ij}\right)}^2}}}}$(13)

where ${\hat{x}}_{ij}=k_j{\mathrm{C}\mathrm{A}}_i+q_j$ is the model prediction. In the implementation, s_j is constrained to a small positive minimum to avoid division by 0. For each biomarker, the Pearson correlation coefficient r_j between x_j and CA is determined.

The algorithm first computes the BA estimate (BA_E) for each participant, which aggregates the contributions of all biomarkers. A weight w_j is assigned to each biomarker, defined as the square of the signal-to-noise ratio:

$w_j={\left(\frac{k_j}{s_j}\right)}^2$(14)

For participant i, the KDM BA estimate BA_E,i was then computed as the weighted average of the ages predicted by each individual biomarker:

${\displaystyle {\displaystyle {\displaystyle {\mathrm{B}\mathrm{A}}_{\mathrm{E},i}=\frac{\sum _{j=1}^m{w}_j\left(\frac{x_{ij}-q_j}{k_j}\right)}{\sum _{j=1}^m{w}_j}}}}$(15)

In addition, our program implements a variance-based adjustment to calculate the corrected biological age (BA_EC). This requires the estimation of the characteristic correlation (r_char) of the biomarker set. Based on the individual correlations r_j, the characteristic correlation is computed by:

${\displaystyle {\displaystyle {\displaystyle r_{\text{char}}=\frac{\sum _{j=1}^m\frac{r_j^2}{\sqrt{1-r_j^2}}}{\sum _{j=1}^m\frac{r_j}{\sqrt{1-r_j^2}}}}}}$(16)

with numerical safeguards applied when |r_j| is extremely small or close to unity.

This metric was then used to estimate the population variance of BA, $S_{\mathrm{B}\mathrm{A}}^2$. In the implementation, 2 components are considered. First, the variance of the difference between the KDM estimate and CA is computed:

${\displaystyle {\displaystyle {\mathrm{Var}}_{\mathrm{\Delta }}=\mathrm{Var}\left({\mathrm{B}\mathrm{A}}_E-\mathrm{C}\mathrm{A}\right)}}$(17)

Second, the variance of CA is approximated from its range:

${\displaystyle {\displaystyle \widehat{\mathrm{V}\mathrm{a}\mathrm{r}}\left(\mathrm{C}\mathrm{A}\right)\approx \frac{{\left(\underset{i}{max}{\mathrm{C}\mathrm{A}}_i-\underset{i}{min}{\mathrm{C}\mathrm{A}}_i\right)}^2}{12}}}$(18)

Following Klemera and Doubal, a correction term that depends on r_char and the number of biomarkers m is then subtracted:

${\displaystyle {\displaystyle S_{\mathrm{B}\mathrm{A}}^2=max\left(S_{min}^2,{\mathrm{Var}}_{\mathrm{\Delta }}-\frac{1-r_{\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}}^2}{r_{\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}}^2}\cdot \frac{\widehat{\mathrm{Var}}\left(\mathrm{C}\mathrm{A}\right)}{m}\right)}}$(19)

where $S_{\mathrm{m}\mathrm{i}\mathrm{n}}^2$ is a small positive floor (set to 0.1 in the present implementation) to avoid unrealistically small or negative variance estimates.

Finally, for each participant i, the corrected KDM BA BA_EC,i is computed by combining the biomarker-based BA estimate (BA_E,i) with CA (CA_i) in a Bayesian fashion:

${\displaystyle {\displaystyle {\displaystyle {\mathrm{B}\mathrm{A}}_{\mathrm{E}\mathrm{C},i}=\frac{{\mathrm{B}\mathrm{A}}_{\mathrm{E},i}\cdot \left(\sum _{j=1}^m{w}_j\right)+\left(\frac{{\mathrm{C}\mathrm{A}}_i}{S_{\mathrm{B}\mathrm{A}}^2}\right)}{\left(\sum _{j=1}^m{w}_j\right)+\left(\frac{1}{S_{\mathrm{B}\mathrm{A}}^2}\right)}}}}$(20)

This formulation effectively reduces the prediction error by anchoring the estimate to the CA when the biomarker data are noisy or incoherent, while still allowing large deviations when the biomarker pattern strongly supports them.

Analogous to the PCA-Dubina pipeline, the KDM module reports a BA acceleration measure based on the corrected KDM estimate:

${\displaystyle {\displaystyle \text{Acce}{\text{l}}_{{\mathrm{B}\mathrm{A}}_{\text{EC,i}}}={\mathrm{B}\mathrm{A}}_{\mathrm{E}\mathrm{C},i}-{\mathrm{C}\mathrm{A}}_i}}$(21)

Program Output and Computational Scheme

The program generates three main output files for each user-defined group: (1) the row-wise predictions, (2) the regression statistics, and (3) the explicit linear equations and PCA loadings. The predictions file stores, for each participant, the CA, the intermediate BAS, the PCA-derived BA and BA_c, the KDM-derived BA_E and BA_EC, and the corresponding acceleration measured (BA_c – CA and BA_EC – CA). The statistics file contained global goodness-of-fit metrics for each method (PCA-Dubina and KDM), including the coefficient of determination (R²), Pearson correlation coefficient (r), RMSE, and mean absolute error (MAE) between each BA estimate and CA.

For a given method, let y_i denote the CA of participant i and ${\hat{y}}_i$ the corresponding estimated BA (eg, BA_c,i or BA_EC,i). The coefficient of determination was calculated using:

${\displaystyle {\displaystyle {\displaystyle R^2=1-\frac{\sum _{i=1}^N{\left(y_i-{\hat{y}}_i\right)}^2}{\sum _{i=1}^N{\left(y_i-\overline{y}\right)}^2}}}}$(22)

where $\overline{y}$ is the mean of y_i. The RMSE and MAE are given by:

${\displaystyle {\displaystyle {\displaystyle \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(y_i-{\hat{y}}_i\right)}^2}}}}$(23)

${\displaystyle {\displaystyle {\displaystyle \mathrm{M}\mathrm{A}\mathrm{E}=\frac{1}{N}\sum _{i=1}^N\left|y_i-{\hat{y}}_i\right|}}}$(24)

Pearson correlation coefficient between y_i and ${\hat{y}}_i$ is reported:

${\displaystyle {\displaystyle {\displaystyle r=\frac{\sum _{i=1}^N\left(y_i-\overline{y}\right)\left({\hat{y}}_i-\overline{\hat{y}}\right)}{\sqrt{\sum _{i=1}^N{\left(y_i-\overline{y}\right)}^2\sqrt{\sum _{i=1}^N{\left({\hat{y}}_i-\overline{\hat{y}}\right)}^2}}}}}}$(25)

with $\overline{\hat{y}}$ being the mean of the estimated ages.

The generation of explicit equations allows users to apply the trained models externally without rerunning the optimization pipeline, which is useful in the development of portable apps on mobile devices where computational resources are rather limited. For each group, the program writes a human-readable summary that includes the PCA-based BAS and BA equations, the Dubina-corrected BA_c equation as a function of CA and biomarker values, and the full set of KDM parameters (k_j, q_j, s_j, and w_j) together with the final BA_E and BA_EC formulas. In addition to tabular outputs, the program generated plots for graphical illustration. For each group, scatterplots of BA_c vs CA and BA_EC vs CA are generated, together with identity lines and fitted regression lines with 95% confidence bands.

Beyond predictive accuracy, the computational efficiency of BACalc was a primary design priority. By leveraging the Python NumPy library for high-performance vectorized operations and implementing compact linear models, BACalc achieves significant speed advantages. Even when processing extensive clinical datasets (n>10,000), the system is capable of training and evaluation in a matter of seconds on standard hardware, making it highly scalable for large-scale epidemiological studies. The BACalc source code and executable can be downloaded freely from [13].

Model Validation and Performance Evaluation

To validate the robustness and accuracy of the proposed model, we used the NHANES dataset as curated by the BioAge R package, with no additional inclusion or exclusion criteria [14]. This validation process incorporated a comprehensive panel of 13 physiological biomarkers: albumin, alkaline phosphatase, blood urea nitrogen, creatinine, hemoglobin A_1c (HbA_1c), C-reactive protein (CRP), lymphocytes, mean corpuscular volume, glucose, red cell distribution width, total cholesterol, white blood cell (WBC) count, and SBP. Using this standardized feature set, we benchmarked the BA estimates generated by our implementation of the KDM and the PCA-Dubina method against the established KDM projections from the BioAge R package. Furthermore, we have ensured median imputation in the BioAge R package so that it matches the default median imputer implemented in BACalc. To empirically evaluate the impact of imputation choice, we performed a sensitivity analysis on the NHANES dataset using 4 commonly used imputation techniques: median, mean, k-nearest neighbors (k=5), and iterative imputation (MICE [multiple imputation by chained equations]) [15]. The NHANES dataset comprised 55,081 participants, taken from the BioAge R package repository.

To further validate our tool, we conducted a comprehensive Bland-Altman agreement and residual analysis against the established BioAge R package, applying both tools to the same unweighted cohort to ensure a direct, one-to-one methodological comparison. We compared the BACalc KDM implementation (BA_EC) against the standard KDM in the R package, and similarly compared the BACalc PCA-Dubina model (BAc) against the standard. To analyze the sensitivity of the proposed method to age and sex, we then evaluated the calibration across distinct age bins (<40, 40‐60, and >60 y) with sub-group analysis between genders.

Additional Datasets for Predictive Validity Evaluation

To validate the clinical relevance of BACalc estimates, we assessed the association between BA acceleration (delta age=BA-CA) and longitudinal clinical outcomes in 2 independent datasets.

We used a retrospective population-based clinical cohort comprising 2781 patients hospitalized with heart failure from the First Affiliated Hospital of Dalian Medical University (2015‐2023). To comprehensively capture multisystem aging, BACalc estimates were generated using a 20-biomarker clinical panel. This included markers of cardiovascular/hemodynamic stress (SBP, DBP, B-type natriuretic, and high-sensitivity cardiac troponin I), renal function (creatinine, urea, and uric acid), metabolic/lipid profiles (glucose, total cholesterol, triglycerides, low-density lipoprotein-C, and high-density lipoprotein-C), liver function (alanine aminotransferase, aspartate aminotransferase, and total bilirubin), and hematological/immune status (WBC, neutrophils, hemoglobin, red blood cells [RBC], and platelets). We used Cox proportional hazards models, controlling for CA at baseline, sex, smoking history, and drinking history, with administrative right-censoring applied to surviving patients. We calculated BA acceleration and linked these estimates to longitudinal clinical outcomes using fully adjusted Cox proportional hazards models and Kaplan-Meier survival analyses.

In addition, we evaluated the predictive validity on the large-scale UK Biobank (UKBB) dataset. We included 501,936 participants with available biomarker data at the initial assessment center visit (assessment instance 0) and linked mortality follow-up. The biomarkers include renal function (creatinine, urea, urate, and cystatin C), metabolic/lipid profiles (glucose, total cholesterol, triglycerides, low-density lipoprotein [direct], high-density lipoprotein-C, and HbA_1c), liver function (alanine aminotransferase, aspartate aminotransferase, alkaline phosphatase, CRP, albumin, total bilirubin, and total protein), and hematological/immune status (WBC, hemoglobin, RBC, nucleated RBC, platelets, monocytes, neutrophils, lymphocytes, eosinophil, hematocrit percentage, mean corpuscular volume, and mean corpuscular hemoglobin concentration). BA acceleration was calculated using 29 blood count and biochemistry biomarkers measured at this visit, and 54,507 all-cause deaths were observed during follow-ups.

Ethical Considerations

This study used existing deidentified and anonymized data in the public domain directly downloadable from the NHANES website. Thus, according to the Common Rule, this study was exempt from the institutional review board review and the informed consent requirement. NHANES was conducted by the Centers for Disease Control and Prevention National Center for Health Statistics. NHANES survey procedures and protocols were approved by the National Center for Health Statistics ethics review board for each survey cycle [16].

The UKBB analysis was performed under application number 79146 as part of a wider study for predictive modeling for cardiovascular disease. The UKBB was approved by the North West Multi-Centre Research Ethics Committee (number 11/NW/0382) and conducted in accordance with the Declaration of Helsinki.

For the dataset from the First Affiliated Hospital of Dalian Medical University, this study was a retrospective analysis of existing clinical data and did not involve any interventions or interactions with participants. The study was conducted in accordance with the Declaration of Helsinki and was approved by the Ethics Committee of the First Affiliated Hospital of Dalian Medical University (approval number PJ-KS-KY-2024‐616). As this was a retrospective study using anonymized data, informed consent was waived by the Ethics Committee.

Results

Initial Validation of the BACalc Model

Benchmarking studies demonstrate that the BACalc model generates BA estimations in close alignment with the BioAge R package (Figure 1). This high level of congruency validates the algorithmic integrity of the BACalc implementation. Specifically, when evaluating the KDM model, the BACalc estimations exhibit a near-perfect correlation with the BioAge R package results (r=0.9937), accompanied by an MAE of only 2.02 years. Furthermore, the PCA-Dubina model within the BACalc framework shows strong statistical agreement with the R package reference, yielding a correlation coefficient of 0.8871. When analyzing the variance across these methodologies, the data distribution of the BACalc KDM model more closely mirrors that of the BioAge R package than the PCA-based approach. Notably, the accelerated age residuals representing the divergence between BA and CA remained largely within a ±20-year window for both the BioAge R package and the BACalc PCA-Dubina model, with the highest data density concentrated within this range.

‎

Sensitivity Analysis of Imputation Methods

This dataset reflects a realistic epidemiological scenario with an overall missingness of 9.79% (ranging from ~9% for blood pressure to 29.69% for CRP). As shown in Figure 2, BA predictions remained highly robust regardless of the imputation method used. The predictions across all 4 methods exhibited near-perfect correlation (Pearson r>0.991 for all pairwise comparisons). The density distributions of predicted ages (BAc and BA_EC) and the variance of age acceleration were visually and statistically indistinguishable across the median, mean, k-nearest neighbors, and MICE techniques.

‎

KDM vs PCA-Dubina: Comparison With CA

The comparative performance of the 2 methods, PCA-Dubina and KDM, implemented in BACalc is shown in Figure 3. When evaluated against CA, the KDM approach demonstrated superior predictive precision, achieving a higher coefficient of determination (R²=0.597) and correlation coefficient (r=0.844) compared to the PCA method (R²=0.174; r=0.740).

‎

There is a strong, positive correlation between the 2 BA estimates (R²=0.803; r=0.896). The Bland-Altman agreement analysis further elucidates this relationship, showing a mean difference of 0.00 years, indicating no systematically fixed bias between the 2 estimators. However, the limits of agreement are notably wide (±21.68 y), implying substantial random variation in individual predictions. Furthermore, the discrepancy between the 2 methods tends to increase as the average BA increases.

KDM vs PCA-Dubina: Comparison With Standard KDM

The Bland-Altman analysis demonstrated exceptional calibration between the BACal and standard KDM models, yielding an overall mean bias of 0.000 years and highly constrained limits of agreement ([−10.68, 10.68] y; Figure 4A). The overall MAE was 4.14 years, with an RMSE of 5.45 years, confirming high precision between the 2 models. As shown in Figure 4C, the residuals are normally distributed and centered around 0.

Between BACalc PCA-Dubina model (BAc) and the standard KDM, the limits of agreement widen to [−21.51, 21.51] years (Figure 4B), with an RMSE of 10.97 years. This illustrates the fundamental mathematical divergence between a variance-based paradigm (PCA) and a regression-based paradigm (KDM) when assessing physiological decline. Crucially, despite this expected individual-level variance, the 2 algorithms maintain an overall mean bias of 0.000 years.

Regarding the calibration across distinct age ranges (<40, 40‐60, and >60 y). For both methods (KDM and PCA in BACalc), the mean bias against standard KDM in BA estimation remained near-absolute zero across all age groups (eg, maximum age-group bias for PCA vs KDM was merely 0.42 y). Subgroup analysis by gender revealed only a marginal mean difference (males: −2.05 y; females: +1.90 y). These robust agreement metrics confirm that BACalc replicates biological aging methodologies with high mathematical fidelity and remains highly stable across all demographic subgroups without systematic drift.

‎

Explicit Linear Equations

The explicit linear equations derived from the BACalc model provide a quantitative framework for assessing biological aging through the integration of 13 key biomarkers. Equation (26) illustrates the multivariate regression model. Notably, biomarkers associated with metabolic dysregulation (HbA_1c and glucose) and inflammation (CRP) exhibit positive coefficients, contributing to an acceleration of the BA estimate. Conversely, markers of synthetic liver function, such as albumin, display negative coefficients, consistent with the age-associated decline in hepatic protein synthesis.

Table 1 further decomposes these relationships using the KDM parametrization. The specific weights assigned to each biomarker reflect their relative reliability in predicting CA within the NHANES cohort. SBP and blood urea nitrogen emerged as highly weighted parameters, suggesting that vascular and renal stiffness are most associated with aging in this specific model.

In contrast, although inflammatory markers, such as CRP, show a steep slope (ie, k), they are assigned lower weights, likely due to higher interindividual variance. Finally, equation 28 applies a corrective transformation to the calculated BA (ie, BA_E). This step integrates a CA anchor to mitigate the characteristic regression to the mean bias often observed in biological clocks, thereby improving prediction accuracy at the extremes of the age range.

${\displaystyle {\displaystyle \begin{array}{rl}{\mathrm{B}\mathrm{A}}_{\mathrm{c}}=& -60.999\\ & +(0.58257\times \mathrm{C}\mathrm{A})\\ & +(-9.8184\times \mathrm{a}\mathrm{l}\mathrm{b}\mathrm{u}\mathrm{m}\mathrm{i}\mathrm{n})\\ & +(0.15875\times \mathrm{a}\mathrm{l}\mathrm{k}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{s}\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{s}\mathrm{e})\\ & +(0.61715\times \mathrm{b}\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n})\\ & +(9.804\times \mathrm{c}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{e})\\ & +(6.6414\times \mathrm{h}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{i}\mathrm{n}{\mathrm{A}}_{1\mathrm{c}})\\ & +(14.637\times \mathrm{l}\mathrm{o}\mathrm{g}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{d}\mathrm{C}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n})\\ & +(-0.27053\times \mathrm{l}\mathrm{y}\mathrm{m}\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{e}\mathrm{s})\\ & +(-0.34334\times \mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{p}\mathrm{u}\mathrm{s}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{r}\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e})\\ & +(0.18261\times \mathrm{g}\mathrm{l}\mathrm{u}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{e})\\ & +(3.3212\times \mathrm{r}\mathrm{e}\mathrm{d}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{w}\mathrm{i}\mathrm{d}\mathrm{t}\mathrm{h})\\ & +(0.0032436\times \mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{l})\\ & +(0.17593\times \mathrm{s}\mathrm{y}\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{c}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{d}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e})\\ & +(1.1794\times \mathrm{w}\mathrm{h}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{d}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t})\end{array}}}$(26)

${\displaystyle {\displaystyle {\displaystyle \mathrm{B}{\mathrm{A}}_{\mathrm{E}}=\frac{\sum \left({\mathrm{W}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}_i\times \frac{{\mathrm{B}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{r}\mathrm{k}\mathrm{e}\mathrm{r}}_i-q_i}{k_i}\right)}{\sum \left({\mathrm{W}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}_i\right)}}}}$(27)

${\displaystyle {\displaystyle {\mathrm{B}\mathrm{A}}_{\mathrm{E}\mathrm{C}}=\frac{(\mathrm{B}{\mathrm{A}}_{\mathrm{E}}\times 0.0024406)+\left(\frac{\mathrm{C}\mathrm{A}}{397.1}\right)}{0.0024406+\left(\frac{1}{397.1}\right)}}}$(28)

The BACalc model generated explicit linear equations that successfully capture the multisystem nature of biological aging. The model parameterization (Table 1) reveals that renal function and vascular health are the most heavily weighted predictors of age in the NHANES dataset. The resulting algorithm integrates these weighted physiological inputs with a CA correction to produce a BA estimate that is both physiologically interpretable and mathematically stable despite variance in individual biomarkers.

Predictive Validity of BACalc Estimates for Clinical Outcomes

Both models generated by BACalc demonstrated highly significant predictive validity for all-cause mortality using 2 additional datasets. Thus, on the UKBB dataset, both biological age acceleration measures generated by BACalc were significantly associated with all-cause mortality in unadjusted Cox models in the UKBB. A 1-year increase in PCA-Dubina–derived age acceleration was associated with a 1.77% higher hazard of all-cause mortality (hazard ratio [HR] 1.018, 95% CI 1.017‐1.018; P<.001). Similarly, KDM-derived age acceleration was associated with a 1.59% higher hazard of mortality per year of accelerated aging (HR 1.016, 95% CI 1.016‐1.016; P<.001; Figure 5).

‎

In the Dalian heart failure cohort, a 1-year increase in age acceleration derived from the PCA-Dubina model (BAc) independently predicted mortality with a 2.84% increased hazard per year (HR 1.028, 95% CI 1.022-1.035; P<.001). Similarly, KDM-derived age acceleration (BA_EC) was associated with a 1.26% increased hazard of mortality per year of accelerated aging (HR 1.013, 95% CI 1.007-1.018; P<.001). Kaplan-Meier survival analysis demonstrated clear risk stratification for mortality based on BACalc estimations. Across both the PCA-Dubina and KDM models, participants in the highest quartile of age acceleration (Q4, fastest aging) exhibited significantly lower overall survival probabilities compared to those in the lowest quartile (Q1, slowest aging; Figure 6). These distinct survival trajectories for all-cause mortality were highly statistically significant (log-rank P<.001), visually underscoring the prognostic capability of the toolkit.

‎

On the UKBB dataset, both BA acceleration measures generated by BACalc were significantly associated with all-cause mortality in unadjusted Cox models. A 1-year increase in PCA-Dubina–derived age acceleration was associated with a 1.77% higher hazard of all-cause mortality (HR 1.018, 95% CI 1.017‐1.018; P<.001). Similarly, KDM-derived age acceleration was associated with a 1.59% higher hazard of mortality per year of accelerated aging (HR 1.016, 95% CI 1.016‐1.016; P<.001; Figure 6).

Discussion

KDM vs PCA-Dubina

The BACalc model generates BA estimations in close alignment with the BioAge R package. When evaluating the KDM model, the BACalc estimates exhibit a near-perfect correlation with the BioAge R package results (r=0.9937), accompanied by an MAE of only 2.02 years. The PCA-Dubina model within the BACalc framework also shows strong statistical agreement with the R package reference, yielding a correlation coefficient of 0.8871. The accelerated age residuals, representing the divergence between BA and CA, remained largely within a ±20-year window for both the BioAge R package and the BACalc PCA-Dubina model, with the highest data density concentrated within this range.

While both KDM and PCA-Dubina methods show a positive linear relationship with CA, the PCA-Dubina model exhibits significantly higher variance and dispersion around the regression line, particularly in older age groups (Figure 3A-D), suggesting less stability in its predictions relative to the KDM model.

The distribution of age acceleration highlights distinct characteristics in how these models quantify aging. The PCA method produces a platykurtic distribution with a wider spread of acceleration values compared to the more leptokurtic distribution of the KDM (Figure 3E). Importantly, the acceleration bias check reveals systematic bias in the PCA method; as CA increases, the PCA-derived age acceleration shows a positive upward trend. In contrast, the KDM model displays a consistent, flat distribution across the age spectrum, indicating that it has reduced residual bias and likely serves as a more robust metric for biological aging independent of time.

PCA models are simple and straightforward and can therefore be used for screening aging markers or solving covariance problems, whereas the assumption of linearity limits the accuracy [17]. In comparison, the KDM models improve accuracy by adjusting the weights of variances. Although KDM models have a higher computational burden, they are more accurate and stable with less dispersion and distortion at either end of the age spectrum [18]. Based on our results and the conclusions of existing studies, we recommend that KDM models are preferable to PCA ones especially when BA has a wide spectrum and heterogeneous biomarkers are included.

Regarding the predictive validity, the aging accelerations derived from both KDM models showed high performance on the datasets of patients with cardiovascular diseases and normal population. The BA-related biomarkers in these datasets cover the liver, renal, and immune functions and comprehensively reflect the health status of the participants [19]. The survival curves on the Dalian dataset of patients who are hospitalized differ from the UKBB ones in shape. The former indicates that the mortality rate was high since the observation, while the latter reflects the common trend that deaths happen at higher ages in the general population. Therefore, both KDM models show the potential to predict long-term health status in both patients and healthy participants. An earlier study on a smaller UKBB subset showed that the biomarkers in KDM-based BA models significantly improved the prediction of mortality and hospital admissions over CA alone [20]. A recent study on 2 large Chinese cohorts showed that KDM-based BA could predict mortality independently from CA and other potential confounders [5]. Therefore, by lowering the technical barrier for BA estimation, our proposed tool may provide a potential tool of not only BA estimation but also clinical risk prediction, particularly for health researchers and public health practitioners who may not be familiar with R-based tools.

Analysis of Aging Factors: Neglected Long-Term Health Issues

Some well-established biomarkers such as blood urea nitrogen and creatinine are reflected in our results [21], while their relationships with aging might be nonlinear [22]. Meanwhile, the results revealed some aging-related biomarkers that are often neglected, including HbA_1c, mean corpuscular volume, and glucose. HbA_1c and glucose levels reflect the risk of metabolic disorder and diabetes, which is a lasting and increasing public health problem especially in lower-income countries [23], where the pathological changes during prediabetes and multiple complications consist of a major challenge in the aging population. Mean corpuscular volume has recently been proposed as a predictive marker for colorectal or esophageal cancer, chronic renal disease, and ischemic stroke or coronary procedures [24]. Further investigation on these factors in different cohorts can indicate the aging acceleration for different organs and systems.

Limitations and Future Directions

Despite the good performance of BACalc, there are several methodological constraints: (1) the current implementation is limited to cross-sectional datasets and relies on linear imputation using the median for missing values, which may introduce bias when missingness is systematic or when biomarkers exhibit skewed distributions; (2) both the PCA-Dubina and KDM pipelines assume linear relationships between individual biomarkers and CA, which can obscure biologically relevant nonlinear or interaction effects; and (3) the present model currently supports only a single modal set of blood-based chemistry variables, restricting its applicability to cohorts that have incorporated emerging omics layers (eg, DNA methylation [25], proteomics [26], and metabolomics [27], or imaging-derived phenotypes [28]).

Future development will address these gaps by extending BACalc to longitudinal analyses, enabling the estimation of individual aging trajectories and the computation of rate of change metrics that are less susceptible to regression-to-the-mean artifacts. Incorporating flexible nonlinear modeling frameworks such as generalized additive models [29], penalized splines [30], machine learning, or neural networks [31] will allow the capture of curvature and interaction effects while preserving interpretability. We also plan to integrate a modular multiomics interface, wherein users can supply heterogeneous biomarker matrices and BACalc will automatically harmonize scales, perform feature selection, and combine modalities through stacked ensemble learning [32]. Finally, to improve robustness to missingness, advanced imputation strategies (MICE and matrix completion) and uncertainty propagation will be implemented.

Conclusions

We developed an open-source and versatile dataset-agnostic model named BACalc; this model employs PCA-Dubina and KDM models for predicting BA from biomarkers and CA. We have validated our model using previous methods and obtained good agreement. BACalc has the ability to accurately estimate BA from biomarkers, produce detailed plots, and provide explicit linear equations for a specific user-supplied cohort. The present toolkit and its open-source nature will be useful in reducing discrepancies in BA estimation through a transparent and reproducible framework.