diff --git a/.Rbuildignore b/.Rbuildignore
index fb47cd7..7e92b8b 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,3 +1,8 @@
 ^ClockSHAP\.Rproj$
 ^\.Rproj\.user$
 ^LICENSE\.md$
+^README\.Rmd$
+^\.github$
+^data-raw$
+^input$
+^dummy-waterfall.*\.png$
diff --git a/DESCRIPTION b/DESCRIPTION
index 75cd694..bb4d843 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -30,4 +30,5 @@ Imports:
     utils
 Depends:
     R (>= 3.6.0)
+LazyData: true
 Config/testthat/edition: 3
diff --git a/R/data.R b/R/data.R
new file mode 100644
index 0000000..4793c5e
--- /dev/null
+++ b/R/data.R
@@ -0,0 +1,39 @@
+#' Synthetic ClockSHAP example dataset
+#'
+#' A small, fully synthetic dataset used in the package examples and vignette.
+#' It contains a cohort of 120 samples with eight genes carrying genuine linear
+#' age trends, plus a ridge-fit linear clock that predicts age from those genes.
+#' Two samples are planted deviators: \code{S001} reads older than expected for
+#' its age and \code{S002} reads younger, both driven by a shared
+#' senescence / inflammaging axis (\code{CDKN2A}, \code{GDF15}, \code{IL6}).
+#'
+#' @format A list with three elements:
+#' \describe{
+#'   \item{\code{features}}{A 120 x 8 numeric matrix of synthetic gene
+#'     expression. Rows are samples (\code{S001}-\code{S120}); columns are genes
+#'     (\code{CDKN2A}, \code{GDF15}, \code{EDA2R}, \code{IL6}, \code{LMNB1},
+#'     \code{ELN}, \code{COL1A1}, \code{FOXO3}).}
+#'   \item{\code{age}}{Numeric vector of length 120, chronological age in years.}
+#'   \item{\code{clock}}{A \code{linear_clock} object (intercept, per-gene
+#'     weights, centering and scaling) fit by ridge regression so predicted age
+#'     tracks chronological age with realistic regression to the mean
+#'     (slope < 1).}
+#' }
+#'
+#' @details The data are entirely synthetic and carry no real biological signal;
+#'   gene symbols are used only to make the example legible. Regenerated by
+#'   \code{data-raw/make_clockshap_example.R}.
+#'
+#' @examples
+#' ref <- fit_reference_profile(
+#'   features = as.data.frame(clockshap_example$features),
+#'   age      = clockshap_example$age
+#' )
+#' cs <- clockshap(
+#'   features  = clockshap_example$features,
+#'   age       = clockshap_example$age,
+#'   clock     = clockshap_example$clock,
+#'   reference = ref
+#' )
+#' cs$deviation
+"clockshap_example"
diff --git a/README.Rmd b/README.Rmd
new file mode 100644
index 0000000..c60cc4b
--- /dev/null
+++ b/README.Rmd
@@ -0,0 +1,201 @@
+---
+output: github_document
+---
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(
+  collapse  = TRUE,
+  comment   = "#>",
+  fig.path  = "man/figures/README-",
+  out.width = "100%",
+  dpi       = 150
+)
+library(ClockSHAP)
+```
+
+# ClockSHAP
+
+ClockSHAP is an interpretability framework for **linear-additive aging clocks** (e.g., transcriptomic, epigenetic, or proteomic clocks fit with linear, ridge, or elastic-net regression).
+
+It answers a single question:
+
+> **Why is this sample predicted to be biologically older or younger than expected for its chronological age?**
+
+The method was developed for tissue-anchored transcriptomic aging clocks but is designed to support comparative aging analyses across diverse clock contexts.
+
+![ClockSHAP method overview](man/figures/clockshap-overview.png)
+
+*ClockSHAP applied to a lung-tissue aging clock on real GTEx-LUNG data. (A) Predicted vs. chronological age across the cohort, with one positive and one negative deviator highlighted. (B) A single feature's contribution, computed as its clock weight times its distance from the age-expected value. (C, D) The full per-feature decompositions for the two deviators, with contributions summing exactly to each sample's age-matched deviation. Figure from [Bock et al., in preparation].*
+
+---
+
+## What you get out of ClockSHAP
+
+For each sample, ClockSHAP returns its **age-matched deviation** (how far the predicted clock age sits above or below what is expected at that chronological age) together with **per-feature contributions that sum exactly to that deviation**.
+
+The figure above shows this on real lung-tissue data for a pathway-level transcriptomic clock. The positive deviator (Panel C) reads about 4.7 years older than expected, attributed to programs such as epithelial–mesenchymal transition (+1.9 years) and fatty-acid metabolism (+0.9); the negative deviator (Panel D) reads about 4.9 years younger, driven largely by low EMT (−4.0). In both cases the contributions sum *exactly* to the deviation, a guarantee of the method rather than an approximation.
+
+This output feeds downstream analyses: comparing which features drive deviation across a cohort, deep dives on individual samples, or hypothesis generation about which biological programs are differentially engaged.
+
+Because the reference here is the GTEx-LUNG cohort itself, the deviations shown are *exactly* the commonly used **relative age acceleration** (RAA), the residual of predicted age regressed on chronological age, which is why Panel A's mean deviation is zero. See [Terminology: deviation vs. AAA vs. RAA](#terminology-deviation-vs-aaa-vs-raa) below for the details.
+
+---
+
+## Installation
+
+```r
+# Recommended
+install.packages("pak")
+pak::pak("montilab/ClockSHAP")
+
+# Alternative
+install.packages("remotes")
+remotes::install_github("montilab/ClockSHAP")
+```
+
+---
+
+## Quick start
+
+ClockSHAP needs three things: a **feature matrix** (samples in rows, features in columns), a vector of **chronological ages** (one per sample), and a **clock** describing how those features map to predicted age. The bundled `clockshap_example` supplies all three, so we can look at the expected shapes first:
+
+```{r inputs}
+library(ClockSHAP)
+
+# Features: rows are samples, columns are the clock's features.
+# (If your matrix is genes x samples, transpose it first.)
+clockshap_example$features[1:3, 1:4]
+
+# One chronological age per sample
+head(data.frame(
+  sample = names(clockshap_example$age),
+  age    = unname(clockshap_example$age)
+))
+
+# The clock: an intercept (alpha), per-feature weights (beta), and the
+# training mean/SD used to standardize each feature (mu, sigma)
+str(clockshap_example$clock)
+```
+
+For your own model you build this object with `linear_clock(alpha, beta, mu, sigma)` from its fitted coefficients; the feature columns just have to be named to match. With those inputs in hand, the workflow is three calls: fit an age-conditioned reference from the data, run the decomposition, and inspect a sample.
+
+```{r quick-start, fig.alt = "Waterfall decomposition of one synthetic sample's age-matched deviation into per-gene contributions"}
+# Age-conditioned reference profile (defines "expected for age"), fit from the data
+ref <- fit_reference_profile(
+  features = as.data.frame(clockshap_example$features),
+  age      = clockshap_example$age
+)
+
+# Decompose each sample's age-matched deviation into per-feature contributions
+cs <- clockshap(
+  features  = clockshap_example$features,
+  age       = clockshap_example$age,
+  clock     = clockshap_example$clock,
+  reference = ref
+)
+
+# Inspect one sample's decomposition as a waterfall
+plot_clockshap_waterfall(cs, sample = 1, top_n = 6)
+```
+
+The returned `clockshap` object is a list with five fields:
+
+- `predicted`: clock prediction $\hat{y}$
+- `expected`: age-matched expectation $\mathrm{E}[\hat{y} \mid \mathrm{Age}]$ under the reference
+- `deviation`: `predicted - expected` ($\Delta$)
+- `phi`: per-feature contribution matrix (rows sum exactly to `deviation`)
+- `age`: input chronological ages
+
+---
+
+## How it works (in two lines)
+
+ClockSHAP is built around two identities. The first defines what it explains: the age-matched deviation of a sample's predicted age from its age-matched expectation,
+
+$$\Delta_i = \hat{y}_i - \mathrm{E}[\hat{y} \mid \mathrm{Age}_i]$$
+
+The second is the exact additive decomposition of that deviation into per-feature contributions,
+
+$$\sum_k \phi_{i,k} = \Delta_i$$
+
+Here $\hat{y}$ is the clock's predicted age, $\mathrm{E}[\hat{y} \mid \mathrm{Age}]$ is the age-matched expectation under a reference, and $\phi$ is a single feature's contribution. Everything below is about how the expectation and the contributions are defined, and why the decomposition is exact.
+
+---
+
+## The building blocks
+
+The rest of this README explains the two identities above:
+
+- what `expected` means (the reference profile),
+- why the per-feature contributions $\phi_k$ sum exactly to the deviation (the SHAP property),
+- how the deviation differs from other clock-acceleration metrics (AAA, RAA).
+
+### How `expected` is defined: the reference profile
+
+Aging-clock features change systematically with chronological age, so a single global baseline is inappropriate for comparing samples of different ages. ClockSHAP instead uses an **age-conditioned reference**. For each feature $k$ it fits an age trend on a reference cohort:
+
+$$\mathrm{E}[X_k \mid \mathrm{Age}] = \gamma_{0,k} + \gamma_{1,k}\,\mathrm{Age}$$
+
+Passing this age-expected feature vector through the clock $f$ gives the age-matched expectation, and the deviation is the gap from it:
+
+$$\mathrm{E}[\hat{y} \mid \mathrm{Age}] = f\big(\mathrm{E}[X \mid \mathrm{Age}]\big), \qquad \Delta = \hat{y} - \mathrm{E}[\hat{y} \mid \mathrm{Age}]$$
+
+So `expected` is the clock's *typical output for reference individuals of the same chronological age*, and the deviation is how far above or below that a sample lies.
+
+**Reference note.** Because `expected` depends on the chosen reference (e.g., GTEx tissue-of-origin), the deviation should be interpreted **relative to that reference**; it is not a universal biological-age score.
+
+**Interpretation note.** The deviation is a cross-sectional **state shift** relative to an age-matched expectation; it does not by itself imply a within-individual aging *rate*.
+
+### Why the contributions sum exactly to the deviation: the SHAP property
+
+For **linear-additive models**, SHAP-style attributions have a closed form. ClockSHAP's per-feature contribution is the feature's clock weight $w_k$ times its distance from the age-expected value, the form shown in Panel B of the figure above:
+
+$$\phi_{i,k} = w_k\big(x_{i,k} - \mathrm{E}[X_k \mid \mathrm{Age}_i]\big)$$
+
+Summed over features, the contributions recover the deviation exactly:
+
+$$\sum_k \phi_{i,k} = \Delta_i$$
+
+**Implementation note.** In practice the features are standardized to the clock's training distribution (mean and standard deviation) before the decomposition is computed. Changing the reference changes the expectation, and therefore the decomposition target, so every component is reference-anchored.
+
+### Terminology: deviation vs. AAA vs. RAA
+
+ClockSHAP explains $\Delta = \text{predicted} - \text{expected}$, not raw prediction error. Two other quantities appear in the aging-clock literature:
+
+**Absolute age acceleration (AAA, sometimes called ΔAge):**
+
+$$\mathrm{AAA} = \hat{y} - \mathrm{Age}$$
+
+Simple, but age-biased when the clock slope is below 1: younger samples are over-predicted and older samples under-predicted.
+
+**Relative age acceleration (RAA, regression residuals):**
+
+```r
+RAA <- residuals(stats::lm(predicted ~ age))
+```
+
+A post-hoc calibration that removes systematic age trends in prediction error within a chosen cohort.
+
+**ClockSHAP deviation.** ClockSHAP's $\Delta$ generalizes RAA. The age-conditional expectation is built through a **feature-level reference profile** ($\gamma_0 + \gamma_1\,\mathrm{Age}$ per feature, passed through the clock) rather than by regressing $\hat{y}$ on age directly. This has an exact special case: **for a linear-additive clock, when the reference is fit on the cohort of interest itself, $\Delta$ equals RAA sample-for-sample.** Because $\hat{y}$ is an exact linear function of the features and OLS projection onto age is linear, passing the per-feature age trends through the clock reproduces the OLS fit of $\hat{y} \sim \mathrm{Age}$ exactly, so $\text{predicted} - \text{expected} = \text{predicted} - \text{fitted} = \text{residual} = \mathrm{RAA}$.
+
+An additional usage of ClockSHAP is what happens when the reference is set to a *different*, biologically meaningful cohort (e.g., GTEx tissue-of-origin): `expected` is then anchored to normal aging in that reference rather than to the target cohort's own internal trend, and the deviation measures departure from that external expectation.
+
+We recommend the name **"deviation"** or **"age-matched deviation"** to keep this distinct from AAA and RAA.
+
+---
+
+## Practical notes
+
+- **Deviation is reference-dependent.** The choice of reference cohort defines what "expected for age" means.
+- **"Years" are clock-relative.** A deviation of +5 under one clock is not comparable to +5 under a different clock, training set, or feature space.
+- **Upstream preprocessing matters a lot.** ClockSHAP assumes the feature matrix is analysis-ready and consistent with the clock (transformations, normalization, feature matching, batch correction, QC). Mismatched feature spaces or batch effects produce technical, not biological, deviations. The importance of this step cannot be overstated.
+
+For broader discussion of calibration, age bias, and computational challenges in clock analysis, see [Epigenetic ageing clocks: statistical methods and emerging computational challenges](https://doi.org/10.1038/s41576-024-00807-w) (Nature Reviews Genetics).
+
+---
+
+## Status
+
+Under active development.
diff --git a/README.md b/README.md
index ed6fd0b..b41f618 100644
--- a/README.md
+++ b/README.md
@@ -1,64 +1,66 @@
-# ClockSHAP
-
-ClockSHAP is an interpretability framework for **linear-additive aging clocks** (e.g., transcriptomic, epigenetic, or proteomic clocks fit with linear, ridge, or elastic-net regression).
-
-It answers a single question:
-
-> **Why is this sample predicted to be biologically older or younger than expected for its chronological age?**
-
-The method was developed in the context of normal tissue transcriptomic aging clocks applied to primary tumors but is designed to support comparative aging analyses across diverse clock contexts.
-
----
-
-## What you get out of ClockSHAP
-
-For each sample, ClockSHAP returns a per-feature decomposition of how far that sample lies above or below the clock's age-matched expectation. A worked example (illustrative numbers):
 
-```
-Sample ID:            TCGA-XX-1234
-Chronological age:    58 years
-Predicted clock age:  66.5 years   <- y_hat
-Expected for age 58:  61 years     <- y_exp(age) from reference profile
-Deviation:            +5.5 years   <- predicted - expected
-
-Top feature contributions (phi, sum exactly to deviation):
-  CDKN2A      +3.1
-  COL1A1      +2.4
-  IL6         +1.8
-  SOX2        -1.0
-  Other (remaining features) -0.8
-  ------------------------
-  Total       +5.5
-```
-
-The interpretation is direct: this sample is 5.5 clock-years older than expected for a 58-year-old in the reference cohort, and CDKN2A, COL1A1, and IL6 are the dominant drivers.
+<!-- README.md is generated from README.Rmd. Please edit that file -->
 
-This output drops into downstream analyses, e.g.: cohort-level comparisons of which features drive deviation, deep dives on individual samples, or hypothesis generation about which biological features are differentially engaged between samples or cohorts. 
+# ClockSHAP
 
-`plot_clockshap_waterfall(cs, sample = i)` renders this decomposition as a waterfall plot.
+ClockSHAP is an interpretability framework for **linear-additive aging
+clocks** (e.g., transcriptomic, epigenetic, or proteomic clocks fit with
+linear, ridge, or elastic-net regression).
 
-Worked-example visualization (same illustrative numbers as above):
+It answers a single question:
 
-![ClockSHAP worked example waterfall](man/figures/readme-worked-example-waterfall.png)
+> **Why is this sample predicted to be biologically older or younger
+> than expected for its chronological age?**
 
----
+The method was developed for tissue-anchored transcriptomic aging clocks
+but is designed to support comparative aging analyses across diverse
+clock contexts.
 
-## How it works (in two lines)
+![ClockSHAP method overview](man/figures/clockshap-overview.png)
 
-ClockSHAP is built around two equations:
+*ClockSHAP applied to a lung-tissue aging clock on real GTEx-LUNG data.
+(A) Predicted vs. chronological age across the cohort, with one positive
+and one negative deviator highlighted. (B) A single feature’s
+contribution, computed as its clock weight times its distance from the
+age-expected value. (C, D) The full per-feature decompositions for the
+two deviators, with contributions summing exactly to each sample’s
+age-matched deviation. Figure from \[Bock et al., in preparation\].*
 
-```text
-deviation     = predicted - expected
-sum_k phi[k]  = deviation
-```
+------------------------------------------------------------------------
 
-Everything else in this README is about how `expected` and `phi` are defined and why the additive decomposition holds exactly.
+## What you get out of ClockSHAP
 
----
+For each sample, ClockSHAP returns its **age-matched deviation** (how
+far the predicted clock age sits above or below what is expected at that
+chronological age) together with **per-feature contributions that sum
+exactly to that deviation**.
+
+The figure above shows this on real lung-tissue data for a pathway-level
+transcriptomic clock. The positive deviator (Panel C) reads about 4.7
+years older than expected, attributed to programs such as
+epithelial–mesenchymal transition (+1.9 years) and fatty-acid metabolism
+(+0.9); the negative deviator (Panel D) reads about 4.9 years younger,
+driven largely by low EMT (−4.0). In both cases the contributions sum
+*exactly* to the deviation, a guarantee of the method rather than an
+approximation.
+
+This output feeds downstream analyses: comparing which features drive
+deviation across a cohort, deep dives on individual samples, or
+hypothesis generation about which biological programs are differentially
+engaged.
+
+Because the reference here is the GTEx-LUNG cohort itself, the
+deviations shown are *exactly* the commonly used **relative age
+acceleration** (RAA), the residual of predicted age regressed on
+chronological age, which is why Panel A’s mean deviation is zero. See
+[Terminology: deviation vs. AAA
+vs. RAA](#terminology-deviation-vs-aaa-vs-raa) below for the details.
+
+------------------------------------------------------------------------
 
 ## Installation
 
-```r
+``` r
 # Recommended
 install.packages("pak")
 pak::pak("montilab/ClockSHAP")
@@ -68,151 +70,236 @@ install.packages("remotes")
 remotes::install_github("montilab/ClockSHAP")
 ```
 
