Development of next-generation reference interval models to establish reference intervals based on medical data: current status, algorithms and future consideration

Abstract

Evidence derived from laboratory medicine plays a pivotal role in the diagnosis, treatment monitoring, and prognosis of various diseases. Reference intervals (RIs) are indispensable tools for assessing test results. The accuracy of clinical decision-making relies directly on the appropriateness of RIs. With the increase in real-world studies and advances in computational power, there has been increased interest in establishing RIs using big data. This approach has demonstrated cost-effectiveness and applicability across diverse scenarios, thereby enhancing the overall suitability of the RI to a certain extent. However, challenges persist when tests results are influenced by age and sex. Reliance on a single RI or a grouping of RIs based on age and sex can lead to erroneous interpretation of results with significant implications for clinical decision-making. To address this issue, the development of next generation of reference interval models has arisen at an historic moment. Such models establish a curve relationship to derive continuously changing reference intervals for test results across different age and sex categories. By automatically selecting appropriate RIs based on the age and sex of patients during result interpretation, this approach facilitates clinical decision-making and enhances disease diagnosis/treatment as well as health management practices. Development of next-generation reference interval models use direct or indirect sampling techniques to select reference individuals and then employed curve fitting methods such as splines, polynomial regression and others to establish continuous models. In light of these studies, several observations can be made: Firstly, to date, limited interest has been shown in developing next-generation reference interval models, with only a few models currently available. Secondly, there are a wide range of methods and algorithms for constructing such models, and their diversity may lead to confusion. Thirdly, the process of constructing next-generation reference interval models can be complex, particularly when employing indirect sampling techniques. At present, normative documents pertaining to the development of next-generation reference interval models are lacking. In summary, this review aims to provide an overview of the current state of development of next-generation reference interval models by defining them, highlighting inherent advantages, and addressing existing challenges. It also describes the process, advanced algorithms for model building, the tools required and the diagnosis and validation of models. Additionally, a discussion on the prospects of utilizing big data for developing next-generation reference interval models is presented. The ultimate objective is to equip clinical laboratories with the theoretical framework and practical tools necessary for developing and optimizing next-generation reference interval models to establish next-generation reference intervals while enhancing the use of medical data resources to facilitate precision medicine.

Keywords:

• Reference interval
• next-generation
• models
• splines
• indirect sampling technique

Introduction

Laboratory medicine underpins patient care and is used to inform diagnosis, monitor treatment, and predict outcomes. At the heart of this practice are reference intervals (RIs) that are critical for interpreting test results and guiding clinical decisions. Traditional RIs face significant limitations and may lead to inaccurate clinical decisions when demographic variables such as age and sex influence test results.

In recent years, collaborative efforts across disciplines and the emergence of real-world studies have paved the way for advances in the establishment of reference intervals. The advances in computational power have catalyzed a paradigm shift in RI development [1–5], with recent methodologies harnessing the power of these technologies to establish more accurate and demographically tailored RIs [6–11]. Drawing upon the contributions of scholars in the field and our own research, we propose the term “next-generation reference interval models” for the first time. Such advanced models are characterized by a curve function for test results for different age and sex categories and offer dynamically tailored intervals in contrast to traditional static or broad age-range intervals. By auto-selecting the appropriate reference interval based on each patient’s age and sex, they equip clinicians with tools for precise, individualized decision-making. This not only refines the accuracy of disease diagnosis and treatment but also elevates the overall quality of health management.

Despite their potential, the path to widespread adoption of such models has many challenges. First, the process of collecting a robust, representative sample of healthy individuals to establish these intervals is costly, time-consuming, and, for certain types of specimens (e.g. cerebrospinal fluid), practically unfeasible. Moreover, the precursors of such models, known as continuous reference intervals, were tailored predominantly for pediatric cohorts because of variations in their physiological parameters [12,13]. The lack of information system capabilities has hindered the widespread application of these findings in clinical practice. Recognizing these challenges, our team, among others, has leveraged the use of medical data derived from the electronic health record (EHR) to validate the feasibility of indirect sampling techniques for various test parameters and to establish more comprehensive next-generation RIs [14–17].

Zierk et al. [10] established a model to develop RIs for hematology tests using an indirect sampling technique (kosmic algorithm) and smooth spline algorithms, which they referred to as the next-generation reference interval in their publication. Such efforts have set the groundwork for broader adoption and application. Yet, significant hurdles remain: the precise definition of the next-generation reference interval model is still nebulous, there are limited models available, and the methodologies behind their creation and refinement lack clarity and consensus.

This review seeks to demystify the development of next-generation reference interval models, to provide a comprehensive overview of their current status, and to describe their advantages and outline ongoing challenges. We describe the complexities of model development and discuss advanced algorithms and the tools that are integral to their development and optimization. Furthermore, we discuss the potential future of these models in utilizing medical data derived from the EHR with the ultimate aim of providing clinical laboratories with the necessary theoretical and practical knowledge to develop their own next-generation RIs. This would also allow the clinical laboratory to use its own medical data, thereby promoting precision medicine.

Current status of next-generation reference interval models

Definition of next-generation reference interval model

The reference interval (RI) is crucial for clinical diagnosis, treatment monitoring, and prognosis evaluation. Its significance is highlighted by the range of test results observed in healthy individuals, typically represented with a confidence level of 95% (Figure 1). The traditional methodology for setting up RIs is the direct approach in which individuals are chosen from a reference population based on specific criteria and the samples are collected and analyzed for the selected measurands. This process bifurcates into the a priori and a posteriori selection mechanisms. The a priori method involves the selection of individuals who meet pre-defined inclusion criteria, while the a posteriori method includes specimens based on supplementary details like clinical data or other test results, and not all collected samples may be used for further examination [18,19]. While the ideal scenario for the direct method involves random selection from the reference population, in practice, this is not often achieved. The selection of the group for testing is often swayed by factors such as cost and convenience. True randomization, which ensures a wholly representative group, demands substantial planning and resources. Another pivotal aspect of the direct approach is the inclusion of each test result in the statistical assessment, necessitating the exclusion of outliers. This exclusion can influence the intervals that are determined [18,19]. Challenges with direct studies include defining health, the presence of unnoticed diseases, and potential biases from smaller cohorts.

