Train an ElasticNet DNA-methylation age clock on public, age-labeled blood data; benchmark it against the published clocks biolearn already ships; validate it on a truly external cohort it never saw; and document where it breaks and why.
The trained model is not the point. The failure analysis is. A clean leaderboard number on a held-out slice of your own training distribution is the part that lies to you; the external number and the residual structure are the honest parts.
TL;DR. A cohort-specific ElasticNet clock beats the general-purpose published clocks on its own cohort (held-out MAE 2.8–2.9 yr vs Horvath 4.5, Hannum 5.3). That's expected and not the story. The story is that it generalizes to a new cohort as a calibration failure: nearly unbiased internally, it over-predicts ages by ~3 years externally and over-predicts the young by ~6 years, while its rank-correlation barely moves. Where a clock "breaks" is mostly bias, not noise.
(Numbers below are written from the committed run — see results/headline.json
and results/*.csv. Reproduce with python scripts/run_pipeline.py.)
| Model | MAE (yr) | median AE | RMSE | Pearson r | bias |
|---|---|---|---|---|---|
| MyClock (clean) | 2.87 | 2.28 | 3.68 | 0.970 | +0.04 |
| MyClock (leaky twin) | 2.90 | 2.33 | 3.72 | 0.969 | +0.09 |
| Horvath v1 | 4.51 | 3.65 | 5.73 | 0.939 | −2.29 |
| Hannum | 5.31 | 4.91 | 6.47 | 0.958 | +4.68 |
| PhenoAge | 9.60 | 8.62 | 11.22 | 0.898 | −9.07 |
| Model | MAE (yr) | median AE | RMSE | Pearson r | bias |
|---|---|---|---|---|---|
| MyClock (clean) | 4.00 | 3.43 | 4.88 | 0.952 | +3.21 |
| Horvath v1 | 4.16 | 3.45 | 5.24 | 0.897 | +0.53 |
| Hannum | 7.39 | 7.19 | 8.47 | 0.920 | +7.10 |
| PhenoAge | 5.32 | 4.36 | 6.65 | 0.877 | −3.38 |
Read these two tables together. MyClock has the lowest MAE and highest correlation on both cohorts, but its bias jumps from +0.04 yr (internal) to +3.21 yr (external), while Horvath — barely worse on MAE — stays near-calibrated (+0.53). A clock can win on accuracy and still be the one that drifts most under distribution shift: internally MyClock beats Horvath by 1.6 yr MAE; externally that margin collapses to 0.16 yr.
Data (Phase A). All loaded through biolearn, which caches the underlying GEO downloads.
| Role | Cohort | Platform | n | Age | Notes |
|---|---|---|---|---|---|
| Train + held-out test | GSE40279 (Hannum 2013) |
450K | 656 | 19–101 (mean 64) | whole blood |
| External validation | GSE42861 |
450K | 689 | 18–70 (mean 52) | whole blood; RA + controls |
| Diversity stress-test | BoAChallengeData |
EPIC | (500) | — | loader returns empty matrix — see below |
The external cohort is deliberately different: younger, narrower age range, and a disease population. That mismatch is what makes its number honest.
Model (Phase C). ElasticNet (L1+L2), the standard clock architecture, fit the right way:
- Split first — an 80/20 train/test split (stratified on age quantiles) before any fitting decision. The test samples never inform feature selection, imputation, or scaling.
- Feature-select inside the training fold only — rank CpGs by |Pearson r| with age computed on train only, keep the top 10,000. (450K arrays have ~470k sites; you don't fit on all of them.)
- Tune by cross-validation —
ElasticNetCVselectsalphaandl1_ratioby internal 5-fold CV on the training fold. The L1 penalty then sparsifies the 10,000 pre-filtered CpGs down to the ~2,000 that actually carry the clock. - Evaluate on the untouched test fold, and on the external cohort.
Code: src/agingclock/model.py. The pipeline
(SimpleImputer → StandardScaler → ElasticNetCV) is fit on the training fold
only, so imputation medians and scaling come from train alone.
Baselines (Phase B). biolearn runs the published clocks (Horvath, Hannum, PhenoAge) with their original coefficients — we don't re-implement them. PhenoAge predicts phenotypic age, so a larger MAE against chronological age is expected, not a bug. Caveat we flag rather than hide: Hannum was originally trained on GSE40279, so its number on that cohort is an in-sample fit, not a fair held-out comparison.
biolearn's own training_simple_model.py template does bootstrap feature
selection on the entire dataset, then train_test_splits afterward. That is
data leakage: the held-out samples helped choose the features, so the reported
error is optimistic. We do it correctly (select inside the training fold) and
also keep a leaky twin that selects on the full data before splitting, to
measure the difference.
Finding — the leak is negligible here (|gap| ≈ 0.04 yr, sign not even
stable). The leaky twin's held-out MAE (2.90) is within four-hundredths of a
year of the corrected clock's (2.87) — and the difference is so small its sign
flips between configurations (the leaky twin came out 0.03 yr better at
max_iter=5000, 0.04 yr worse at the converged max_iter=20000). In other
words, the leakage signal is smaller than run-to-run noise. Why? The leak only
touches the pre-filter, and a correlation ranking over 524 training samples is
already stable — adding the 132 test samples barely reorders the top 10,000
CpGs, and ElasticNetCV does the real selection honestly either way.
The lesson is not "leakage doesn't matter" — it's that the magnitude depends on the mechanism. Leak the final model selection, or fit the scaler/imputer on a small sample including test, and the optimism is large; leak a stable correlation pre-filter on 500+ samples and it vanishes into noise. The methodological fix is still correct and worth making — but the honest report is that here it buys you ~0 years, not a dramatic correction. Quantifying that, rather than asserting a scary leakage number, is the point.
Each point is a claim with evidence. Figures are in figures/;
per-group numbers in results/.
1. Generalization is a calibration shift, not a collapse of signal. Moving from the held-out GSE40279 slice to the external GSE42861 cohort, MAE rises from 2.87 → 4.00 yr (+1.13). But Pearson r barely moves (0.970 → 0.952): the clock still ranks people by age almost as well. What breaks is the zero-point — bias goes from +0.04 yr internally to +3.21 yr externally. The external scatter (fig 3) sits as a near-parallel band above the y=x line. The clock didn't get noisier; it got shifted. ~3 of the 4 external MAE-years are a constant offset that a single per-cohort recalibration could remove.
2. The over-prediction is distribution shift, not disease. GSE42861 is half
rheumatoid-arthritis patients — an obvious hypothesis is that RA accelerates
epigenetic age and inflates the prediction. The data says no: RA subjects are
over-predicted less than healthy controls (mean residual +3.00 vs +3.44 yr;
results/bias_by_disease_external.csv). If accelerated aging were driving the
bias, the sick group would be more over-predicted, not less. The offset
tracks the cohort, not the pathology — pointing at the train/test age-distribution
mismatch (training mean 64 yr vs external mean 52 yr) as the cause.
3. It fails hardest at the young end — regression to the mean. Residuals vs
age (figs 2, 4; results/residuals_by_age_*.csv):
| Age band | Internal mean resid | External mean resid |
|---|---|---|
| <40 | +1.96 | +5.99 |
| 40–55 | +1.24 | +3.94 |
| 55–70 | +0.59 | +1.48 |
| >70 | −1.69 | (none in external) |
Internally the pattern is textbook regression to the mean: over-predict the young, under-predict the old, pivoting near the training mean (~64). Externally the young are over-predicted by ~6 years — the same downward-pull, amplified because the external cohort is full of under-40s the training set (mean age 64) barely contained. A clock learns the center of its training age distribution and drags the extremes toward it.
4. A sex differential in the external bias (~1.3 yr). On GSE42861 the
majority class (sex=0, n=492) is over-predicted by +3.58 yr vs +2.30 yr for
sex=1 (n=197); see fig 5 / results/bias_by_sex_external.csv. Internally the
sexes are balanced (−0.22 vs +0.29). The two ethnic groups present in GSE40279
(Caucasian-European, Hispanic-Mexican) show no meaningful internal bias gap
(−0.15 vs +0.45). So the one subgroup disparity that survives is sex-on-the-
external-cohort — though it is confounded with that cohort's sex × age structure
and needs a population with crossed sex/age to confirm. The bigger honest gap
is what's missing: the EPIC diversity cohort that would have tested
cross-platform and cross-ancestry generalization did not load (below), so the
populations most likely to expose bias went untested.
The cleanest, most actionable failure is #1+#2: cross-cohort generalization loss in this clock is dominated by a calibration shift (a near-constant age offset that tracks the cohort's age distribution), not by a loss of rank-correlation or by disease state.
Can cross-cohort calibration drift be removed without labels? If a methylation clock's external error is mostly a distribution-dependent bias (here ~3 of 4 MAE-years), can it be corrected using only unlabeled target methylation — e.g. aligning predictions to the target cohort's age distribution, or a domain-adaptation recalibration — rather than retraining on new labeled data? And does that correction generalize, or does each new cohort need its own?
This is a real research-version project: it predicts a concrete intervention (label-free recalibration), has a clear success metric (does the +3.2 yr external bias shrink while r is preserved?), and bears directly on whether published clocks can be trusted across populations. Two secondary questions fall out of #3 and #4 — can a clock be built to hold at the age extremes where training data is sparse? and are clocks systematically mis-calibrated across sex/ancestry, and is that correctable? — each worth its own study, and best attacked once the EPIC cross-platform data (the part that didn't load here) is in hand.
uv venv --python 3.11 .venv && source .venv/bin/activate
uv pip install -r requirements.txt
python scripts/run_pipeline.py # writes results/ and figures/ (~15 min)The first run downloads GEO data through biolearn (a few hundred MB, cached
afterward). Outputs: results/metrics_*.csv, results/headline.json, and
figures/*.png. A narrative walk-through is in
notebooks/aging_clock.ipynb.
src/agingclock/
data.py # cohort loaders (GSE40279, GSE42861, BoAChallengeData)
baselines.py # published clocks via biolearn.run_predictions
model.py # leakage-free ElasticNet clock + leaky twin
evaluate.py # metrics, scatter/residual/subgroup plots
scripts/run_pipeline.py # runs every phase end-to-end
notebooks/aging_clock.ipynb
results/ figures/
The reported numbers use biolearn 0.9.1 with pandas 2.2.3 / numpy 2.1.3
(pinned in requirements.txt). pandas 3.0 breaks biolearn's on-disk object cache
(a StringDtype deserialization error), so newer pandas is not yet usable with
this stack.
- One external cohort. GSE42861 is a single study; cross-cohort numbers would be firmer averaged over several. It is also a 450K-array, whole-blood, largely-European cohort — the same regime as training — so the +1 yr gap is a lower bound on real-world generalization loss.
- The BoA EPIC stress-test did not run (upstream loader bug, diagnosed and reported as biolearn#208; see below), so the cross-platform (450K→EPIC) and cross-population checks are missing — exactly the conditions most likely to break the clock.
- Hannum's baseline is in-sample on GSE40279 (noted above).
The spec's diversity stress-test cohort does not load in biolearn 0.9.1:
DataLibrary().get("BoAChallengeData").load() returns a methylation matrix of
shape (930659, 0) — 930k CpGs and zero samples — with scrambled,
mostly-NaN metadata (1003 rows for 500 subjects). So predictions cannot be
generated for it and the EPIC diversity test could not be run as specified.
Tracing it down: the dataset's biomarkers-challenge-2024 parser hardcodes
GEO series-matrix line numbers that are off by one against the current
GSE246337_series_matrix.txt. It reads the betas-column key from line 53
(!Sample_scan_protocol) instead of line 54 (!Sample_description, the sentrix
array IDs), so the betas-column → GSM mapping comes back empty and all 500
sample columns are pruned. The same offset shifts every metadata field onto the
wrong row. Reproducing the mapping at the corrected lines (54 / 32) yields the
full 500 sentrix-ID → GSM pairs — a one-line-numbers fix.
Reported upstream with reproduction, evidence, and a verified fix: bio-learn/biolearn#208. In its place this project reports subgroup bias on GSE42861 (disease, sex) and GSE40279 (ethnicity, sex) — see Failure modes.
MIT — see LICENSE.