----
+------------------------------------------------------------------------
 
 ## Quick start
 
-```r
+ClockSHAP needs three things: a **feature matrix** (samples in rows,
+features in columns), a vector of **chronological ages** (one per
+sample), and a **clock** describing how those features map to predicted
+age. The bundled `clockshap_example` supplies all three, so we can look
+at the expected shapes first:
+
+``` r
 library(ClockSHAP)
-set.seed(1)
 
-features <- matrix(rnorm(20), nrow = 5, ncol = 4)
-colnames(features) <- paste0("F", 1:4)
-age <- c(40, 50, 60, 70, 80)
+# Features: rows are samples, columns are the clock's features.
+# (If your matrix is genes x samples, transpose it first.)
+clockshap_example$features[1:3, 1:4]
+#>        CDKN2A    GDF15    EDA2R      IL6
+#> S001 7.034590 7.851806 5.223082 6.046552
+#> S002 5.641907 6.397235 5.266366 4.738019
+#> S003 5.293459 6.050086 3.950593 4.233756
+
+# One chronological age per sample
+head(data.frame(
+  sample = names(clockshap_example$age),
+  age    = unname(clockshap_example$age)
+))
+#>   sample  age
+#> 1   S001 80.3
+#> 2   S002 81.5
+#> 3   S003 45.7
+#> 4   S004 75.7
+#> 5   S005 65.3
+#> 6   S006 58.6
+
+# The clock: an intercept (alpha), per-feature weights (beta), and the
+# training mean/SD used to standardize each feature (mu, sigma)
+str(clockshap_example$clock)
+#> List of 4
+#>  $ alpha: num 59
+#>  $ beta : Named num [1:8] 3.43 3.08 2.65 1.67 -2.08 ...
+#>   ..- attr(*, "names")= chr [1:8] "CDKN2A" "GDF15" "EDA2R" "IL6" ...
+#>  $ mu   : Named num [1:8] 5.34 6.23 4.46 4.6 5.68 ...
+#>   ..- attr(*, "names")= chr [1:8] "CDKN2A" "GDF15" "EDA2R" "IL6" ...
+#>  $ sigma: Named num [1:8] 1.034 1.074 0.933 1.126 0.998 ...
+#>   ..- attr(*, "names")= chr [1:8] "CDKN2A" "GDF15" "EDA2R" "IL6" ...
+#>  - attr(*, "class")= chr "linear_clock"
+```
 
-clock <- linear_clock(
-  alpha = 10,
-  beta  = setNames(runif(4), colnames(features)),
-  mu    = setNames(rep(0, 4), colnames(features)),
-  sigma = setNames(rep(1, 4), colnames(features))
-)
+For your own model you build this object with
+`linear_clock(alpha, beta, mu, sigma)` from its fitted coefficients; the
+feature columns just have to be named to match. With those inputs in
+hand, the workflow is three calls: fit an age-conditioned reference from
+the data, run the decomposition, and inspect a sample.
 
+``` r
+# Age-conditioned reference profile (defines "expected for age"), fit from the data
 ref <- fit_reference_profile(
-  features = as.data.frame(features),
-  age = age
+  features = as.data.frame(clockshap_example$features),
+  age      = clockshap_example$age
 )
 