Figure 1. The concept of reference interval.

In comparison, the indirect sampling technique, also termed data mining, taps into existing laboratory data. These methods use statistical tools to parse reference values from large datasets. Critics of this method highlight the limited oversight of pre-analytical conditions and the dependence on statistical tools for filtering out unhealthy subjects. Its supporters appreciate its clinical relevance and its straightforwardness, especially for specific populations like neonates and the elderly. In recent years, the advent of the digital age has increased the interest in real-world studies to derive RIs from mixed datasets of physiological and pathological test results [11,20–22]. Real-world data is information collected from actual, everyday clinical events and interactions as opposed to data generated in controlled or experimental settings. This development reduces the economic cost of establishing RI significantly and addresses the situation in which RIs from other sources are difficult to apply to various test systems and specific populations. However, challenges remain for test parameters heavily influenced by age and sex. Utilizing a single RI or discrete RIs grouped by age results in increased misjudgments that impact clinical decision-making (Figure 2). A model for RIs that offers continuous age- and sex-specific intervals tailored to different age and sex categories has been introduced [10,11].

Figure 2. The limitations of the single reference interval. The blue dots represent individual test results; taking the upper limit of the reference interval as an example, there are a large number of results in the area between the upper limit of a single reference interval and the continuous upper limit calculated by the next-generation reference interval model (the grey shaded part).

Automatically selection of the appropriate reference interval according to the patient’s sex and age facilitates the determination of normal or abnormal test results, aids clinical decision-making and reduces misjudgment rates associated with static RIs. These models are termed “next-generation reference interval models” in this review.

Advantages of and challenges in developing next-generation reference interval models

A next-generation reference interval model embedded in the clinical laboratory information system can provide clinicians with more accurate decision-making tools and enhance convenience for utilization. In terms of methodology, the subjects used in establishing next-generation reference interval models can be acquired through both direct and indirect sampling techniques. Like the establishment of a single reference interval, the next-generation reference interval model based on reference individuals obtained by a direct sampling technique is accurate and reliable. However, our previous studies have found that establishing a single stable reference interval for tests with relatively large variation, such as thyroid stimulating hormone, would require at least 850 samples [23,24]. Thus, development of a next-generation reference interval model by recruiting healthy individuals free of disease would require a vast number of samples to achieve a relatively stable and accurate model. This method is both costly and time-consuming. Moreover, it is challenging to obtain specimens like aqueous humor and cerebrospinal fluid through a direct approach [18,19]. Compared with this, mining massive medical data generated in clinical laboratories provides a new way for establishing next-generation reference intervals and compensates for the limitations of direct sampling techniques. This approach significantly reduces the cost of establishing next-generation reference intervals and simplifies the procedure. Moreover, this method enables the establishment of next-generation reference intervals for specific types of hard-to-obtain specimens. Next-generation reference interval models can also present results in percentiles or Z-score charts [25], enabling normalization across different systems and significantly enhancing convenience in interpreting test results to assist clinical decision-making.

However, challenges still remain. It has not been determined which tests would benefit by having their RIs established using a next-generation reference interval model. Additionally, such a model has been developed for only a few tests. Lastly, given the proliferation of modeling algorithms in the field, along with the many parameters they entail, learning about this domain may be challenging.

Process and algorithms for developing next-generation reference interval models and establishing next-generation reference intervals based on medical data

In recent years, there has been a notable surge in interest regarding next-generation reference interval models. The methodology for formulating such models consists of two fundamental stages: the selection of an appropriate reference population followed by the development of the model. Methods for constructing such reference interval models can be classified into three primary categories based on the selection of the reference population and the subsequent model development. In this paper, they are referred to as category I, II and III methods. Specifically, these are:

• Direct sampling techniques with modeling (category I methods);

• Indirect sampling techniques (1) with modeling (category II) methods; and

• Indirect sampling techniques (2) with curve fitting (category III methods).

Table 1 provides a concise overview of these three methodological approaches, placing special emphasis on category III methods, which have been the subject of more intensive research compared the other categories. These studies have centered around biochemical analytes, hematological parameters, and hormones integral to endocrine metabolism. Several studies [39,40] have also addressed tests in pregnancy.

Table 1. Methods and content of previous studies on establishment of next-generation reference intervals models.

