Direct, age- and gender-specific reference intervals: applying a modified M-estimator of the Yeo-Johnson transformation to clinical real-world data

Direct parametric reference interval estimation

The aforementioned guideline proposed the estimation of RIs directly from the observed values or from a mathematical “normalisation” of said values, if either of them follow a Gaussian distribution [4]. The Box-Cox (BC) transformation family was proposed to find the most appropriate transformation for the observed distribution to a Gaussian one [8]. Under the mild assumption that the underlying data are non-uniformly distributed, it has been an established transformation method for handling data in direct RI studies [9], [10], [11], [12], [13]. The BC transformation is used to transform the data through:

The transformation parameter λ (lambda) can be estimated using maximum likelihood estimation (MLE), however this estimation can be highly affected by outlying observations from the main bulk of the distribution, frequent in real-world routine data. Various methods have been suggested to address the influence of outliers before transformation, such as evaluating the Dixon’s Q statistic, the Grubbs G statistic or the Tukey quantile range on suspected values [14], [15], [16]. These methods have a variety of limitations, such as limited power with large sample sizes, where the location and prevalence of outliers is not known a priori. Further, the distinction between a data point from the tail and an outlier in skewed distributions is even harder, biasing RI estimates [17, 18]. The re-evaluation of suspected outliers can help reduce their influence on the transformation to normality.

Aim

In this work, we propose using a robust transformation algorithm to achieve central normality in univariate blood analyte test distributions from routine laboratory data from which RIs can be estimated parametrically. Using an established reweighted M-estimator (rewME) approach for the Box-Cox and Yeo-Johnson transformations, the accuracy of the estimated transformation parameter lambda can be improved. This renders the direct parametric estimation of RIs from the resulting transformed distribution more robust to outliers that may be present in routine patient data. We can thus directly infer clinically relevant RIs as “expectation ranges” from gender- and age-stratified reference populations.

Reweighted reference interval estimation

The proposed algorithm for robustly estimating RIs assumes that the reference population from routine patient data contains a majority of non-pathological test result values. The distribution can be closely approximated with a Box-Cox (BC) or a Yeo-Johnson (YJ) transformed Gaussian normal distribution. The YJ transformation is a generalisation of the BC transformation and also employs a single transformation parameter lambda as follows:

Patient data from the clinical routine contain an unknown prevalence of outliers, i.e. outlying observations that deviate from the majority of values in the distribution, that can skew the optimal estimated transformation if the transformation parameter lambda is estimated by MLE. However, the implemented algorithm for direct RIs (rewRI) uses a reweighted M-estimator (rewME) for the transformation parameter lambda to achieve central normality, as established by Raymaekers et al. [19]. The implemented algorithm adaptively adjusts the transformation parameter lambda with associated weights, incorporates machine learning principles to minimise the impact of outliers and achieve a more robust, data-driven normalisation of the real-world data. With this, the influence of outliers on the direct RI inference can be reduced compared to the traditional one-step Box-Cox transformation approaches. This is achieved by minimising a robust M-estimate criterion for the parameter λ for the BC or the YJ transformation. In robust linear regression, M-estimators are usually computed by iteratively applying classical weighted regression, i.e. using the standard maximum likelihood estimator incorporating the weights, where the weights are recalculated in each iteration step. A similar scheme is adapted here to estimate the parameters for the Box-Cox or YJ transformation and the parameters of the normal distribution of the transformed non-pathological values. In the scheme proposed here, binary instead of the continuous weights in robust regression are used. A detailed description of the heuristic algorithm can be found in the Supplementary material. The resulting optimal λ can be used to transform the data to robustly achieve a close approximation of a Gaussian distribution, from which the appropriate (1 − α) RI with the lower reference limit (LL) and the upper reference limit (UL) can be estimated. Using bootstrap, confidence intervals (CIs) are estimated to infer precision by drawing n-sized samples with replacement from the reference population, repeating this methodology for a reasonable number of times (e.g. 400), and calculating CIs from the quantiles of the bootstrap distribution. The full method has been implemented in an R package called rewRI. The package supports the fully automated estimation of direct RIs from clinical cohort data, with an accompanying graphical output to further help to assess the accuracy of the inferred estimates.

