Age-Related Reference Interval
Many laboratory analytes change continuously with age — haemoglobin, creatinine, IGF-1, and bone markers are well-known examples. A single reference interval applied across all ages will be too wide in some age groups and too narrow in others, leading to misclassification of both healthy individuals and patients. The solution is to model the reference limits as smooth continuous functions of age rather than computing a separate interval within each age stratum.
The age-stratified approach — dividing the age range into discrete classes such as 20–30, 30–40, 40–50 years and establishing a reference interval within each — is intuitively simple but has three practical disadvantages. First, each stratum requires at least 120 subjects for non-parametric confidence intervals, so the total sample size grows with the number of classes. Second, the reference limits are discontinuous at the class boundaries: a patient aged 39 may receive a different reference limit from one aged 40 despite no underlying biological discontinuity. Third, where to place the boundaries is subjective. The regression-based approach described here avoids all three problems: it uses all subjects simultaneously, produces smooth continuous limits, and — because the model borrows strength across the full age range — needs fewer subjects in total.
The Wright & Royston Method
MedCalc implements the approach of Altman (1993), Altman & Chitty (1993), and Wright & Royston (1997). The core idea is straightforward: fit a weighted polynomial regression of the analyte on age to estimate mean(age), then fit a second weighted polynomial regression on the absolute residuals to estimate SD(age). The 95% reference limits at any age follow directly as mean(age) ± 1.96 × SD(age). The six-step procedure is detailed below.
Step 1: Transformation
If the distribution of measurements shows age-dependent skewness, a logarithmic or Box-Cox power transformation is applied before fitting. The Box-Cox formula is:
$$ x^{(\lambda)} = \begin{cases} \dfrac{(x + c)^\lambda - 1}{\lambda} & \lambda \neq 0 \\[6pt] \log(x + c) & \lambda = 0 \end{cases} $$where $c \geq 0$ is a shift parameter that handles zero or negative values, and $\lambda$ is estimated from the data by minimising age-dependent skewness in the residuals. After transformation, the normality of residuals is reassessed at each age. If no transformation achieves an approximately symmetric distribution, the non-parametric reference interval method described in the Reference interval chapter should be considered for each age stratum.
Step 2: Modelling the Mean
The transformed measurements are modelled on age using weighted polynomial regression:
$$ \widehat{\text{mean}}(\text{age}) = b_0 + b_1 f_1(\text{age}) + b_2 f_2(\text{age}) + \cdots $$The user selects which polynomial terms to include. Available terms include integer powers of age (age, age², age³, …) and the logarithm of age (log(age)). Specifying powers {0, 1, 2} produces:
$$ \widehat{\text{mean}}(\text{age}) = b_0 + b_1 \log(\text{age}) + b_2 \cdot \text{age} + b_3 \cdot \text{age}^2 $$A polynomial of degree 1 in age (a straight line) is often adequate for analytes that change monotonically. Adding a squared term allows for a curve that rises and then falls, or vice versa. The log(age) term is appropriate when the rate of change is rapid in early life and slows in adulthood.
Step 3: Residuals
The residual for each subject is the difference between the observed (transformed) value and the fitted mean at that subject's age:
$$ r_i = y_i - \widehat{\text{mean}}(\text{age}_i) $$These residuals represent the individual deviations from the age-dependent mean and carry the information needed to estimate how variability changes with age.
Step 4: Modelling the Standard Deviation
The absolute residuals, multiplied by $\sqrt{\pi/2}$, are modelled on age using a second weighted polynomial regression to produce SD(age). The scaling factor corrects for the systematic underestimation inherent in the mean absolute deviation: for a normally distributed variable with standard deviation $\sigma$, the expected value of $|r|$ is $\sigma \sqrt{2/\pi}$, so multiplying by $\sqrt{\pi/2}$ restores the correct scale and gives an approximately unbiased estimator of $\sigma$:
$$ \widehat{\text{SD}}(\text{age}) = c_0 + c_1 g_1(\text{age}) + c_2 g_2(\text{age}) + \cdots $$The polynomial terms for the SD model are chosen independently from those for the mean model. In many analytes the standard deviation increases with the mean (proportional variability), in which case a log transformation of the measurements followed by a constant SD model is equivalent to modelling the original data with SD proportional to the mean. When the SD varies substantially with age — as in the illustrative example in the figure, where older subjects show greater variability — including an age term in the SD model is necessary.
Step 5: Reference Limits
For any age of interest, the reference limits are:
$$ \text{Lower RL}(\text{age}) = \widehat{\text{mean}}(\text{age}) - 1.96 \times \widehat{\text{SD}}(\text{age}) $$ $$ \text{Upper RL}(\text{age}) = \widehat{\text{mean}}(\text{age}) + 1.96 \times \widehat{\text{SD}}(\text{age}) $$If a transformation was applied in Step 1, the limits are back-transformed to the original measurement scale. Because back-transformation is a non-linear operation, the resulting reference limits are asymmetric around the mean — they correctly reflect the skewed shape of the original distribution rather than forcing a symmetric interval.
Step 6: Model Validation
After fitting, each subject receives a standardised residual (z-score):
$$ z_i = \frac{r_i}{\widehat{\text{SD}}(\text{age}_i)} $$A well-fitted model produces z-scores that are:
- Approximately normally distributed across all ages (Shapiro-Francia test).
- Randomly scattered around zero with no systematic trend across the age range — a fan-shaped pattern indicates that the SD model is misspecified; a curved pattern indicates that the mean model needs additional polynomial terms.
- About 5% lying beyond ±1.645 (the expected fraction outside a 95% interval).
If the z-scores reveal a pattern, the polynomial structure of the mean or SD model is revised and the procedure is repeated. Model adequacy should be confirmed before the reference limits are reported.
Confidence Intervals for the Reference Limits
Because the reference limits at each age depend on two regression models, no simple closed-form expression for their standard errors exists. MedCalc therefore uses the bootstrap resampling procedure of Wright & Royston (1997): subjects are resampled with replacement, the entire six-step procedure is repeated for each replicate, and the 5th and 95th percentiles of the bootstrap distribution of each reference limit form the 90% confidence interval at that age. The resulting CI curves widen where the reference sample is sparse and narrow where it is dense.
Selecting Polynomial Terms
There is no automatic procedure for choosing which polynomial terms to include in the mean and SD models; the choice requires biological judgement and iterative model checking. Start with the simplest biologically plausible model — often a single linear term in age, or log(age) if the relationship is more rapid early in life — and check the z-score scatter plot. Add complexity only if the residuals reveal a systematic pattern that the simpler model cannot explain. A model with many polynomial terms may fit the sample well but will be unstable out of sample; the reference limits should vary smoothly and the fitted curves should be biologically interpretable.
MedCalc's Age-related reference interval procedure fits the mean and SD models with user-specified polynomial terms, produces the continuous reference limit curves with 90% bootstrap confidence intervals, and generates the z-score diagnostic plot for model validation.