Reference	Subjects	Analytes	Category	Approach type	Modeling reasons	Grouped by sex	CI	Validation (method)
Wilson [26]	Healthy children and adolescents	6 biomarkers, 4 vitamins, minerals, and/or blood proteins, 4 endocrine and/or fertility markers	1	Direct method, LMS model	Clinically relevant	Yes	95% CI	NR
Yan [27]	Healthy children and adolescents	13 biochemical analytes, 8 hematological analytes	1	Direct method, GAMLSS model	Visual examination, decision tree algorithm	Yes	NR	Yes (out of range)
Tybirk [28]	Adults (age ≧21 years)	Glial fibrillary acidic protein	1	Direct method, non-parametric quantile regression	NR	No	95% CI	Yes, (n-1)-fold cross-validation
Hall [29]	Healthy children and adolescents	19 endocrine, fertility and immunochemical markers	1	Direct method, non-parametric quantile regression	NR	Yes	95% CI	NR
Cao [30]	Healthy children and adolescents (≦19 years)	Insulin-like growth factor 1	1	Direct method, LMS model	Clinically relevant	Yes	NR	NR
Wilson [31]	Healthy children and adolescents (≦18 years)	19 hematological analytes	1	Direct method, non-parametric quantile regression	Robustness of next-generation reference intervals models, age partitioning not required	Yes	95% CI	Yes, k-fold cross-validation
Holmes [32]	Healthy children and adolescents(≦18 years)	Testosterone, sex hormone binding globulin, free testosterone, bioavailable testosterone	1	Direct method, non-parametric quantile regression	Robustness of next-generation reference intervals models, age partitioning not required	Yes	95% CI	10-fold cross-validation
Cai [33]	Healthy neonates and children(≦18 years)	Cystatin C	1	Direct method, non-parametric quantile regression	Clinically relevant	Yes	95% CI	NR
Li [34]	Healthy children and adolescents (0 - 19 years)	Alkaline phosphatase, creatinine	1	Direct method, LMS model, GAMLSS model, radial smoothing method, Arzideh method and smoothing spline function	Comparison of algorithms	Yes	NR	Yes (out of range)
Asgari [35]	Healthy children and adolescents (6 months − 19 years)	38 biochemical markers	1	Direct method, non-parametric quantile regression	NR	Yes	95% CI	K-fold cross-validation
Peitzsch [36]	Healthy children and adolescents (age ≦18 years)	Free normetanephrine, metanephrine, 3-methoxytyramine and 3-O-methyldopa	1	Direct method, improved Altman method by Royston and Wright	Clinically relevant	Yes	NR	NR
Vogel [37]	Non-diseased children (age ≦18 years)	High density lipoprotein cholesterol	1	Direct method, combination of LMS model and resampling methods	NR	Yes	95% CI	NR
Ma [38]	Adults (age ≧18 years)	Homocysteine	2	Indirect method, GAMLSS model	Visual examination, variance component model	Yes	95% CI	Yes (far outside reference interval)
Ma [39]	Pregnant women	Thyroid hormones	2	Indirect method, Altman method	Clinically relevant	UN	NR	NR
Markus [40]	Pregnant women	Hematology parameters	2	Indirect method, LMS model	Clinically relevant	UN	NR	Random split validation
Ma [11]	Apparently healthy adults (age ≧18 years)	22 biochemistry analytes	2	Indirect method, Altman method	Multifactor linear regression	Yes	NR	NR
Zhu [41]	Adults (age ≧20 years)	Thyroid hormones	2	Indirect method, LMS model	Z score, decision tree	Yes	NR	NR
Mokhtar [42]	Out-patients (age ≧18 years)	Thyroid stimulating hormone	2	Indirect method, quantile regression	NR	No	95% CI	NR
McCudden [43]	Patients (age ≧18 years)	Cerebrospinal fluid total protein	2	Indirect method, quantile regression	Significant correlation, visual identification	Yes	NR	NR
Zhang [44]	Apparently healthy adults (age ≧20 years)	Alpha-fetoprotein, carcinoembryonic antigen	2	Indirect method, Altman method	Correlation coefficient	Yes	NR	NR
Ammer [25]	Patients (age ≦18 years)	Alkaline phosphatase, creatinine, hemoglobin	3	Indirect method, GAMLSS(special)	Test new method	Yes	95% CI	Comparing results to continuous RIs from PEDREF and CALIPER, as well as to discrete RIs provided in manufacturers’ package inserts
Zierk [45]	In- and out-patients	15 biochemical analytes	3	Indirect method, kosmic method, smoothing spline function	Clinically relevant, biological variable	Yes	NR	NR
Haeckel [46]	In- and out-patients	High-sensitivity cardiac troponin T	3	Truncated minimum chi-square (R)	NR	Yes	NR	NR
Haeckel [47]	Hospitalized patients and outpatients	Uric acid, cystatin, creatinine	3	Truncated minimum chi-square (R)	NR	Yes	NR	NR
Wosniok [48]	Outpatients	Sodium, total protein, alanine aminotransferase, thrombocytes, γ-glutamyl transferase, high-sensitivity cardiac troponin T, hemoglobin	3	Indirect method, truncated minimum chi-square, cubic smoothing spline function	Biological variable	Yes	NR	NR
Zierk [10]	Outpatients	9 Hematology tests	3	Indirect method, Arzideh method, smoothing spline function	Biological variable	Yes	NR	NR
Weidhofer [49]	Patients (age ≦18 years)	Activated partial thromboplastin time, prothrombin time, thrombin clotting time, fibrinogen and antithrombin	3	Indirect method, Arzideh method, smoothing spline function	Clinically relevant, biological variable	Yes	NR	NR
Zierk [50]	Inpatients and outpatients (age ≤18 years)	Alkaline phosphatase	3	Indirect method, Arzideh method, smoothing spline function	Clinically relevant, biological variable	Yes	NR	Yes, comparison to results from CALIPER
Zierk [13]	Inpatients and outpatients (age ≤18 years)	22 Hematologic and biochemical analytes	3	Expanded version of algorithm developed by Arzideh, smoothing spline function	Clinical relevance and methodologic requirements	Yes	NR	Yes, comparison to results from KiGGS and CALIPER
Zierk [12]	Patients (from birth to adulthood)	Blood count quantities and sodium	3	Algorithm developed by Arzideh, smoothing spline function	Clinical relevance and methodologic requirements	Yes	NR	Yes, comparison to results from KiGGS and instrument specific RIs

CI: confidence interval; NR: no report; UN: unnecessary.

While continuous RIs serve as the precursor for next-generation reference interval models, they are not innovative per se. The novel aspect is the application of innovative algorithms specifically tailored to utilize medical data in the EHR, that is, real-world data. It is through these novel computational methods that we can accurately establish and refine continuous reference intervals. The subsequent sections provide detailed descriptions of the methodology and algorithms employed for establishing next generation RIs using next-generation reference interval models.

Direct sampling techniques with modeling (category I methods)

This method utilizes data primarily from multicenter studies originally designed to establish reference intervals using direct sampling. Selected reference subjects are employed to construct an age-test curve model. In these studies, reference subjects are selected using direct sampling, where individuals are sampled randomly from various centers. Based on when exclusion criteria and partitioning are applied, methods can be divided into a priori and a posteriori. Exclusion criteria encompass variables like blood pressure, body mass index, smoking, and alcohol consumption that are sourced from questionnaires and onsite physical examinations. CLSI document EP28-A3c, Establishing, and Verifying Reference Intervals in the Clinical Laboratory [51], provides a comprehensive questionnaire detailing these criteria. Upon obtaining reference individual test results, a curve model illustrating age-test result relationships is developed. Several methodologies exist for this purpose.

Altman method [52–54]

This method, based on Altman et al. starts with performing a Box-Cox power transformation or log transformation on the data. Then a weighted polynomial model is developed for the mean of the test (y₁) as a function of age, and it gives the standard deviation of the test (y₂) as a function of age. Finally, the models are used to establish reference intervals (y₁±Z × y₂) for every age within the age range; for RIs with a 95% confidence interval, the Z value is equal to 1.96. The calculated results are back-transformed to the original scale.