Method evaluation on simulated data

A direct parametric RI approach from clinical data tries to leverage the mean and standard deviation of the estimated Gaussian distribution from the underlying data. Both of these metrics can become inaccurate if the presence of outliers and residual skewness is not addressed properly. Generally, the presence of either erroneous or significantly differing (potentially pathological) test results in routine patient data can vary. To evaluate the performance of the “rewRI”, a rewME of the best fitting BC/YJ transformation to estimate RIs, different approaches have been performed to evaluate the robust RI estimation.

Firstly while minimising the rewME for the transformation parameter λ, local convergence can be obtained with the resulting λ̂ estimate not necessarily being the global minimum. Especially in lower sample sizes, where the influence of individual data points on the overall variance is much bigger, this can lead to substantially biased estimates. To assess the accuracy of the rewME of the optimal lambda for the reweighted YJ and BC transformation algorithms in absence of significant outliers, simulated datasets were created based on a normal distribution X∼N(µ=5, σ²=1), which were transformed with the inverse of the corresponding transformation using a grid of lambda values λ ∈ {0, 0.001, 0.1, 0.2, 0.5, 1} at varying sample sizes n [18]. The distributions were then used to estimate λ̂ by bootstrapping the algorithm b=400 times with the reweighting CI τ=0.99 across the chosen grid of expected lambda values for both the BC and YJ transformation.

Secondly to evaluate the method’s robustness in presence of outliers, different prevalences of outlying observations have been simulated. Again, simulated datasets were created based on the inversely transformed distribution based on X∼N(µ=5, σ²=1) using the grid of lambda values λ ∈ {0, 0.001, 0.1, 0.2, 0.5, 1} of various sample sizes n. To introduce outliers in both the lower and upper tail of the distribution (X _left and X _right ), a number of outliers representing 1, 5 and 10 % of the overall sample sizes were randomly drawn. Higher prevalence of outliers (>15 %) was not considered during the simulations, as more than 15 % of outlying observations would significantly alter the overall skewness of the data, making it impossible to robustly optimise the criterion [19]. Outliers were added before the inverse transformation. The amount of values added to each side was randomly assigned, totalling in the percentage set, and the values were drawn from uniform distributions, as:

These outliers therefore solely contribute to the tails of the inverse transformed distribution, for which an estimate on their effect on the inferred RIs was calculated. For the application in robust RI estimation, this was quantified as the change in the estimated LL and UL, as well as the width of the accompanying CIs. With an increased presence of outliers in the tails of the distributions, the reweighting step of the rewME (τ) was carefully evaluated at various sample sizes relevant for robust direct estimation.

Method evaluation on reference data

The rewRI algorithm with the YJ transformation was further tested on an available reference data, the stratified HCV dataset. The dataset contains test results of 238 female and 377 male participants, where in addition to test laboratory results the person’s age (in years), gender (female and male) and the patient’s eligibility for blood donation (categorical value) is recorded. The non-qualification for blood donation is likely due to the underlying medical condition, which resulted in potentially pathological measurement for some analytes. The HVC dataset is beneficial for evaluating the efficiency of the rewRI approach as we can first estimate 95 % RIs as target ranges directly from the seemingly healthy blood donors by nonparametric quantile bootstrapping of the data. The obtained RIs can be further compared to RIs obtained using an indirect RI estimation method based on truncated QQ-plots [20, 21]. An effective implementation of the latter method is done in the reflimR package, which allows the indirect estimation of RIs from routine patient data [22, 23].

Direct gender- and age-stratified RIs