-cs <- clockshap(features = features, age = age, clock = clock, reference = ref)
-
-cs$phi
-cs$deviation
-cs$predicted
-cs$expected
+# Decompose each sample's age-matched deviation into per-feature contributions
+cs <- clockshap(
+  features  = clockshap_example$features,
+  age       = clockshap_example$age,
+  clock     = clockshap_example$clock,
+  reference = ref
+)
 
-plot_clockshap_waterfall(cs, sample = 1, top_n = 10)
+# Inspect one sample's decomposition as a waterfall
+plot_clockshap_waterfall(cs, sample = 1, top_n = 6)
 ```
 
+<img src="man/figures/README-quick-start-1.png" alt="Waterfall decomposition of one synthetic sample's age-matched deviation into per-gene contributions" width="100%" />
+
 The returned `clockshap` object is a list with five fields:
 
-- `predicted` — clock prediction `y_hat`
-- `expected` — age-matched expectation `y_exp(age)` under the reference
-- `deviation` — `predicted - expected`
-- `phi` — per-feature contribution matrix (rows sum exactly to `deviation`)
-- `age` — input chronological ages
+-   `predicted`: clock prediction $\hat{y}$
+-   `expected`: age-matched expectation
+    $\mathrm{E}[\hat{y} \mid \mathrm{Age}]$ under the reference
+-   `deviation`: `predicted - expected` ($\Delta$)
+-   `phi`: per-feature contribution matrix (rows sum exactly to
+    `deviation`)
+-   `age`: input chronological ages
 
----
+------------------------------------------------------------------------
 
-## The building blocks
+## How it works (in two lines)
 
-The rest of this README explains the two equations above:
+ClockSHAP is built around two identities. The first defines what it
+explains: the age-matched deviation of a sample’s predicted age from its
+age-matched expectation,
 
-- What `expected` means (the reference profile)
-- Why per-feature `phi` values exist and sum exactly to `deviation` (the SHAP property)
-- How `deviation` differs from other clock-acceleration metrics (AAA, RAA)
+$$\Delta_i = \hat{y}_i - \mathrm{E}[\hat{y} \mid \mathrm{Age}_i]$$
 
-### How `expected` is defined: the reference profile
+The second is the exact additive decomposition of that deviation into
+per-feature contributions,
 
-Aging-clock features change systematically with chronological age. A single global mean baseline is inappropriate for comparing samples of different ages, so ClockSHAP uses an **age-conditioned reference**.
+$$\sum_k \phi_{i,k} = \Delta_i$$
 
-For each feature `k`, we fit an age trend on a reference cohort:
+Here $\hat{y}$ is the clock’s predicted age,
+$\mathrm{E}[\hat{y} \mid \mathrm{Age}]$ is the age-matched expectation
+under a reference, and $\phi$ is a single feature’s contribution.
+Everything below is about how the expectation and the contributions are
+defined, and why the decomposition is exact.
 
-```text
-x_ref[k](age) = gamma0[k] + gamma1[k] * age
-```
+------------------------------------------------------------------------
 
-`gamma0[k]` and `gamma1[k]` are stored in the `reference_profile` object. Given a clock `f` mapping features to predicted age:
+## The building blocks
 
-```text
-y_exp(age)  = f( x_ref(age) )
-deviation   = y_hat - y_exp(age)
-```
+The rest of this README explains the two identities above:
 
-So `expected` is the clock's *typical output for reference individuals of the same chronological age*, and `deviation` is how far above or below that age-matched expectation a sample lies.
+-   what `expected` means (the reference profile),
+-   why the per-feature contributions $\phi_k$ sum exactly to the
+    deviation (the SHAP property),
+-   how the deviation differs from other clock-acceleration metrics
+    (AAA, RAA).
 
-**Reference note.** Because `expected` depends on the chosen reference (e.g., GTEx tissue-of-origin), `deviation` should be interpreted **relative to that reference**. It is not a universal biological-age score.
+### How `expected` is defined: the reference profile
 
-**Interpretation note.** `deviation` is a cross-sectional **state shift** relative to an age-matched expectation; it does not by itself imply a within-individual aging *rate*.
+Aging-clock features change systematically with chronological age, so a
+single global baseline is inappropriate for comparing samples of
+different ages. ClockSHAP instead uses an **age-conditioned reference**.
+For each feature $k$ it fits an age trend on a reference cohort:
 
-### Why per-feature contributions sum exactly to `deviation`: the SHAP property
+$$\mathrm{E}[X_k \mid \mathrm{Age}] = \gamma_{0,k} + \gamma_{1,k}\,\mathrm{Age}$$
 
-For **linear-additive models**, SHAP-style additive attributions have a closed form. Writing the clock on standardized features:
+Passing this age-expected feature vector through the clock $f$ gives the
+age-matched expectation, and the deviation is the gap from it:
 
-```text
-y_hat = alpha + sum_k beta[k] * z[k]
-z[k]  = (x[k] - mu[k]) / sigma[k]
-```
+$$\mathrm{E}[\hat{y} \mid \mathrm{Age}] = f\big(\mathrm{E}[X \mid \mathrm{Age}]\big), \qquad \Delta = \hat{y} - \mathrm{E}[\hat{y} \mid \mathrm{Age}]$$
 
-Given a standardized baseline `z_ref`, linear-model SHAP values reduce to:
+So `expected` is the clock’s *typical output for reference individuals
+of the same chronological age*, and the deviation is how far above or
+below that a sample lies.
 
-```text
-phi[k]        = beta[k] * ( z[k] - z_ref[k] )
-sum_k phi[k]  = y_hat(x) - y_hat(z_ref)
-```
+**Reference note.** Because `expected` depends on the chosen reference
+(e.g., GTEx tissue-of-origin), the deviation should be interpreted
+**relative to that reference**; it is not a universal biological-age
+score.
 
-ClockSHAP sets the baseline using the age-conditioned reference profile:
+**Interpretation note.** The deviation is a cross-sectional **state
+shift** relative to an age-matched expectation; it does not by itself
+imply a within-individual aging *rate*.
 
-```text
-x_ref[k](age)     = gamma0[k] + gamma1[k] * age
-z_ref[k](age)     = ( x_ref[k](age) - mu[k] ) / sigma[k]
+### Why the contributions sum exactly to the deviation: the SHAP property
 