LMS (lambda-median-sigma) [55–57]

This well-established and widely-used technique, developed by Cole and Green, constructs age-related reference percentile curves using LMS [55]. In this method, with the three numerical parameters, L (degree of skewness), M (central tendency) and S (dispersion), continuous normal distributions are defined with age. As with the Altman method, this method requires Box-Cox transformation or log transformation on the data, and back-transformation to the original scale after calculating the age-specific distribution and the reference limits of the test.

GAMLSS model [58–60]

The full name of the GAMLSS model is General Additive Model for Location Scale and Shape method. Functions for GAMLSS were introduced by Rigby and Stasinopoulos in 2005. In this method, models are developed using distributional regression in which all parameters of the conditional distribution for the response variable are modeled using the explanatory variables. In addition to focusing on the mean (or location) of the distribution, other parts of the distribution like variance, quantiles and skewness and kurtosis or tails may also need to be considered to establish appropriate models. More than 100 continuous, discrete, and mixed distributions are available in GAMLSS for modeling the response variable, which makes modeling flexible and controllable.

Quantile regression [61–65]

Quantile regression is a non-parametric method. Unlike linear regression, which estimates parameters of the model by the least square method, the quantile regression estimates model parameters using the weighted least absolute deviation method to minimize the sum of the absolute residuals for each sample. Linear regression also relies on assumptions of normality and homogeneity of variances for residuals, which are often challenging to meet in practice, it is sensitive to outliers, and it fails to capture the global distribution of response variables in many scenarios. In contrast, quantile regression addresses these limitations effectively by not requiring residual distribution and by providing robust results that are less influenced by outliers.

The development of next-generation reference interval models through quantile regression often presents challenges:

1. flexibility: each fitted quantile curve should possess sufficient flexibility to capture the nonlinear patterns observed across ages;

2. comparison: the absence of an overall (likelihood-based) measure of fit, such as Generalized Akaike Information Criterion, in the fitted quantile regression model, poses challenges when comparing competing models;

3. non-crossing: ensuring non-crossing behavior is crucial for the fitted quantile curves, particularly when dealing with nonlinear relationships, limited sample sizes, and sparse data; and

4. instability of edges: the quantile curves near the extremes (i.e. for quantiles close to 0 or 100) vary more (i.e. are more erratic) than those in the center of the distribution of y, because those curves are supported by fewer observations.

Nonparametric estimation with radial smoothing

Kernel function is a function introduced in local regression to describe the local features of data. The radial basis function is a basis function in a space of functions. It satisfies the condition that for a fixed-point, c, the values of the function are the same for x equidistant around the fixed-point. There are many common radial functions, and the Gaussian function is one of them. The kernel method maps a data matrix to a new design matrix by using a kernel or a set of M basis functions. In this way, a complex nonlinear model can be fitted indirectly in the original input space. In the study of Wan et al. [66], the radial basis function was used to establish the mean, variance, skewness, and kurtosis models. Then, the Edgeworth-Cornish-Fisher expansion was applied to approximate quantile functions. Finally, rearrangement was used to transform the original quantile estimate to a monotonic quantile estimate.

Double-kernel-based method [67,68]

In 2010, Li et al. developed the double kernel-based method and automatic bandwidth selection procedure [67]. The method is nonparametric and based on the double-kernel percentile estimator proposed by Yu and Jones [68]; it can be obtained using the following steps:

1. Estimate the conditional distribution function of y given x using local-linear weighting and applying a kernel smoother with a specified bandwidth along the x-axis.

2. Simultaneously, apply another kernel smoother with its own bandwidth along the y-axis.

3. Finally, obtain the estimated percentile from the estimated conditional distribution.

The double-kernel estimator was favored over the single-kernel estimator by Yu and Jones because of its improved mean-squared-error properties and smoother appearance. In the method of Li et al., sample weight was incorporated into the double-kernel percentile estimator to fit the data.

Indirect sampling techniques (1) with modeling (category II methods)

These methods closely mirror category I methods, direct sampling techniques combined with a modeling method, in terms of their procedural approach. The primary distinction lies in its use of indirect sampling techniques to ascertain the underlying Gaussian distribution. As depicted in Table 2, data for this approach is extracted from individuals participating in routine physical exams during periodic health screenings as well as from blood bank donors, that is, individuals who are mostly healthy. Then, other data within the EHR are used to filter out subjects who may exhibit abnormalities [69]. To address outliers, methods such as Tukey’s method [70,71] and the latent abnormal value exclusion [72] are applied. Once the outliers have been identified and excluded, a reference population is determined. Subsequently, the next-generation reference interval model is constructed following the modeling methodologies delineated under category I methods.

Table 2. Comparison of three methodologies for establishing next-generation reference interval models.

	Category I. Direct sampling technique with modeling	Category II. Indirect sampling technique (1) with modeling	Category III. Indirect sampling technique (2) with curve fitting
Data Source	Multicenter study (usually)	Real-world data in the clinical laboratory	Real-world data in the clinical laboratory
Subjects	Apparently healthy individuals	Individuals undergoing routine physical examinations for periodic health screening in hospitals or from blood bank donors	Outpatient and inpatient data (usually), individuals undergoing routine physical examinations for periodic health screening in hospitals, blood bank donors
Reference individual acquisition	Direct sampling technique	Indirect sampling technique (1)	Indirect sampling technique (2)
Data used for curve fitting	All raw data	All original data	Reference limits obtained by indirect sampling techniques
Cost	High	Low	Low
Process	Cumbersome	Simple	Simple
Can it be used for new markers?	Yes	No, not unless the database has enough data	No, not unless the database has enough data

Indirect sampling techniques (2) with curve fitting (category III methods)

Over the last decade, there has been a growing interest in constructing next-generation reference interval models utilizing outpatient and inpatient data. The methodological flow of these investigations encompasses three principle steps that are depicted in Figure 3 [10,12]:

Figure 3. Procedure for developing a next-generation reference interval model based on indirect sampling (2) and curve fitting.