Finally, the performance of the rewRI algorithm was reviewed on a retrospective dataset of routine patient data from a tertiary care hospital. The inferred RIs can be compared against estimates obtained using a non-parametric quantile bootstrapping approach and from the reflimR algorithm. The data used were collected between January 2014 and December 2022 at the University Hospital of Bern, CH, from adult inpatients, who have given general consent for the re-use of their health data. In the extracted data, each laboratory test result is accompanied with the demographic patient factors (gender and age) as well as factors from the (pre-)analytical phase (device/test kit information), resulting in a highly standardised and interoperable clinical dataset [24]. The dataset has been prepared for the inference of RIs, with only the first measurement per administrative case included during the extraction. The dataset has been anonymized according to the ethical waiver granted by the Bern cantonal Ethics Committee (Req-2020-00630).

Reweighted lambda estimation

The M-estimate of the transformation parameter λ for the reweighted YJ and BC transformation algorithms exhibited high accuracy in the obtained λ due to the focus on obtaining central normality. With a sample size of n=1,000 values, the resulting lambda estimates from the bootstrap replicates retained various degrees of variability for both the reweighted BC and YJ transformation (Figure 1A and C respectively). However, this variability did not significantly hinder the inference of the respective RIs by the rewRI algorithm, as the estimated LL/UL were contained in the CIs of the expected “true” lower limits (LL) and “true” upper limits (UL), or the “true” CIs overlapped with the estimated CIs (Figure 1B and D). In absence of distinct outlying observations, the BC did resolve more accurate RIs from the distributions with various skewness, as the estimated RIs were more closely to their true value (Figure 1B). In contrast, the rewME for the YJ transformation appeared to be better equipped to distinguish between significant positive skewness (λ<0.5) and mild skewness (λ≥0.5) in the transformed distributions compared to the rewME of the BC transformation, as there was no overlap between the obtained 90 % CIs (Figure 1C and A respectively).

Figure 1:

Estimated lambdas of reweighted M-estimator and its influence on the 95 % reference interval inference. (A) Estimated lambdas for the grid of expected lambda values (red: median value, black: 90 % confidence interval) for the reweighted M-estimator for the Box–Cox (BC) transformation. (B) Estimated limits of the 95 % reference interval (LL: lower limit, UL: upper limit) for the various transformed distributions using the rewRI algorithm using the reweighted BC transformation. The expected LL and UL limit (with 90 % CI) are shown as vertical blue lines. The median estimate is shown as a dot point and the 90 % CI as vertical lines. (C) Estimated lambda values for the reweighted M-estimator for the Yeo-Johnson (YJ) transformation. (D) Estimated limits of the 95 % reference interval (LL: lower limit, UL: upper limit) from the rewRI algorithm using the reweighted YJ transformation. Sample size n=1,000, number of bootstrap replicates b=400, τ=0.99.

In larger sample sizes the overall precision of the estimated lambda increases, allowing a continuous distinction between variously skewed distributions for both the YJ and BC transformation. However, in sample sizes below n=1,000, the distinction between mild and more severe skewness in the transformed distribution for both the YJ and BC became less pronounced, as seen in sample sizes of n=200 values, which is around the recommended minimal number for RI inference (Figure 2A and C respectively). Due to the low sample sizes, individual values present in the transformed distributions disproportionately influence the estimation of lambda during the rewME. However the reweighting step during the rewME focuses on achieving central normality, making the RI estimation with the rewRI algorithm still possible at sample sizes of n=200 – although with an increase in the observed CIs (Figure 2B and D). Again, in absence of clear outlying observations the BC did resolve more accurate RIs (i.e. closely to their true values), yet all estimated CIs overlapped with the true CIs (Figure 2B and D).

Figure 2:

Estimated lambdas of reweighted M-estimator and its influence on the 95 % reference interval inference. (A) Estimated lambdas for the grid of expected lambda values (red: median value, black: 90 % confidence interval) for the reweighted M-estimator for the Box–Cox (BC) transformation. (B) Estimated limits of the 95 % reference interval (LL: lower limit, UL: upper limit) for the various transformed distributions using the rewRI algorithm using the reweighted BC transformation. The expected LL and UL limit (with 90 % CI) are shown as vertical blue lines. The median estimate is shown as a dot point and the 90 % CI as vertical lines. (C) Estimated lambda values for the reweighted M-estimator for the Yeo-Johnson (YJ) transformation. (D) Estimated limits of the 95 % reference interval (LL: lower limit, UL: upper limit) from the rewRI algorithm using the reweighted YJ transformation. Sample size n=200, number of bootstrap replicates b=400, τ=0.99.

