Maximizing insights from longitudinal epigenetic age data: simulations, applications, and practical guidance

Abstract

Background

Epigenetic age (EA) is an age estimate, developed using DNA methylation (DNAm) states of selected CpG sites across the genome. Although EA and chronological age are highly correlated, EA may not increase uniformly with time. Departures, known as epigenetic age acceleration (EAA), are common and have been linked to various traits and future disease risk. Limited by available data, most studies investigating these relationships have been cross-sectional, using a single EA measurement. However, the recent growth in longitudinal DNAm studies has led to analyses of associations with EA over time. These studies differ in (1) their choice of model; (2) the primary outcome (EA vs. EAA); and (3) in their use of chronological age or age-independent time variables to account for the temporal dynamic. We evaluated the robustness of each approach using simulations and tested our results in two real-world examples, using biological sex and birthweight as predictors of longitudinal EA.

Results

Our simulations showed most accurate effect sizes in a linear mixed model or generalized estimating equation, using chronological age as the time variable. The use of EA versus EAA as an outcome did not strongly impact estimates. Applying the optimal model in real-world data uncovered advanced GrimAge in individuals assigned male at birth that decelerates over time.

Conclusion

Our results can serve as a guide for forthcoming longitudinal EA studies, aiding in methodological decisions that may determine whether an association is accurately estimated, overestimated, or potentially overlooked.

Background

Chronological age, the passage of time since birth, does not fully capture an individual’s state or pace of biological aging [1]. Genetics, along with various environments, behaviors, and diseases faced throughout the life course appear to be potent causes of these disparities. While the exact biochemical mechanisms mediating these effects remain largely unknown, there is emerging evidence of a role for epigenetics, a biological process that can induce changes in gene expression without changing the underlying DNA sequence [2]. The adaptable nature of epigenetic modifications is therefore of utmost interest when evaluating the effects that exposures have on health and lifespan.

DNA methylation (DNAm) is the most studied and easiest-to-measure epigenetic modification. DNAm most commonly occurs at cytosine nucleotides that are followed by guanine, known as CpG sites. A decade ago, Horvath [3] as well as Hannum and colleagues [4] introduced algorithms, also known as epigenetic clocks, that identified a set of CpG sites whose methylation state can be used to accurately estimate chronological age. The resulting measure is called epigenetic age (EA). Additional epigenetic clocks have since been published; some aim to best estimate chronological age [5], while others focus on health and mortality [6,7,8,9,10]. Over the past decade, the scientific community has thoroughly investigated EA-related associations, largely indicating relations between advanced EA and adverse health outcomes [11,12,13,14].

In most of those studies, EA was investigated in a cross-sectional setting. However, the increasing accessibility of longitudinal epigenetic cohort data [15] has created growing interest in studying EA over time [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44]. Unfortunately, there are many disparities in modeling strategies across these studies. Our study investigated whether these disparities impacted findings to an extent that might lead to false conclusions, through significantly inflating effects or leaving true associations unnoticed. We evaluated the robustness of methods using simulations. To test our results in real-world data, we applied the same methods in two examples from the Avon Longitudinal Study of Parents and Children (ALSPAC) [45, 46] involving biological sex and birthweight as predictors of longitudinal EA. Our aim was to provide readers with practical guidance in modeling choices that fit their data and maximize insights in epidemiological relations.

Methods

Modeling longitudinal epigenetic age

Numerous approaches exist to model effects that exposures have on longitudinal outcomes like EA. Approaches differ in (1) the choice of model, (2) outcome, and (3) time variables, as well as (4) the number of repeated measures included in those methods.

The three most frequently applied models are linear mixed effect models (LME) [16,17,18,19,20, 23, 25, 30, 31, 33,34,35,36, 38, 41], generalized estimating equations (GEE) [21, 26, 29, 43], and Δ aging [22, 27, 37, 39, 40, 42, 47, 48]. LME models and GEE both analyze repeated measures, like tracking changes in EA over time in an individual’s life, and model differences in mean trends between groups, such as those exposed to a factor compared to those who were not. In both methods, the fixed effect ($\beta_{2}$), which stays consistent throughout all measures, as well as the interaction effect ($\beta_{4}$), which accumulates over time, are typically modeled as:

$$\begin{aligned} Outcome_{{{\text{ij}}}} & = \beta_{0i} + \beta_{1} Time Variable_{ij} + \beta_{2} Exposure_{i} + \beta_{3} Sex_{i} \\ & \quad + \beta_{4} \left( {Time Variable_{ij} \times Exposure_{i} } \right) + \epsilon_{ij} \\ \end{aligned}$$

(1)

where $Outcome_{ij}$ is either EA or epigenetic age acceleration (EAA, the residual resulting from regressing EA on chronological age) measured in individual i, at timepoint j. To accommodate longitudinal changes, both models account for time. Time variables commonly used, summarized in Fig. 1, range from chronological age at the time of measurement to age-independent variables, including the duration in days or years between measurements, numerical ranks (e.g., 1, 2, 3) or factorized values (e.g., F07, F09, F15).

Fig. 1

Overview of common time variables included in longitudinal epigenetic age studies and their application in ARIES cohort data. Age variables in our example should be interpreted as a single study participant’s age at measurement, since chronological age differs between individuals

Another method to analyze variations in temporal changes is a two-step approach involving “Δ aging” (“delta aging”). The Δ aging method is limited to studies with two repeated measures, as it quantifies the difference between measurements, often referred to as the “Δ aging” score, and then compares these scores between groups, for example using linear regression. Firstly, Δ aging is typically calculated as the difference between a follow-up and a baseline measure of either EA or EAA, with or without adjustment for the duration of time between measurements:

$$\begin{aligned} & \Delta\, Aging = Outcome_{2} - Outcome_{1} \\ & \Delta\, Aging \left( {\text{age adjusted}} \right) = \frac{{Outcome_{2} - Outcome_{1} }}{{Age_{2} - Age_{1} }} \\ \end{aligned}$$

(2)

where $Outcome_{1}$ and $Age_{1}$ represent EAA or EA and age at the initial measure, while $Outcome_{2}$ and $Age_{2}$ correspond to the follow-up measure. Secondly, models such as linear regression can be used to compare trends in Δ aging between the different groups:

$$\Delta Aging_{{\text{i}}} = {\upbeta }_{0} + \beta_{1} Exposure_{i} + \beta_{2} Sex_{i} + \epsilon _{i}$$

(3)

where $\Delta Aging_{{\text{i}}}$ is the above-described difference score between two measures for individual i.

Study population

This study used longitudinal DNAm data generated as part of the Avon Longitudinal Study of Parents and Children (ALSPAC) [45, 46]. Initially, ALSPAC recruited 14,541 pregnant women, resident in Avon, UK, with expected dates of delivery between April 1991 and December 1992. Of the initial 14,541 pregnancies, 14,062 resulted in live births and 13,988 children were alive at 1 year of age. To bolster the initial sample, 1000 additional children were included after the initial participants were approximately 7 years old, increasing the total sample of data collected after age 7–14,901. As part of the Accessible Resource for Integrated Epigenomic Studies (ARIES) [15], a sub-sample of ALSPAC mother–child pairs have undergone genome-wide DNAm analysis. DNAm wet-lab and pre-processing analyses were performed at the University of Bristol as part of the ARIES project [15, 49]. Here, we included up to three within-person DNAm measures, drawn from blood, at age 7, 9, and 15 or 17. A summary of all variables used in our analyses is presented in Table 1. Out of 1,162 individuals, 942 have DNAm measured for at least two of these time points, and 178 individuals have DNAm measured at all three.

Table 1 Characteristics of the ARIES sample across all measurement waves.

Simulation study

To better understand the extent to which methodological choices influence the robustness of results, we compared commonly used methods (introduced under Modeling Longitudinal Epigenetic Age). We conducted a series of simulations (n = 1000), in which we manipulated longitudinal EA data from the ALSPAC cohort, applied all models introduced above, and compared the accuracy of effect estimates across methods (Fig. 2).

Fig. 2