1. Segment the input data into overlapping age categories.

2. For each age group, establish discrete reference intervals using indirect sampling techniques.

3. Transition from discrete reference intervals to continuous percentile curves.

Category III methods are generally based on hospitalized patients and use data from mixed datasets that are extracted from hospital EHRs, which have a higher proportion of pathological results than the datasets used for category II methods, which are based on mostly health individuals. Initially, category III methods did not incorporate extensive data filtration. As research progressed, more rigorous steps were integrated to refine the data. In category III methods, advanced algorithms are used to extract the Gaussian distribution. This advanced approach is gaining traction because of the maturation of the tools used and the availability of large patient datasets. As indicated in Table 1, this methodology has been used extensively in the development of pediatric reference intervals.

Indirect sampling (2)

The common approaches for indirect sampling (2) used in the establishment of next-generation reference interval models consist of Hoffmann’s approach, the Bhattacharya approach, the expectation-maximization algorithm, the kosmic approach, the refineR algorithm, and the truncated minimum chi-square (TMC) algorithm. Among them, the TMC and kosmic algorithms are used most often. These methods have been summarized in our previous studies [73–75]. The algorithms are briefly described here.

Hoffmann

The Hoffmann method is an algorithm based on graphs [76,77]. First, the measurement values are arranged in ascending order and the cumulative probability for each value is calculated. The probabilities (y-axis) against the measurement values (x-axis) are plotted on normal probability paper. The y-axis scale also can be marked with the Z value of the standard normal distribution [77]. In the plot, the distinct linear zones correspond to different subgroups representing health and disease. The upper limit and lower limit of reference intervals are determined by extending the linear region of the healthy subgroup and calculating their corresponding x-axis values using 2.5 (-1.96) and 97.5 (1.96), respectively, in a specific formula, to obtain percentiles. The Hoffmann algorithm’s dependence on visual inspection may introduce subjectivity into the results of the analysis and impedes the potential for automated analysis. Although two studies [78,79] have enhanced Hoffmann’s method to enable automatic selection of linear regions, some scholars [77,80] contend that this modification has altered the fundamental essence of the Hoffmann approach and consequently has impacted the confidence level associated with the 95% reference interval. Thus the improved Hoffmann method remains controversial [81–83]. Consequently, incorporating this algorithm into Category III techniques of constructing next-generation reference interval models is challenging.

Bhattacharya

This method involves applying a logarithm transformation to the variables $\frac{\mathit{dln}\left(Y\left(x\right)\right)}{d x}$ and x, followed by calculating the derivative of the resulting normal distribution density function (Equation 1(1) $$\frac{\mathit{dln}\left(Y\left(x\right)\right)}{d x}=-\frac{1}{{\sigma }^2}\times x+\frac{\mu }{{\sigma }^2}$$(1) ). (1) $$\frac{\mathit{dln}\left(Y\left(x\right)\right)}{d x}=-\frac{1}{{\sigma }^2}\times x+\frac{\mu }{{\sigma }^2}$$(1)

To determine the mean and standard deviation of the healthy distribution, plot a scatter plot of x and $\frac{\mathit{dln}\left(Y\left(x\right)\right)}{d x}$ manually and fit an analytical formula to represent the linear region observed in the scatter plot. Then, utilize the slope and intercept values from this straight line to calculate these parameters. Subsequently, use these calculated mean and standard deviation values to compute the reference intervals [84]. It is important to note that both Hoffmann’s and Bhattacharya’s methods have been successfully applied for analytes with a normal distribution. However, caution should be exercised when dealing with analytes exhibiting skewed distributions as they may yield unreliable RIs. Additionally, as for the Hoffmann algorithm, the Bhattacharya algorithm relies on visual inspection that can introduce subjectivity into the analysis and that hinders automated analysis. Therefore, it is difficult to embed this algorithm in category III techniques for building next-generation reference interval models.

Expectation-maximization

The expectation-maximization algorithm is an iterative algorithm that consists of two steps: the expectation (E) step and the maximization (M) step) [85]. The E step calculates the model parameters and the probability for each individual. Using this probability, the M step then estimates the model parameters. After each iteration, the E step reevaluates the parameters until convergence con­ditions are met. Finally, at the end of the process, the RI is established based on the final parameters obtained from the analysis of the data distribution.

Compared to Hoffmann and Bhattacharya, the expectation-maximization algorithm allows a higher level of automation; however, it does necessitate parameter adjustment. If one does not actively engage in the parameter adjustment process, calculated reference intervals may become highly unreasonable. Consequently, incorporating this method into category III techniques for developing next-generation reference interval models poses a challenge.

kosmic

The proposed approach involves estimating the distribution of physiological test results by analyzing a mixture of normal and abnormal distributions [86]. Initially, the data is transformed using the Box-Cox transformation, followed by fitting a Gaussian distribution to a truncated section of the dataset. Subsequently, the Kolmogorov-Smirnov distance between this Gaussian distribution and the observed truncated distribution is calculated. The distribution with the smallest Kolmogorov-Smirnov distance is then chosen as representative data for individuals with good health. Finally, this estimated distribution is utilized to calculate the reference interval. The algorithm allows a high level of automation and is robust, making it one of the prominent methods in establishing next-generation reference interval models.

refineR

The refineR algorithm employs an inverse modeling approach and can be divided into three steps [87]. Firstly, the parameter search region and the primary peak are carefully selected. Secondly, a comprehensive multi-level grid search is performed to determine the optimal model parameters, which include λ (power), μ, σ, and P (scaling factor). Thirdly, the optimal model is used to calculate RIs. This algorithm is a relatively new method that exhibits superior stability and accuracy compared to the kosmic algorithm. Despite limited modeling studies utilizing this approach, its potential for practical application in establishing next-generation reference interval models is promising. It has performed well in several comparative studies [74,75,88]. Recently, Ammer et al. developed a new method for establishing next-generation reference intervals model based on refineR and GAMLSS [25]. This novel method can automatically and objectively optimize the smoothing strength. Furthermore, in comparison to other methodologies of category III techniques using nonparametric method fitting curves, this approach employs distribution parameters that are estimated through an indirect technique. This facilitates a weighted contribution from all data points and effectively controls the influence of patients results on the estimation of smooth curves. By utilizing the GAMLSS approach, the semi-parametric model developed enables computation of any percentile for any given explanatory variable. Because this is a recent model and research using it is limited, this methodology is classified as a special approach in category III methods of establishing next-generation reference intervals.