Presence of outliers

The reweighting step of the rewME (τ) was carefully evaluated at various sample sizes to address the increased presence of outliers in the tails of the distributions. The robustness of the rewRI algorithm with the YJ transformation is demonstrated at simulated sample sizes of 1,000 and 200 (Figure 3). In the presence of >1 % of outliers, the estimates of the LL and UL obtained by the reweighted YJ transformation were more precise than the reweighted BC transformation for the distributions with the various skewnesses. It was generally observed that <10 % outliers, the CI of the reweighting set to τ=0.99 resulted in accurate UL and LL estimates. With 10 % of outliers, and especially in bigger sample sizes (n≥1,000), the estimates obtained with τ=0.95 were more precise. In lower sample sizes (n=200) with 1 % outliers, it was observed that the LL/UL estimates were more accurate with a wider CI during reweighting (τ=0.995).

Figure 3:

Yeo-Johnson transformation: Influence of different amounts (1 , 5 and 10 %) of outliers on the estimation of the 95 % reference limits. The rewRI algorithm with the YJ transformation was bootstrapped 400 times, which resulted in a median reference limit estimate (LL=lower limit, UL=upper limit) (dot point) and a 90 % CI (vertical line). Confidence interval of the reweighting (τ) top: sample size of 1,000 values with 10 (1 %), 50 (5 %), and 100 (10 %) values as outliers in the transformed distributions. Bottom: Sample size of 200 values with 2 (1 %), 10 (5 %), and 20 (10 %) values as outliers in the transformed distributions.

The rewRI algorithm did perform better in larger sample sizes with the rewME for the YJ transformation compared to the BC transformations (Supplementary Figure 1). This could be seen in distributions with a higher prevalence of outliers (>5 %), as lower estimates for both the LL and UL with bigger CIs were obtained with the reweighted BC transformation, resulting in right-shifted RIs. Also a lower precision was observed for both the LL and UL.

Evaluation on reference data

The rewRI algorithm with the YJ transformation was further tested on the stratified HCV dataset (Table 1). In almost all cases, the obtained LL and UL of the 95 % RI estimated with the rewRI algorithm were consistent with the target values, meaning that the CI of the target values included the inferred values. Further, the estimates obtained were comparable with estimates obtained with the reflimR method. The rewRI algorithm performed the best with τ=0.99 during the reweighting step, which labelled and downweighted an average of 2.9 % (IQR: 1.72–3.75 %) of the data as outliers in the defined patient strata. In most considered populations, there was seemingly a clear separation between in- and out-lying observation possible. In rare cases, where a clear separation was not possible during reweighting, this usually affected only the UL. Such is the case in the UL of Alanine Aminotransferase (AST) in male patients, where 14.3 % of the total values in this stratum exceeded the target UL. In contrast, only 9.2 % of the total values in female patients were above the target UL, which enabled the robust estimation of the UL with the rewRI algorithm. Furthermore, the UL obtained for both male and female patients for Gamma-Glutamyl Transferase (GGT) was highly biased, as there was a significant overlap between values from apparently healthy individuals and pathological values, which the rewRI algorithm could not clearly separate during the reweighting, resulting in much higher ULs.

Direct gender- and age-stratified RIs

The rewRI algorithm was finally applied to a set of laboratory analytes from a real-world patient dataset, which was collected during routine monitoring in a tertiary hospital. The patient data for the following analytes was used: Creatinine (LOINC: 14682-9), Leukocytes (LOINC: 6690-2), Sodium (LOINC: 2951-2), and Potassium (LOINC: 2823-3). With a preliminary outlier exclusion, extreme and likely implausible values were excluded from the dataset. This resulted in the various sizes for the gender-stratified reference populations used for the analysis (Table 2).