Overview of simulation study. Simulations were based on longitudinal ARIES cohort data [15] available at ages 7, 9, and 15–17. Epigenetic age (EA) was calculated using the Horvath clock [3] or GrimAge [9] in separate simulations. The original EA measure was then altered based on a simulated exposure. In each binary exposure simulation, a random n = 100 individuals had their original EA increased by 2 years (fixed effect), which accumulated by 0.1 year of EA per year of life (interaction effect). In each continuous exposure, all individuals were assigned a value (N(3.5, 0.5²)) which impacted their original EA by 0.1 years, times the level of exposure (fixed effect), and caused an interaction between the exposure and age by 0.02 (interaction effect). In the next step of our simulation, a series of methods was applied to model the simulated effects. Models are linear mixed effect models (LME), generalized estimating equations (GEE), and regression on difference between two epigenetic age (EA) measures (Δ aging). Outcome variables included are epigenetic age acceleration (EAA, residual from regressing EA on age), or EA itself. We ran n = 1000 simulations for each exposure type (binary vs. continuous), epigenetic clock (Horvath vs. GrimAge2), and different scopes of data (two measures: Age 7 and 15–17 vs. three measures: Age 7, 9, and 15–17)

Simulations were based on longitudinal DNAm data from ARIES [15], condensed into EA measures. To compare EA derived from conceptually different epigenetic clocks, we investigated EA calculations from the Horvath clock [3], GrimAge [9], and their principal component (PC) versions [50] in separate simulation cycles. In our binary exposure simulations, 100 participants were randomly selected as “exposed” (n = 918 “unexposed”). To create an association between the exposure and the average outcome, EA was increased by a fixed effect of 2 years in those 100 “exposed” individuals ($\beta_{2}$). Additionally, to create an association between the exposure and change in the average outcome over time, an interaction with chronological age by 0.1 years EA per year of age was included ($\beta_{4}$). In our continuous exposure simulations, all participants were randomly assigned a value ($N\left( {3.5, 0.5^{2} } \right)$), which affected EA with a fixed effect coefficient of 0.1 ($\beta_{2}$) and an interaction coefficient of 0.02 ($\beta_{4}$).

We then iterated (n = 1000) through all (i) models, (ii) outcomes and (iii) time variables discussed above, including (iv) different numbers of repeated measures (two timepoints vs. three timepoints). First, we evaluated three models (i): LME models (with and without random slope term), GEE, and regression on Δ aging. Second, we assessed two outcomes within these models (ii): EA an EAA. Third, we investigated four time variables (iii): age, years between measures, numerical ranks {1, 2, 3}, and factorized values {F07, F09, F15}. Δ aging was calculated with and without adjusting for the time between initial and follow-up measure. Fourth, we evaluated the impact of repeated measures (iv) on model accuracy by fitting LME models and GEE with two or three timepoints. Δ aging was limited to only two timepoints.

To compare robustness across methods and variables, we extracted fixed ($\beta_{2}$) and interaction ($\beta_{4}$) effect estimates from all models and evaluated how accurate they met the simulated exposure effect. We measured each model’s performance by comparing whether the 95% confidence interval (CI) contained the simulated effect size. Models resulting in CI that were fully above or below the simulated effect were labeled as “inflated” or “deflated,” while models resulting in CI that contained the simulated effect were defined as performed best.

Real-world example

To apply our simulation results in a real-world example, we used sex and birthweight as accessible biological parameters and examples of binary and continuous predictors of longitudinal EA. Longitudinal EAA or EA at ages 7 and 15–17, derived from the Horvath clock [3], GrimAge [9], and their PC versions [50], were modeled as the outcome. We applied all models and time variables discussed above and included sex and cell type proportions as covariates. Cell counts were estimated using the Houseman algorithm [51] applied to ALSPAC DNAm data with a peripheral blood reference [52].

In LME models and GEE, the fixed ($\beta_{2}$) and interactive (Formula 4: $\beta_{9}$, Formula 5: $\beta_{10}$) effects were estimated as:

Binary exposure (Biological Sex):

$$\begin{aligned} Outcome_{ij} & = \beta_{0i} + \beta_{1} Time Variable_{ij} + \beta_{2} Sex_{i} + \beta_{3} NK_{ij} \\ & \quad + \beta_{4} Granu_{ij} + \beta_{5} Mono_{ij} + \beta_{6} CD4T_{ij} + \beta_{7} CD8T_{ij} \\ & \quad + \beta_{8} Bcell_{ij} + \beta_{9} \left( {Time Variable_{ij} \times Sex_{i} } \right) + \epsilon_{ij} \\ \end{aligned}$$