Truncated minimum chi-square (TMC)

The TMC algorithm [48] is well-suited for analyzing data distributions that exhibit high skewness and contain a substantial proportion of data at or below the detection limit. The approach begins with utilizing a quantile-quantile plot to obtain initial estimates for the parameters (λ, μ, σ) of the assumed power normal distribution (Note that while the parameters λ, μ and σ here in the TMC method serve a similar role to the LMS parameters in that they describe a distribution’s skewness, median, and variation, they are derived and applied differently in the two methods.). These estimates are then refined iteratively using the TMC estimation. The resulting parameters are used to calculate a 95% reference interval. Confidence intervals for the reference limits are determined using an asymptotic formula for quantiles, while tolerance limits are obtained through bootstrapping techniques. The TMC method [46,47] is employed widely in establishing next-generation reference interval models because of its inherent stability and ease of automation.

Fitting percentile curves

Curve fitting, the third step of category III techniques, involves determining the correlation between the reference limits and age. These approaches can be a parametric method or a smoother. Smoothers use nonparametric approaches to determine the relationship between explanatory variables (x) and the predictor (y). Here, the commonly used approaches are described

Polynomial regression

This is a linear regression model employed to capture nonlinear relationships, thereby overcoming the limitations of simple linear regression (which accommodates solely linear associations between explanatory variables (x) and the predictor (y)). The polynomial regression method incorporates additional polynomial features, where $\beta$ denotes the parameter that requires training of the model. The model parameters of polynomial regression can be estimated through either the method of least squares or gradient descent (Equation 2(2) $$f\left(x\right)={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+{\beta }_3{x}^3+\dots +{\beta }_p{x}^p$$(2) ). (2) $$f\left(x\right)={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+{\beta }_3{x}^3+\dots +{\beta }_p{x}^p$$(2)

In particular, when the power p of $x^p$ in a polynomial regression is negative, $x^p$ is a fractional term. This approach, pioneered by Royston and Altman in 1994 [89], aimed to achieve more accurate curve fitting with the minimum number of polynomials.

However, the use of polynomial regression has several challenges. Firstly, as the number of terms increases, while the data may be increasingly well-fitted, overfitting becomes a concern. In other words, while the current data may exhibit a perfect fit, there is a high likelihood of poor performance when it is applied to new or modified data. Secondly, collinearity among the terms is an issue because correlations between different variables can lead to difficulties in accurately estimating the individual effect of each variable. Lastly, polynomial regression employs a single function to model global data; however, in many cases, distinct patterns within different subsets of the dataset emerge that necessitate the use of diverse functions for accurate fitting.

Regression spline

Data may exhibit sudden changes at specific values. To address such a data modeling scenario, piecewise polynomial regression can be employed for processing. The x value corresponding to each segment is referred to as a knot or breakpoint. Piecewise polynomial regression, also known as regression spline according to Smith [90], necessitates the continuity of all derivatives at these nodes. Commonly utilized splines include B-splines and natural splines. Splines are defined as Equation 3(3) $$f\left(x\right)=\sum _{j=0}^D{\beta }_{0j}{x}^j+\sum _{i=1}^K{{\beta }_i\left(x-b_i\right)}^{D}H(x-b_i)$$(3) , where D is the degree of the polynomial in x and K is the number of knots. When D = 3, the cubic spline is obtained. Cubic splines are particularly useful for smoothing techniques because of their continuous first and second derivatives at the breakpoints, resulting in exceptionally smooth curves. B-splines [91,92], or basic splines, are a set of smooth, piecewise polynomials that are defined over a certain range and are used extensively in numerical analysis and data fitting. They are particularly useful because of their stability and minimal support with respect to a given partition of the domain. B-splines of a given degree are continuous and have continuous derivatives up to a certain order. Each B-spline basis function is non-zero only over a small interval of the domain, which makes them computationally attractive for modeling complex shapes and functions. (3) $$f\left(x\right)=\sum _{j=0}^D{\beta }_{0j}{x}^j+\sum _{i=1}^K{{\beta }_i\left(x-b_i\right)}^{D}H(x-b_i)$$(3)

$H\left(x-b_i\right)$ is the Heaviside function: $$H\left(t\right)=\left\{\begin{array}{c}0,t<0\\ 1,t\ge 0\end{array}\right.$$

The use of B-splines for truncated piecewise polynomials is analogous to that of orthogonal polynomials for polynomial regression. The number of knots has a direct impact on the model’s degrees of freedom and subsequently affects its complexity. A higher number of knots leads to increased model complexity, larger calculations, and reduced curve smoothness. Therefore, careful consideration should be given to selecting the appropriate number of knots, typically ranging from 3 to 5. One limitation with B-splines is that they tend to exhibit high variance in prediction at both ends of the x value range, resulting in wide confidence intervals. To address this issue, natural splines that impose boundary restrictions on regression splines are used; the assumption of linearity in the leftmost and rightmost segments enhances accuracy in predicting boundaries.

Smoothing spline

For the development of next-generation models, the main considerations include model goodness of fit and model generalization; this means that the fit is good and that the model represents the original data as much as possible. In the best case, the model can accurately predict every data point. However, the biggest problem is that the prediction ability of the model for a new data set may be poor, that is, that the model generalization ability will be poor. This situation still exists in the development of the next-generation reference interval models. Thus, smooth splines [93] have been introduced (Equation 4(4) $$\sum _{i=1}^n{\left[y_i-f\left(x\right)\right]}^2+\lambda \int {\left[f^{″}\right(x\left)\right]}^{2}dx$$(4) ). (4) $$\sum _{i=1}^n{\left[y_i-f\left(x\right)\right]}^2+\lambda \int {\left[f^{″}\right(x\left)\right]}^{2}dx$$(4)

The smaller the first term$,\hspace{0.17em}{{\displaystyle \sum }}_{i=1}^n{\left[y_i-f\left(x\right)\right]}^2$, in Equation 4(4) $$\sum _{i=1}^n{\left[y_i-f\left(x\right)\right]}^2+\lambda \int {\left[f^{″}\right(x\left)\right]}^{2}dx$$(4) , the lower the fitting error; however, it also increases the risk of overfitting. The second term, $\lambda \int {\left[f^{″}\right(x\left)\right]}^{2}dx$, quantifies the function’s unsmoothness, while $\lambda$ serves as a penalty term for model complexity. The penalty becomes more severe and the curve becomes smoother as $\lambda$ increases. The $\lambda$ is typically determined through cross-validation. The fundamental principle of smooth splines involves integrating a collection of B-spline basis functions with a smoothing penalty to control the smoothness of the resulting curve. B-splines in this framework typically have a degree of three, hence cubic splines are frequently used in smooth splines. When the degree is set to three, smooth splines are characterized as piecewise cubic polynomials, and linear functions are often utilized for the outermost segments to ensure stability at the data’s boundaries.

Local regression

This non-parametric curve fitting method involves selection of a small local neighborhood, also known as bandwidth, around each independent variable. Points within this neighborhood are used to fit a local zeroth, linear and quadratic polynomial. The polynomial in the local neighborhood, a zeroth polynomial, is referred to as kernel regression [94]. This curve-fitting approach is called local polynomial regression when the local fit involves linear and quadratic polynomials. Specifically, in local regression, the number and position of fitting points are determined first. Subsequently, the m nearest points are determined with x_i as the central point for fitting, and the weight of these m points is computed using weight functions. Finally, linear or quadratic polynomial fitting is performed through weighted regression in the local neighborhood and the value of the regression function at x_i is obtained as the dependent variable. The final curve can be drawn after repeating the above steps for all the original data points. In kernel regression, after obtaining weights for x_i, y_i can be calculated directly using Equation 5(5) $$y_i=w_1{x}_1+w_2{x}_2+\dots +w_m{x}_m$$(5) . (5) $$y_i=w_1{x}_1+w_2{x}_2+\dots +w_m{x}_m$$(5) w_m is the weight of m raw data in the neighborhood centered on x_i, where the closer w_m is to x_i, the greater the weight. The selection of bandwidth plays a crucial role in model performance. A wider bandwidth results in a smoother curve, potentially leading to underfitting. Conversely, a narrower bandwidth can lead to a more intricate curve and increase the risk of overfitting.

Comparison of the three methods

The three methods for constructing next-generation reference interval models differ primarily in their choice of reference individuals and curve-fitting approaches.

The primary distinction between category II and category III techniques is their approach to data selection and refinement. In category I methods, reference subjects are selected using direct sampling, Category II methods rely primarily on datasets with a comparatively low proportion of pathological results (individuals undergoing routine physical examinations or blood donors), and these methods focus on outlier identification. Category 3 methods use mixed datasets that have a higher proportion of pathological results and use algorithms to obtain a Gaussian distribution. Further distinctions are detailed in Table 2. Compared to the direct sampling technique, the indirect sampling method presents clear advantages in cost-effectiveness and operational feasibility in model development. Notably, category III techniques, which combine the indirect sampling technique (2) with percentile curve fitting, depend on reference limits acquired through indirect sampling. In contrast, the first two methods use original data.

It should be emphasized that the pipeline developed by Ammer et al. [25] is a special method based on the use of refineR and GAMLSS. Ammer’s method has gained widespread acceptance and is now considered a standard method. It utilizes the refineR algorithm to separate non-pathological distributions from mixed data. The results of individuals belonging to non-pathological distributions are used to establish the next-generation reference interval model using GAMLSS. Unlike the other algorithms used in category III techniques, this method does not rely on reference limits obtained from indirect sampling technology to fit curves. However, this method also does not use original data. In summary, considering that this class of methods is relatively new and that an indirect sampling technique (2) is applied, this method is classified as a special subclass of category III methods.

Diagnosis, validation and testing of next-generation reference interval models

The process of model diagnosis and validation is paramount in ensuring the precision and reliability of the finalized model. Diagnostic evaluations are indispensable for certain parametric, semi-parametric, and even non-parametric models. For all models, validation is essential. This segment offers an overview of the methods for diagnosis and validation.

Normalized quantile residuals for diagnostics

The raw residuals as well as the Pearson residuals pose challenges in generalizing to distributions beyond normality. Additionally, some response variables in the models exhibit extreme skewness. To address these issues and enhance parametric model diagnostics, normalized quantile residuals are introduced [95,96]. The process for calculating normalized quantile residuals involves determining the probability density function and cumulative distribution function of the observed value y, followed by computing the corresponding function F(y) on the cumulative distribution function. Finally, the corresponding x-value of F(y) on the standard normal distribution’s cumulative distribution function is calculated to obtain the residual. Fit residuals are computed using a similar approach. In a correct model, both residuals and fit residuals should conform to a standard normal distribution. For visual inspection of normality, the GAMLSS package in R provides several plots that include residual-related plots (Figure S1). Adherence of these plots to a standard normal distribution for residuals demonstrates satisfactory model performance. Worm plots [97,98] (Figure S2), used primarily for assessing appropriate fitted distributions, show good fit when points lie between two curves and align closely with a central horizontal line. The Q-statistic plot (Figure S3) is used to detect deviations from normality, with respect to the mean (z₁), variance (z₂), skewness (z₃), and kurtosis (z₄) values of residuals. The manifestation of squares on this plot is indicative of absolute Z-values surpassing 2, which signals potential misalignment of certain distributional parameters within the model [99].

k-fold cross validation

The k-fold cross-validation, one of multiple methods for internal validation, is used primarily for comparing models competitively. This method initially partitions the dataset into k subsets, with one subset designated as the validation set and the remaining k-1 subsets used for training. The model’s goodness of fit is then computed for each validation set, and after conducting cross-repeated validation k times, the average result from these iterations is employed as the final evaluation metric for model validation [31,32].

Calculating the fraction outside the reference interval

The methodology for determining the fraction outside the reference interval (FOR) [34,38] is illustrated in Figure S4. Data from healthy individuals, distinct from the modeling dataset, were curated to form the test set. From this collection, a random 30% was designated as the test subset. Utilizing the pre-established model, the test results for individuals within this subset were compared against the reference intervals corresponding to sex and age. Subsequently, the percentage of these individuals who fell outside the designated reference set was computed. This evaluation process was iterated 100 times, yielding multiple FOR values. From the distribution of these FOR values, a 95% confidence interval was derived. In the next-generation reference interval model with a two-sided 95% confidence level, the expected percentage of observations falling outside the reference interval stands at 5%. This is split evenly, with 2.5% anticipated beyond each of the upper and the lower bounds. Theoretically, the 95% confidence interval of the FOR value should encompass this 2.5% threshold. Should it not align with this expectation, it is essential to adjust the input parameters in the modeling approach to account for the observed discrepancies. The modeling is iteratively refined until the desired specifications are satisfied.

Tools for establishing big data driven next-generation reference interval models

The tools and code sources commonly employed for establishing next-generation reference interval models are summarized in Table 3. Currently, the predominant modeling approach relies on programming, with R [100] and Python [101] being the primary languages utilized because they have strong functionality and easy-to-use packages and libraries. Additionally, RStudio [102] and PyCharm [103] provide convenient platforms for programming in R and Python. Also, Medcalc or SPSS can be used.

Table 3. Tools for establishment of next-generation reference interval model.

Methods	Example study (title)	References	Uniform resource locator of code or software
Altman method	Age-specific pediatric reference intervals for plasma free normetanephrine, metanephrine, 3-methoxytyramine and 3-O-methyldopa: Particular importance for early infancy	[36]	Software: Medcalc
LMS	LMS-based continuous reference percentiles for 14 laboratory parameters in the CALIPER cohort of healthy children and adolescents	[26]	R package:gamlss
GAMLSS	Establishment of reference interval and aging model of homocysteine using real-world data	[38]	R package:gamlss
Double-kernel-based method	Application of nonparametric quantile regression to body mass index percentile curves from survey data	[67]	R package:Stats
Radial smoothing	Nonparametric estimation of age-specific reference percentile curves with radial smoothing	[66]	Software:SAS (The code is at the end of the study)
Quantile regression	Continuous reference intervals for pediatric testosterone, sex hormone binding globulin and free testosterone using quantile regression	[32]	1. R code (The code is at the end of the study) 2. SPSS
TMC and splines	Age- and sex-dependent reference intervals for uric acid estimated by the truncated minimum chi-square (TMC) approach, a new indirect method	[47]	TMC R Code:https://www.allacronyms.com/TMC/Truncated_Minimum_Chi-square
kosmic and splines	Next-generation reference intervals for pediatric hematology	[10]	https://kosmic.diz.uk-erlangen.de/
refineR and GAMLSS	A pipeline for the fully automated estimation of continuous reference intervals using real‑world data	[25]	https://www.nature.com/articles/s41598-023-40561-3

Future perspectives of the next-generation reference interval models based on data from clinical laboratories

Novel report pattern of the next-generation reference interval model

Next-generation reference interval models can be represented using quantile plots or normalized Z-score plots, which offer the advantage of eliminating the influence of units when comparing test results across different platforms. In future reports, a single RI may be presented for the “non-elderly adult” population or simply adult population—those under 65 or 60 years of age, as depicted in Figure 4. This interval will offer a baseline against which the health metrics of “older adults”, those over 60 or 65 years, can be normalized. By comparing the laboratory values of older adults with this RI, we can delineate the extent of their health status deviations attributable to the aging process. This comparison is designed to enrich clinical decision-making by providing a clear benchmark for identifying age-related changes in health parameters. Furthermore, while sex and age are often used in model construction because they are readily available in patient records, other determinants also play a role. The emphasis on age and sex in our review arises from their universal accessibility, e.g. in the EMR, making models based on them particularly practical in everyday clinical work. Genetic background and other variables are also important, but the challenge lies in the lack of availability of such data for all patients. This may pose limitations in implementing models in real-world scenarios. Nevertheless, it is our belief that acknowledging the many influences on test results is paramount. In the future, to cater to the intricacies of multi-variable scenarios, one approach might be to construct segmented models. Ideally, integrating these variables for model adjustment would be the optimal route and would significantly enhance the precision and dependability of the model.

Figure 4. Novel report pattern of the next-generation reference interval.

Develop an information system that can use next-generation reference interval models

Our team has developed a multi-center big data analysis platform for real-world studies of clinical laboratory data that can be utilized in training, validating, and testing next-generation reference interval models. Additionally, we aim to establish data interfaces between our platform and both the clinical laboratory information system and the hospital information system to transfer quantile values derived from the model. In the future, test reports will include plots of patient results in next-generation reference interval model diagrams, providing clinicians with a visual basis for decision-making. Importantly, upon opening the interface, data will be analyzed and transmitted in real-time, allowing regular re-estimation of next-generation reference interval model parameters using the new data. The dynamic nature of next-generation reference intervals based on next-generation reference interval models will enable more informed evaluation of patient results.

Authors’ contributions

Chaochao Ma wrote and revised this manuscript. Zheng Yu provided guidance on the mathematics. Ling Qiu made suggestions for the revision of the manuscript. All authors reviewed the manuscript and approved the submission.

Consent for publication

All the authors gave their consent for publication.

Abbreviations
EHR	=	electronic health record
FOR	=	fraction outside the reference interval
GAMLSS	=	generalized additive models for location, scale, and shape
RI	=	reference interval
TMC	=	truncated minimum chi-square

Disclosure statement

This review is part of the Real-World Study of Big Data Mining in Laboratory Medicine and has been conducted in accordance with the Declaration of Helsinki. This review was approved by the Ethics Committee of Peking Union Medical College & Chinese Academy of Medical Sciences, Peking Union Medical College Hospital (approval number: S-K1192). All methods were performed in accordance with relevant guidelines and regulations.

Data availability statement

The authors declare that all data generated or analyzed during this study are included in this article.

Additional information

Funding

This study was supported by the National Natural Science Foundation of China (72274218).