For each analyte, test results from patients, ages 18–90 years, were stratified by the administrative gender and a specific age (in years). Specifically, RIs were calculated using a sliding window approach for each three-year age band, as established by Zierk et al. [25]. This meant for example that the inferred RI for 57-year-old male patients included data from males ages 56 to 58. Iteratively, the rewRI algorithm for the YJ transformation (400 bootstrap replicas, τ ∈ {0.99, 0.95}) was run on each gender- and age-stratum, resulting in quasi-continuous RIs. To compare the estimated results for each obtained stratum, the analysis was repeated to estimate RIs by a nonparametric quantile approach with 400 bootstrap replicates, as well as estimating RIs by the reflimR method. Estimates obtained with the non-parametric quantile bootstrapping shows where 95 % of the clinical data fall without considering outliers. This essentially provides a safeguard, as estimated reference limits by the rewRI algorithm falling outside would be considered not improved. All estimated reference limits (UL/LL) were shown superimposed on 2D heatmaps (based on kernel density estimates of the age-stratified data used for the analysis [26]) to visualise their position relative to the distribution and concentration of data points across different age groups. The gender- and age-stratified RIs for creatinine as well as leukocytes from the rewRI algorithm were similar yet more precise than reflimR estimates, and they were consistently narrower than the non-parametric ones (Figure 4).

Figure 4:

Gender- and age-stratified reference intervals for creatinine and leukocytes. Inferred estimates for the upper and lower limits of the RIs are provided with a 90 % confidence band for the rewRI algorithm (with τ ∈ {0.99, 0.95}, 400 bootstrap replicas), the reflimR method, and the nonparametric quantile bootstrapping (np-quantile, 400 bootstrap replicas). 2D kernel density estimates with the colour intensity representing the density of available measurement for each age year (x-axis) and the respective measurement level (y-axis).

For creatinine, the rewRI algorithm with τ=0.99 identified an average of 5.3 % of the data as outliers in each stratum, and with τ=0.95 an average of 8.0 % as outliers per stratum. The latter result was comparable to the reflimR method, which labelled an average of 8.7 % of the data as “pathological” values per stratum. Similarly for leukocytes, the rewRI algorithm with τ=0.99 labelled an average of 5.7 % of the data as outliers and with τ=0.95 an average of 6.2 % were labelled per stratum. In this analyte, the reflimR method labelled an average of 6.6 % of the data as pathological values per stratum. In most cases, the estimated results from the rewRI with τ=0.99 were only marginal improvements over the non-parametric quantile bootstrapping. If however rewRI is supplied with τ=0.95, which optimises the central 95 % of the total values, narrower age-and sex-stratified RIs were obtained. As these RIs are optimised for the underlying central 95 % range, they can be used as personalised expectation ranges, highlighting where values from a comparable homogeneous reference population come from.

Age- and gender-adjusted estimates obtained for sodium and potassium inferred by the rewRI algorithm (with τ=0.95) were comparable to the reflimR estimates (Figure 5). For sodium, the rewRI algorithm (τ=0.95) labelled an average of 11.1 % of the data as outliers per stratum. Results obtained by the reflimR method, estimated an average of 19 % of the data as pathological values per stratum. For potassium, the rewRI algorithm (τ=0.95), labelled an average of 11.7 % of the data as outliers. Estimates obtained by the reflimR method labelled an average of 15.8 % of the data as pathological values.

Figure 5:

Gender- and age-stratified reference intervals for sodium and potassium. Inferred estimates for the upper and lower limits of the RIs are provided with a 90 % confidence band for the rewRI algorithm (with τ=0.95, 400 bootstrap replicas), the reflimR, and the nonparametric quantile bootstrapping (np-quantile, 400 bootstrap replicas). 2D kernel density estimates with the colour intensity representing the density of available measurement for each age year (x-axis) and the respective measurement level (y-axis).