(4)

Continuous exposure (Birthweight):

$$\begin{aligned} Outcome_{ij} & = \beta_{0i} + \beta_{1} Time Variable_{ij} + \beta_{2} Birthweight_{i} + \beta_{3} Sex_{i} \\ & \quad + \beta_{4} NK_{ij} + \beta_{5} Granu_{ij} + \beta_{6} Mono_{ij} + \beta_{7} CD4T_{ij} + \beta_{8} CD8T_{ij} \\ & \quad + \beta_{9} Bcell_{ij} + \beta_{10} \left( {Time Variable_{ij} \times Birthweight_{i} } \right) + \epsilon_{ij} \\ \end{aligned}$$

(5)

where $Outcome_{ij}$ is EA or EAA, measured in individual i, at timepoint j.

Different trends ($\beta_{1}$) in Δ aging (with and without adjusting for time between measures) across groups were modeled using linear regression:

Binary exposure (Biological Sex):

$$\Delta Aging_{{\text{i}}} = {\upbeta }_{0} + \beta_{1} Sex_{i} + \epsilon_{i}$$

(6)

Continuous exposure (Birthweight):

$$\Delta Aging_{{\text{i}}} = {\upbeta }_{0} + \beta_{1} Birthweight_{i} + \beta_{2} Sex_{i} + \epsilon_{i}$$

(7)

where $\Delta Aging_{{\text{i}}}$ is the difference score between two measures for individual i.

Results

Simulation study

Figure 3 and Table 2 summarize the effect estimates for a binary exposure for all models and time variables considered, using a subset of two within-person EA measures derived from the Horvath clock. Additional tables, summarizing the complete range of effect estimates, including the continuous exposure and models involving three within-person measures of EA derived from the Horvath clock as well as GrimAge, can be found in the supplemental material (Supplement Figs.1–24 and Tables 1–6). Across all approaches, the choice of time variable had the most substantial impact on the effect estimate bias, which stayed consistent across models (LME models, GEE and Δ aging) and outcome variables (EA, EAA). Including chronological age in the model gave most accurate estimates, while other age-independent time variables led to inflated results.

Fig. 3

Summary of n = 1,000 simulations from ARIES data (two measurements, age 7 and 15 or 17) [15]. Rows show the distribution (median, 25th and 75th percentile, outliers) of effect size estimates derived from different models and time variables included in those models, respectively. The two columns differentiate between estimates of the interaction term as well as the fixed effect. Simulated effect sizes are marked in red (interaction = 0.1; fixed effect = 2.0). Time variables are chronological age (Age), years between measures (Years), number of measure (Timepoint, i.e., 1, 2), factorized measure (Timefactor, i.e., F07, F15). Models are linear mixed effect models (LME), generalized estimating equations (GEE), and regression on difference between two epigenetic age (EA) measures (Δ aging). All models contained Horvath clock derived EA as outcome [3]

Table 2 Summary of n = 1000 simulations from ARIES data (two measurements, age 7 and 15 or 17) [15]

For both fixed and interaction effects alike, we obtained similar estimates across LME models and GEE, holding the time variable constant. Models that included timepoint as their time variable resulted in a slightly higher proportion of deflated fixed effect estimates than models including age, when using two within-person measures (LME model with random slope, binary simulation: timepoint 5% deflated, age 3% deflated; Supplement Fig. 9). The majority of estimates from models accounting for years between measures or categorical time overestimated the simulated fixed effect (LME model with random slope, binary simulation: > 60% inflated; Supplemental Fig. 9). All estimates from models incorporating time as a categorical variable were inflated.

Regression on Δ aging showed high precision and accuracy for estimating the interaction coefficients, when adjusted for years between measures. Similar performance in approximating the age-exposure interaction was achieved using LME models and GEE accounting for time through either chronological age or years between measures. Including either categorical time or numerical timepoint led to overestimated effects across all models, resulting in up to 52% inflated interaction estimates (LME model with random slope, binary exposure: timepoint, Supplemental Fig. 9). Models based on three repeated within-person measures resulted in exclusively deflated interaction estimates, when including categorical time as the time variable.