-y_exp(age)        = alpha + sum_k beta[k] * z_ref[k](age)
-phi[k](age)       = beta[k] * ( z[k] - z_ref[k](age) )
-sum_k phi[k](age) = deviation
-```
+For **linear-additive models**, SHAP-style attributions have a closed
+form. ClockSHAP’s per-feature contribution is the feature’s clock weight
+$w_k$ times its distance from the age-expected value, the form shown in
+Panel B of the figure above:
+
+$$\phi_{i,k} = w_k\big(x_{i,k} - \mathrm{E}[X_k \mid \mathrm{Age}_i]\big)$$
 
-This preserves exact additivity while making the baseline biologically appropriate for aging-clock interpretation.
+Summed over features, the contributions recover the deviation exactly:
 
-**Implementation note.** Changing the reference changes `expected(age)` and therefore changes the decomposition target. Every component of the output is reference-anchored.
+$$\sum_k \phi_{i,k} = \Delta_i$$
 
-### Terminology: deviation vs. AAA vs. RAA
+**Implementation note.** In practice the features are standardized to
+the clock’s training distribution (mean and standard deviation) before
+the decomposition is computed. Changing the reference changes the
+expectation, and therefore the decomposition target, so every component
+is reference-anchored.
 
-ClockSHAP explains **`predicted − expected`**, not raw prediction error. Two other quantities are commonly used in the aging-clock literature, and it is worth being explicit about which one ClockSHAP targets.
+### Terminology: deviation vs. AAA vs. RAA
+
+ClockSHAP explains $\Delta = \text{predicted} - \text{expected}$, not
+raw prediction error. Two other quantities appear in the aging-clock
+literature:
 
 **Absolute age acceleration (AAA, sometimes called ΔAge):**
 
-```text
-AAA = y_hat - age
-```
+$$\mathrm{AAA} = \hat{y} - \mathrm{Age}$$
 
-Simple, but age-biased when the clock has slope < 1: younger samples get over-predicted and older samples under-predicted.
+Simple, but age-biased when the clock slope is below 1: younger samples
+are over-predicted and older samples under-predicted.
 
 **Relative age acceleration (RAA, regression residuals):**
 
-```r
+``` r
 RAA <- residuals(stats::lm(predicted ~ age))
 ```
 