Real-world examples

Supplemental Tables 7–10 as well as Supplemental Figs. 25–28 show the complete range of effect size estimates for both exposures on longitudinal EA, across models, clocks, and time variables. Again, the choice of time variable had the most substantial impact on effect estimates, which stayed consistent across models.

Figure 4A and Supplemental Table 7 show the fixed effect estimates of biological sex on longitudinal EA between the ages 7 and 15–17 years across methods. Different model choices led to very similar estimates of fixed effects, so we compared the choice of time variable within the LME model (random slope) using Horvath EA as the outcome. Individuals assigned male at birth had 0.12 years lower average EA compared to individuals assigned female at birth (95% CI − 0.72, 0.47) using chronological age as the time variable, which based on our simulation study is an unbiased estimate. With the age-derived result as a reference, using either years between measures or categorical time resulted in estimates three times larger and in the opposite direction (males were on average 0.36 [95% CI 0.05, 0.66] years older). The use of timepoint led to estimates almost three times larger in the same direction (males were on average 0.36 [95% CI − 1.07, 0.36] years younger). While point estimates differed using GrimAge (likely due to the different training outcome measures of these clocks), the choice of time variable had a similar effect as seen when using the Horvath clock. In both Horvath and GrimAge EA measure-based models, including age or timepoint as a time variable led to much wider confidence intervals compared to using categorical time or time between measures.

Fig. 4

Age interaction and fixed effect estimates of the effect of biological sex (A) and birthweight (B) on EA over time. Included were two within-person measures (age 7 and 15 or 17) from the ARIES cohort [15]. Rows contain effect size point estimates as well as 95% confidence intervals derived from different models and time variables included in those models, respectively. Effect estimates are based on individuals assigned female at birth as reference. Significant estimates are marked in red (p < 0.05). Models are linear mixed effect models (LME), generalized estimating equations (GEE), and regression on difference between two epigenetic age (EA) measures (Δ aging). Time variables are chronological age (Age), years between measures (Years), number of measure (Timepoint, i.e., 1, 2), and factorized measure (Timefactor, i.e., F07, F15). Models contained EA measures derived from the Horvath clock [3]. All models were corrected for cell-type proportion, while models estimating the effect of birthweight additionally account for biological sex

Interaction effects between sex and age, leading to an accumulating positive or negative effect on EA over time, are shown in the left column of Fig. 4A and in Supplemental Table 8. The choice of model or outcome led to approximately equal estimates holding the time variable constant. Once again, while effect sizes differed profoundly between clocks, the effect of different time variables was similar using Horvath or GrimAge. Using LME models with age as the time variable, EA in individuals assigned male at birth increased by an extra 0.07 years per year of age, compared to those of individuals assigned female at birth (95% CI 0.02, 0.12). Using either years as the time variable, or Δ aging, led to the same point estimate. However, using either timepoint or categorical time inflated the interaction effect drastically (0.72, 95% CI 0.23, 1.21).

The fixed effect estimates of birthweight on longitudinal EA between the age 7 and 15–17 years of age are shown in the right column of Fig. 4B and in Supplemental Table 9. Due to similar estimates across models and outcome choices, we again compared results for different time variables within LME models (random slope) using Horvath EA as the outcome. Using chronological age as the time variable showed that an increase in birthweight of 1 kg is associated on average with an additional 1.08 years of EA (95% CI: 0.48, 1.69). With the age-derived estimate as a reference (assuming it is unbiased as per our simulation results), using either years between measures or categorical time as the time variable led to an estimate 2.5 times smaller (0.43 [95% CI 0.11, 0.74] years average increase per 1 kg). The use of timepoint increased estimates by 25% (1.3 [95% CI 0.56, 2.03] years average increase per kg). Models using GrimAge measures led to opposite effect estimates, none of which showed statistical significance. Trends across time variables were similar between Horvath and GrimAge EA measure-based models, leading to much wider confidence intervals when including age or timepoint as time variable compared to using categorical time or time between measures.