-A post hoc calibration that removes systematic age trends in prediction error within a chosen cohort.
-
-**ClockSHAP deviation (`predicted − expected`):**
+A post-hoc calibration that removes systematic age trends in prediction
+error within a chosen cohort.
 
-ClockSHAP's `deviation` generalizes RAA. The age-conditional expectation is built through a **feature-level reference profile** (per-feature OLS trends `gamma0 + gamma1 * age`, passed through the clock) rather than by regressing `y_hat` on age directly. This has an exact special case: **for a linear-additive clock, when the reference profile is fit on the cohort of interest itself, `deviation` equals RAA sample-for-sample.** Because `y_hat` is an exact linear function of the features and OLS projection onto age is linear, passing the per-feature age trends through the clock reproduces the OLS fit of `y_hat ~ age` exactly, so `predicted − expected = predicted − fitted = residual = RAA`.
+**ClockSHAP deviation.** ClockSHAP’s $\Delta$ generalizes RAA. The
+age-conditional expectation is built through a **feature-level reference
+profile** ($\gamma_0 + \gamma_1\,\mathrm{Age}$ per feature, passed
+through the clock) rather than by regressing $\hat{y}$ on age directly.
+This has an exact special case: **for a linear-additive clock, when the
+reference is fit on the cohort of interest itself, $\Delta$ equals RAA
+sample-for-sample.** Because $\hat{y}$ is an exact linear function of
+the features and OLS projection onto age is linear, passing the
+per-feature age trends through the clock reproduces the OLS fit of
+$\hat{y} \sim \mathrm{Age}$ exactly, so
+$\text{predicted} - \text{expected} = \text{predicted} - \text{fitted} = \text{residual} = \mathrm{RAA}$.
 
-The purpose of ClockSHAP is what happens when the reference is set to a *different*, biologically meaningful cohort (e.g., GTEx tissue-of-origin). `expected` is then anchored to normal aging in that reference rather than to the target cohort's own internal trend, and `deviation` measures departure from that external expectation. 
+An additional usage of ClockSHAP is what happens when the reference is
+set to a *different*, biologically meaningful cohort (e.g., GTEx
+tissue-of-origin): `expected` is then anchored to normal aging in that
+reference rather than to the target cohort’s own internal trend, and the
+deviation measures departure from that external expectation.
 
-We recommend the name **"deviation"** or **"age-matched deviation"** to keep this distinct from AAA and RAA.
+We recommend the name **“deviation”** or **“age-matched deviation”** to
+keep this distinct from AAA and RAA.
 
----
+------------------------------------------------------------------------
 
 ## Practical notes
 