Estimates of interaction effects between birthweight and age are shown in the left column of Fig. 4B and in Supplementary Table 10. Again, the selection of either the model or the outcome led to similar estimates when keeping the time variable constant. Using LME models with age as the time variable, EA decreased on average by 0.09 years per year of age for each kg increase in birthweight (95% CI − 0.14, − 0.04). Using years between measures, or age-adjusted Δ aging, as the time variable led to a similar point estimate. Regressing on non-age-adjusted Δ aging or including timepoint or categorical time in the model inflated the interaction effect by a factor of 10 (− 0.87, 95% CI − 1.37, − 0.37).

Discussion

The findings of this study present valuable insight into the intricacies of modeling associations with longitudinal epigenetic age, emphasizing the critical influence of different exposure and time variables on the robustness of effect estimates and conclusions.

Simulation and real-world examples

All simulation and real-word analyses showed consistent estimates across models and outcome variables, while the choice of time variable significantly impacted the accuracy and precision of effect estimates. In simulations, including chronological age as the time variable in an LME model or GEE led to the highest number of correct fixed effect and interaction effect estimates. Hence, models accounting for chronological age also produced the most robust results in our real-world analysis of biological sex and birthweight. When estimating interaction effects, i.e., the accumulating effect that an exposure has on the pace of longitudinal EA, models including timepoint or categorical time had the highest number of incorrect effect size estimates (inflated or deflated) in simulations. We observed a difference of similar magnitude in interaction effect estimates in both real-world analyses compared to estimates from models using chronological age. Including years between measures or categorical time led to > 60% inflated fixed effect estimates and narrow confidence intervals in our binary exposure simulation. The resulting bias led to what appears like a false positive finding in the real-world analysis of sex.

Application examples: epidemiological conclusion

The two application examples not only lend support to our simulation findings, but also yield several novel epidemiological findings. First, our study suggests that on average, individuals assigned male at birth have higher GrimAge than individuals assigned female at birth (2.83 additional years of EA, 95% CI 2.41, 3.26). However, individuals assigned male at birth showed a decelerated rate of GrimAge over the time observed. More specifically, individuals assigned male at birth showed 0.08 years EA decrease per year of life on average (95% CI 0.05–0.11 decrease per year), compared to individuals assigned female at birth. Both effects remained consistent and comparable in magnitude when using EA estimates from the PC GrimAge clock. Biological aging studies focusing on other molecular biomarkers to investigate sex-specific differences support the lower baseline aging in women [53, 54]. Although cross-sectional studies across different age groups, ancestries, and epigenetic clocks have identified an association between biological sex and EA [34, 55, 56], no prior study, to our knowledge, has shown that the effect changes longitudinally. While an age interaction effect was detected with Horvath’s clock, it did not remain consistent when using its PC version. The differences in results between these clocks are likely due to their conceptual distinctions: the Horvath clock is designed to best predict chronological age, whereas the GrimAge clock is trained to predict health outcomes and lifespan. The variations in effect direction and magnitude observed across different clocks in our analysis highlight the importance of carefully selecting a clock that is appropriate for the specific study design and population before conducting the analysis [57].

Second, our results suggest a positive association between birthweight and longitudinal EA based on Horvath EA measures (1.08 additional years of EA per increase in kg birthweight, 95% CI 0.48, 1.69), as well as a negative birthweight/age interaction effect over time (− 0.09 years decrease in EA per year per increase in kg birthweight, 95% CI − 0.14, − 0.04). These results indicate that children born with higher birthweight have higher EA on average, while their pace of EA appears to slow down over the period measured, compared to children born with lower birthweight. However, the pattern does not remain consistent when using EA estimates from the PC Horvath clock. Results from models using EA estimates from PC GrimAge indicate lower baseline EA for individuals with higher birthweight, which aligns with other studies that have examined the relationship between birthweight and EA, and generally report that lower birthweight is associated with higher EA [58,59,60,61]. However, Simpkin et al. [34] found that birthweight was positively associated with EA at age 7, but negatively associated at age 17. Most studies analyzed data cross-sectionally, included dichotomized birthweight rather than continuous, and were based on EA measures in adult cohorts. Furthermore, associations were predominantly identified in individuals assigned male at birth [59, 60] or in male-dominated cohorts [58]. Future large-scale studies are needed to clarify longitudinal relationships and explore the effect birthweight has on EA throughout the life course.