-- **Deviation is reference-dependent.** Choice of reference cohort defines what "expected for age" means.
-- **"Years" are clock-relative.** A deviation of +5 under one clock is not comparable to +5 under a different clock, training set, or feature space.
-- **Upstream preprocessing matters — a lot.** As in any clock application scenario, application and/or comparison of aging clock results across cohorts is not trivial. ClockSHAP assumes the feature matrix is analysis-ready and consistent with the clock (transformations, normalization, feature matching, batch correction, QC). Mismatched feature spaces or batch effects produce technical, not biological, deviations. The importance of this step cannot be overstated.
-
-For broader discussion of calibration, age bias, and computational challenges in clock analysis, see [Epigenetic ageing clocks: statistical methods and emerging computational challenges](https://doi.org/10.1038/s41576-024-00807-w) (Nature Reviews Genetics).
-
----
+-   **Deviation is reference-dependent.** The choice of reference cohort
+    defines what “expected for age” means.
+-   **“Years” are clock-relative.** A deviation of +5 under one clock is
+    not comparable to +5 under a different clock, training set, or
+    feature space.
+-   **Upstream preprocessing matters a lot.** ClockSHAP assumes the
+    feature matrix is analysis-ready and consistent with the clock
+    (transformations, normalization, feature matching, batch correction,
+    QC). Mismatched feature spaces or batch effects produce technical,
+    not biological, deviations. The importance of this step cannot be
+    overstated.
+
+For broader discussion of calibration, age bias, and computational
+challenges in clock analysis, see [Epigenetic ageing clocks: statistical
+methods and emerging computational
+challenges](https://doi.org/10.1038/s41576-024-00807-w) (Nature Reviews
+Genetics).
+
+------------------------------------------------------------------------
 
 ## Status
 
diff --git a/data-raw/make_clockshap_example.R b/data-raw/make_clockshap_example.R
new file mode 100644
index 0000000..0b1abd2
--- /dev/null
+++ b/data-raw/make_clockshap_example.R
@@ -0,0 +1,88 @@
+## Generates the bundled synthetic example dataset `clockshap_example`.
+##
+## Run with the package loaded:
+##   devtools::load_all(".")
+##   source("data-raw/make_clockshap_example.R")
+##
+## Produces data/clockshap_example.rda.
+## The data are entirely synthetic; gene symbols are used only for legibility.
+
+set.seed(42)
+
+genes <- c("CDKN2A", "GDF15", "EDA2R", "IL6", "LMNB1", "ELN", "COL1A1", "FOXO3")
+n     <- 120
+age   <- round(runif(n, 30, 85), 1)
+
+age_slope <- c(
+  CDKN2A = 0.040, GDF15 = 0.038, EDA2R = 0.042, IL6 = 0.036,
+  LMNB1 = -0.039, ELN = -0.037, COL1A1 = 0.035, FOXO3 = -0.034
+)
+baseline <- c(
+  CDKN2A = 3, GDF15 = 4, EDA2R = 2, IL6 = 2.5,
+  LMNB1 = 8, ELN = 7, COL1A1 = 5, FOXO3 = 6
+)
+noise_sd <- c(
+  CDKN2A = 0.80, GDF15 = 0.85, EDA2R = 0.80, IL6 = 0.90,
+  LMNB1 = 0.80, ELN = 0.85, COL1A1 = 0.85, FOXO3 = 0.80
+)
+
+X <- sapply(
+  genes,
+  function(g) baseline[g] + age_slope[g] * age + rnorm(n, 0, noise_sd[g])
+)
+rownames(X) <- sprintf("S%03d", seq_len(n))
+colnames(X) <- genes
+names(age) <- rownames(X)
+
+ref_trend <- vapply(genes, function(g) coef(lm(X[, g] ~ age)), numeric(2))
+age_expected <- function(a) ref_trend[1, ] + ref_trend[2, ] * a
+
+mu     <- colMeans(X)
+sigma  <- apply(X, 2, sd)
+Z      <- scale(X, center = mu, scale = sigma)
+lambda <- 40
+beta   <- as.numeric(solve(
+  crossprod(Z) + lambda * diag(length(genes)),
+  crossprod(Z, age - mean(age))
+))
+names(beta) <- genes
+alpha <- mean(age)
+
+clock <- linear_clock(alpha = alpha, beta = beta, mu = mu, sigma = sigma)
+
+shift_genes <- c("CDKN2A", "GDF15", "IL6")
+shift_z     <- c(0.7, 0.7, 0.6)
+
+X["S001", ] <- age_expected(age["S001"])
+X["S001", shift_genes] <- X["S001", shift_genes] + shift_z * sigma[shift_genes]
+X["S002", ] <- age_expected(age["S002"])
+X["S002", shift_genes] <- X["S002", shift_genes] - shift_z * sigma[shift_genes]
+
+clockshap_example <- list(features = X, age = age, clock = clock)
+
+if (requireNamespace("usethis", quietly = TRUE)) {
+  usethis::use_data(clockshap_example, overwrite = TRUE)
+} else {
+  dir.create("data", showWarnings = FALSE)
+  save(clockshap_example,
+       file = file.path("data", "clockshap_example.rda"),
+       compress = "xz",
+       version = 2)
+}
+
+ref <- fit_reference_profile(features = as.data.frame(X), age = age)
+cs  <- clockshap(features = X, age = age, clock = clock, reference = ref)
+
+cat(sprintf(
+  "predicted %.1f-%.1f | cor(pred,age)=%.2f | slope=%.2f | dev sd=%.1f\n",
+  min(cs$predicted), max(cs$predicted), cor(cs$predicted, age),
+  coef(lm(cs$predicted ~ age))[2], sd(cs$deviation)
+))
+cat(sprintf(
+  "S001 deviation = %+.1f | S002 deviation = %+.1f\n",
+  cs$deviation[1], cs$deviation[2]
+))
+cat(sprintf(
+  "additivity max|rowSums(phi) - deviation| = %.1e\n",
+  max(abs(rowSums(cs$phi) - cs$deviation))
+))
diff --git a/data/clockshap_example.rda b/data/clockshap_example.rda
new file mode 100644
index 0000000..0eaff34
Binary files /dev/null and b/data/clockshap_example.rda differ
diff --git a/man/clockshap_example.Rd b/man/clockshap_example.Rd
new file mode 100644
index 0000000..b848c99
--- /dev/null
+++ b/man/clockshap_example.Rd
@@ -0,0 +1,50 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{clockshap_example}
+\alias{clockshap_example}
+\title{Synthetic ClockSHAP example dataset}
+\format{
+A list with three elements:
+\describe{
+\item{\code{features}}{A 120 x 8 numeric matrix of synthetic gene
+expression. Rows are samples (\code{S001}-\code{S120}); columns are genes
+(\code{CDKN2A}, \code{GDF15}, \code{EDA2R}, \code{IL6}, \code{LMNB1},
+\code{ELN}, \code{COL1A1}, \code{FOXO3}).}
+\item{\code{age}}{Numeric vector of length 120, chronological age in years.}
+\item{\code{clock}}{A \code{linear_clock} object (intercept, per-gene
+weights, centering and scaling) fit by ridge regression so predicted age
+tracks chronological age with realistic regression to the mean
+(slope < 1).}
+}
+}
+\usage{
+clockshap_example
+}
+\description{
+A small, fully synthetic dataset used in the package examples and vignette.
+It contains a cohort of 120 samples with eight genes carrying genuine linear
+age trends, plus a ridge-fit linear clock that predicts age from those genes.
+Two samples are planted deviators: \code{S001} reads older than expected for
+its age and \code{S002} reads younger, both driven by a shared
+senescence / inflammaging axis (\code{CDKN2A}, \code{GDF15}, \code{IL6}).
+}
+\details{
+The data are entirely synthetic and carry no real biological signal;
+gene symbols are used only to make the example legible. Regenerated by
+\code{data-raw/make_clockshap_example.R}.
+}
+\examples{
+ref <- fit_reference_profile(
+  features = as.data.frame(clockshap_example$features),
+  age      = clockshap_example$age
+)
+cs <- clockshap(
+  features  = clockshap_example$features,
+  age       = clockshap_example$age,
+  clock     = clockshap_example$clock,
+  reference = ref
+)
+cs$deviation
+}
+\keyword{datasets}
diff --git a/man/figures/README-quick-start-1.png b/man/figures/README-quick-start-1.png
new file mode 100644
index 0000000..384fe71
Binary files /dev/null and b/man/figures/README-quick-start-1.png differ
diff --git a/man/figures/clockshap-overview.png b/man/figures/clockshap-overview.png
new file mode 100644
index 0000000..7d54d5b
Binary files /dev/null and b/man/figures/clockshap-overview.png differ
diff --git a/man/figures/readme-worked-example-waterfall.png b/man/figures/readme-worked-example-waterfall.png
deleted file mode 100644
index 7bea581..0000000
Binary files a/man/figures/readme-worked-example-waterfall.png and /dev/null differ
diff --git a/tests/testthat/test-clockshap-core.R b/tests/testthat/test-clockshap-core.R
index 0f8d0a3..46cb208 100644
--- a/tests/testthat/test-clockshap-core.R
+++ b/tests/testthat/test-clockshap-core.R
@@ -452,6 +452,33 @@ test_that("fit_reference_profile recovers coefficients for perfectly linear data
 })
 
 