Recommendations for future studies

As we recognize the importance of methodological choices, this section offers recommendations and guidance for researchers embarking on similar investigative paths. We highly recommend the use of LME models or GEE, including chronological age as the time variable, for studies working with repeated EA measures. These approaches improve precision and accuracy of fixed and interaction effect estimates. Alternatively, research evaluating interaction effects based on only two within-person measures can yield similar validity by using linear regression on age-adjusted Δ aging. We acknowledge that due to limited data collection or access, it is not always possible to implement the best-possible model. In cases where chronological age is not accessible and the cohort under study was measured in synchronized waves, we recommend the use of numerical timepoint to get more accurate fixed effect estimates. The effect size might appear slightly attenuated compared to models including age but is less susceptible to false positives. In studies aimed at exploring interaction effects, it is advisable to opt for years between measures instead of timepoint in the absence of chronological age. Factorized categorical time should be avoided due to its potential to introduce bias in fixed effect and interaction effect estimates, especially in studies incorporating more than two within-person measurements. When EA itself is not the preferred outcome, we strongly recommend using EAA, the residual from regressing epigenetic age on chronological age, instead of difference scores (EA minus chronological age). Residualizing outcome variables is a well-established practice in epigenetics epidemiology [62,63,64] to capture part of the EA estimate that is not associated with age, making it more informative about factors beyond age that may impact biological aging. Although the underlying interpretation of results remains consistent, the presentation differs. An effect on longitudinal EA indicates how a factor influences the overall biological age estimate over time, while an effect on longitudinal EAA reflects how this factor impacts the deviation between age and EA over time.

Strengths and limitations

The main strength of our study lies in its comprehensive evaluation of longitudinal EA models, thereby offering valuable insight to direct future research toward more reliable results. The increasing accessibility of repeated DNAm measures and growing interest in comprehending the effects of exposures on EA and its subsequent influence on health have shown the necessity for guidelines to address the issue robustly. A further strength of our work is the incorporation of both simulation and application in two real-world examples, which support the credibility and applicability of our results. One limitation worth noting is that the real-world examples were modeled using EA measured at only two time points, which might have impacted the comprehensiveness of our epidemiological findings. However, our simulation based on three repeated EA measures suggests results from models including age as time variable are consistent across different scopes of data. It also reflects the reality faced by many researchers, as most studies typically have access to only two measures. Second, we limited our methodological evaluation to common models and time variables found in recent literature [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44] and excluded rare and simplistic approaches. While there are certainly various other methods available, we assume that our study has addressed the applications most pertinent to most epidemiological studies. We furthermore limited our study to four epigenetic clocks, instead of evaluating all available options. The primary goal of this work is to explore the discrepancies that arise from using different methods to analyze associations with epigenetic age, independent of the specific clock employed. Our selection includes both first- and second-generation clocks, ranging from those designed to predict chronological age to those predicting mortality risk, as well as their principal component versions, which are increasingly favored in longitudinal studies. The choice of an epigenetic clock is generally context-dependent, influenced by study design and available data [57]. Our work aims to provide methodological guidance for researchers after they have selected the appropriate clock for their study.

Conclusions

In conclusion, our study presents a comprehensive evaluation of various methods utilized in modeling exposure effects on EA over time. Through a combination of simulation and real-world analyses, we have demonstrated that the methodological decisions made in longitudinal EA modeling significantly impact the reliability of effect estimates. Findings highlight that LME models or GEE, using chronological age as the time variable, are the optimal approach. Moreover, recognizing the constraints faced by some studies regarding data availability, we have provided practical recommendations to accommodate such limitations. Our thorough assessment serves as a valuable resource for guiding future epidemiological epigenetic aging research endeavors. By optimizing methodological approaches based on the insights from our study, researchers can enhance the depth and accuracy of their investigations, ultimately advancing our understanding of the complex interplay between exposures and epigenetic aging processes.