+test_that("clockshap_example runs end to end", {
+
+  data(clockshap_example, package = "ClockSHAP")
+
+  expect_named(clockshap_example, c("features", "age", "clock"))
+  expect_equal(dim(clockshap_example$features), c(120, 8))
+  expect_identical(names(clockshap_example$age), rownames(clockshap_example$features))
+  expect_s3_class(clockshap_example$clock, "linear_clock")
+
+  ref <- fit_reference_profile(
+    features = as.data.frame(clockshap_example$features),
+    age = clockshap_example$age
+  )
+  cs <- clockshap(
+    features = clockshap_example$features,
+    age = clockshap_example$age,
+    clock = clockshap_example$clock,
+    reference = ref
+  )
+
+  expect_gt(cor(cs$predicted, clockshap_example$age), 0.8)
+  expect_equal(unname(rowSums(cs$phi)), unname(cs$deviation), tolerance = 1e-10)
+  expect_gt(cs$deviation[1], 0)
+  expect_lt(cs$deviation[2], 0)
+})
+
+
 test_that("Accessors return the expected components", {
 
   features <- matrix(rnorm(12), nrow = 4, ncol = 3)
diff --git a/vignettes/clockshap-introduction.Rmd b/vignettes/clockshap-introduction.Rmd
index 36bb479..0cd7a2d 100644
--- a/vignettes/clockshap-introduction.Rmd
+++ b/vignettes/clockshap-introduction.Rmd
@@ -1,125 +1,255 @@
 ---
-title: "ClockSHAP Introduction"
+title: "Introduction to ClockSHAP"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{clockshap-introduction}
+  %\VignetteIndexEntry{Introduction to ClockSHAP}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---
 
-```{r setup, include=FALSE}
-library(ClockSHAP)
-
+```{r setup, include = FALSE}
 knitr::opts_chunk$set(
-  collapse = TRUE,
-  comment = "#>"
+  collapse  = TRUE,
+  comment   = "#>",
+  fig.width = 7,
+  fig.height = 4,
+  out.width = "100%"
 )
+library(ClockSHAP)
 ```
 
+## Overview
 
-## Introduction
+ClockSHAP interprets the output of a linear aging clock. Instead of asking which
+features push a predicted age up or down in absolute terms, it answers:
 
-ClockSHAP is a method for interpreting *biological age deviation* produced by
-linear aging clocks.
+> **Why is this sample predicted to be biologically older or younger than expected for its chronological age?**
 
-Rather than explaining which features increase or decrease a predicted age,
-ClockSHAP explains **why a sample appears biologically older or younger than
-expected for its chronological age**, relative to a reference aging trajectory.
+For a linear clock it provides an exact, closed-form decomposition of that gap
+into per-feature contributions. This vignette walks through the full workflow on
+the package's bundled synthetic dataset. The data, the core objects, the
+decomposition, the accessors, and how the result relates to the more familiar
+age-acceleration metrics.
 
-The method provides an exact, closed-form decomposition for linear models.
+## The method in brief
 
----
+A linear clock predicts age, $\hat{y}$, from features. ClockSHAP defines the
+**age-matched deviation** as the gap between that prediction and the age-matched
+expectation under a reference:
 
-## Conceptual overview
+$$\Delta_i = \hat{y}_i - \mathrm{E}[\hat{y} \mid \mathrm{Age}_i]$$
 
-Consider a linear aging clock trained on reference data:
+The expectation comes from an **age-conditioned reference profile** — a
+per-feature linear age trend fit on a reference cohort:
 
-\[
-\hat{y} = \alpha + \sum_k \beta_k \tilde{X}_k
-\]
+$$\mathrm{E}[X_k \mid \mathrm{Age}] = \gamma_{0,k} + \gamma_{1,k}\,\mathrm{Age}$$
 
-ClockSHAP defines deviation as:
+Each feature's contribution is its clock weight $w_k$ times its distance from
+that age-expected value, and for a linear clock the contributions reconstruct the
+deviation exactly:
 
-\[
-\Delta = \hat{y} - E[\hat{y} \mid \text{Age}]
-\]
+$$\phi_{i,k} = w_k\big(x_{i,k} - \mathrm{E}[X_k \mid \mathrm{Age}_i]\big), \qquad \sum_k \phi_{i,k} = \Delta_i$$
 
-and decomposes this deviation into per-feature contributions:
+In practice features are standardized to the clock's training distribution before
+the decomposition is computed.
 
-\[
-\phi_k = \beta_k \left( \tilde{X}_k - E[\tilde{X}_k \mid \text{Age}] \right)
-\]
+## The example data
 
+The package ships a small, fully synthetic dataset, `clockshap_example`, so every
+example runs without external data. It is a list with three pieces:
 
-such that:
+```{r data-structure}
+str(clockshap_example, max.level = 1)
+```
 
-\[
-\sum_k \phi_k = \Delta
-\]
+`features` is a 120 × 8 matrix of synthetic "gene" values (rows `S001`–`S120`),
+each gene carrying a genuine linear age trend plus noise; `age` is the
+chronological age; and `clock` is a fitted clock object. Two samples are planted
+deviators: **`S001`** reads older than expected for its age and **`S002`** reads
+younger, both driven by a shared senescence / inflammaging axis (`CDKN2A`,
+`GDF15`, `IL6`). The data carry no real biological signal and the gene symbols are
+only there to make the output legible.
 
-where \(\tilde{X}_k\) denotes features standardized using the clock’s training mean and standard deviation.
+## The building blocks: clock and reference
 
-This decomposition is exact for linear clocks.
+ClockSHAP works with two specification objects.
 
----
+A **`linear_clock`** holds the intercept, the per-feature weights, and the
+training mean and standard deviation used to standardize features. You construct
+one with `linear_clock()`; here we inspect the bundled clock:
+
+```{r clock}
+clock <- clockshap_example$clock
+class(clock)
+round(clock$beta, 2)   # per-gene weights learned by the clock
+```
 
-## A minimal example
+A **`reference_profile`** holds, for each feature, the intercept and slope of its
+age trend ($\gamma_0$, $\gamma_1$). You can build one directly with
+`reference_profile(gamma0, gamma1)`, but usually you fit it from a cohort with
+`fit_reference_profile()`, which regresses each feature on age:
 
-This example uses synthetic data to demonstrate the full ClockSHAP workflow.
+```{r reference}
+ref <- fit_reference_profile(
+  features = as.data.frame(clockshap_example$features),
+  age      = clockshap_example$age
+)
+class(ref)
+round(ref$gamma1, 4)   # per-gene age slopes recovered from the cohort
+```
 
-```{r}
-set.seed(1)
+The non-zero slopes confirm the features genuinely move with age, which is what
+makes the age-conditioned expectation meaningful.
 
-features <- matrix(rnorm(20, mean = 0, sd = 20), nrow = 5, ncol = 4)  # sd was 1 → now 20
-colnames(features) <- paste0("F", 1:4)
+## Running ClockSHAP
 
-age <- c(40, 50, 60, 70, 80)
+`clockshap()` ties the pieces together: features, ages, a clock, and a reference.
+It standardizes the features with the clock's parameters, builds the age-matched
+expectation from the reference, and returns the decomposition.
 
-clock <- linear_clock(
-  alpha = 10,
-  beta  = setNames(runif(4, min = -2, max = 2), colnames(features)), # wider than 0..1
-  mu    = setNames(rep(0, 4), colnames(features)),
-  sigma = setNames(rep(1, 4), colnames(features))
+```{r run}
+cs <- clockshap(
+  features  = clockshap_example$features,
+  age       = clockshap_example$age,
+  clock     = clockshap_example$clock,
+  reference = ref
 )
 
-reference <- fit_reference_profile(
-  features = as.data.frame(features),
-  age = age
+cs                # print method: dimensions + deviation range
+summary(cs)$deviation
+```
+
+`clockshap()` also enforces its core invariant internally — it will error rather
+than return a result where the contributions fail to sum to the deviation.
+
+## The output: structure and accessors
+
+A `clockshap` object is a list with five elements, each retrievable either by
+`$` or with an accessor function:
+
+| field | accessor | meaning |
+|---|---|---|
+| `predicted` | `clockshap_predicted()` | clock prediction $\hat{y}$ |
+| `expected`  | `clockshap_expected()`  | age-matched expectation $\mathrm{E}[\hat{y}\mid\mathrm{Age}]$ |
+| `deviation` | `clockshap_deviation()` | $\Delta = \text{predicted} - \text{expected}$ |
+| `phi`       | `clockshap_phi()`       | per-feature contribution matrix (samples × features) |
+| `age`       | —                       | the input chronological ages |
+
+```{r results}
+head(data.frame(
+  age       = cs$age,
+  predicted = round(clockshap_predicted(cs), 1),
+  expected  = round(clockshap_expected(cs), 1),
+  deviation = round(clockshap_deviation(cs), 1)
+))
+```
+
+The contributions in `clockshap_phi()` sum, row by row, exactly to the deviation:
+
+```{r additivity}
+all.equal(rowSums(clockshap_phi(cs)), clockshap_deviation(cs), tolerance = 1e-10)
+```
+
+## Cohort view
+
+Plotting each sample's predicted age against its chronological age reproduces the
+cohort-level picture. The black line is the age-matched expectation; a sample's
+vertical distance from it is its deviation. The two planted deviators sit clearly
+off the line.
+
+```{r cohort, fig.alt = "Predicted vs chronological age across the synthetic cohort, with two deviators highlighted"}
+library(ggplot2)
+
+cohort <- data.frame(
+  age       = cs$age,
+  predicted = clockshap_predicted(cs),
+  role      = ifelse(rownames(clockshap_example$features) == "S001", "S001 (older)",
+               ifelse(rownames(clockshap_example$features) == "S002", "S002 (younger)",
+                      "cohort"))
 )
 
-cs <- clockshap(features, age, clock, reference)
-cs
+ggplot(cohort, aes(age, predicted)) +
+  geom_point(data = subset(cohort, role == "cohort"),
+             colour = "grey70", alpha = 0.7) +
+  geom_smooth(method = "lm", se = FALSE, colour = "black", linewidth = 0.6) +
+  geom_point(data = subset(cohort, role != "cohort"),
+             aes(colour = role), size = 3) +
+  scale_colour_manual(values = c("S001 (older)" = "#B2182B",
+                                 "S002 (younger)" = "#2166AC")) +
+  labs(x = "Chronological age (years)", y = "Predicted age (years)", colour = NULL) +
+  theme_minimal()
+```
+
+## Age-matched deviation vs. AAA and RAA
 
+It helps to see how the deviation relates to two metrics common in the
+aging-clock literature:
+
+- **Absolute age acceleration (AAA):** $\hat{y} - \mathrm{Age}$.
+- **Relative age acceleration (RAA):** the residuals of `lm(predicted ~ age)`.
+
+```{r metrics}
+age  <- cs$age
+pred <- clockshap_predicted(cs)
+
+metrics <- data.frame(
+  AAA       = pred - age,
+  RAA       = as.numeric(residuals(lm(pred ~ age))),
+  deviation = clockshap_deviation(cs)
+)
 ```
 
-## Inspecting results
-ClockSHAP returns a structured object with accessor functions.
+Because this clock has a slope below 1, AAA is **age-biased** — systematically
+positive for younger samples and negative for older ones — while RAA and the
+deviation are not:
 
-```{r}
-clockshap_deviation(cs)
-clockshap_phi(cs)
+```{r metrics-cor}
+round(sapply(metrics, function(z) cor(z, age)), 3)
 ```
 
-By construction, per-feature contributions sum to the total deviation for each
-sample.
+And here the reference was fit on the same cohort we are analyzing. In that
+self-reference case the age-matched deviation is **exactly** relative age
+acceleration:
 
-```{r}
-all.equal(rowSums(clockshap_phi(cs)), clockshap_deviation(cs), tolerance = 1e-10)
+```{r metrics-equal}
+all.equal(metrics$deviation, metrics$RAA, tolerance = 1e-9)
 ```
 
-## Waterfall visualization
+This is not a coincidence. For a linear clock, passing per-feature age trends
+through the clock reproduces the OLS fit of predicted age on age, so
+`predicted - expected` equals the regression residual. ClockSHAP becomes distinct
+from RAA — and starts carrying new biological signal — only when the reference is
+a *different*, biologically meaningful cohort (e.g., a tissue-of-origin
+reference). What it always adds, even in the self-reference case, is the
+**per-feature attribution** of that deviation, which neither AAA nor RAA provides.
+
+## Per-sample decomposition
+
+`plot_clockshap_waterfall()` shows, for a single sample, how each feature carries
+the prediction from the age-matched expectation to the predicted age. The bars
+sum to the deviation. Compare the two planted deviators:
+
+```{r waterfall-pos, fig.height = 3.4, fig.alt = "Waterfall for S001, an older-than-expected sample"}
+plot_clockshap_waterfall(cs, sample = "S001", top_n = 6)
+```
 
-```{r}
-plot_clockshap_waterfall(cs, sample = 1, top_n = 3, show_age_labels = FALSE, show_delta_label = FALSE )
+```{r waterfall-neg, fig.height = 3.4, fig.alt = "Waterfall for S002, a younger-than-expected sample"}
+plot_clockshap_waterfall(cs, sample = "S002", top_n = 6)
 ```
 
+`S001` reads older because its senescence / inflammaging genes sit above their
+age-expected levels; `S002` reads younger for the mirror-image reason. In both,
+the labeled contributions add up to the age-matched deviation.
+
 ## Scope and assumptions
-ClockSHAP is defined for linear aging clocks and requires:
 
-- linear feature contributions
-- a reference age-feature relationship
-- consistent feature definitions between clock and reference
+ClockSHAP is defined for linear aging clocks and requires:
 
-The method does not train clocks or impose assumptions about the source of
-reference data. Any dataset may be used to define a reference aging trajectory.
+- linear, additive feature contributions (a `linear_clock`),
+- an age–feature reference relationship (a `reference_profile`),
+- consistent feature definitions between the clock and the reference.
 
+It does not train clocks, and it does not assume a particular source of reference
+data — any cohort can define the reference aging trajectory. Because the choice of
+reference determines what "expected for age" means, deviations should always be
+interpreted relative to it.