Age- and Sex-Adjusted Reference Intervals in Tear Cytokine Levels in Healthy Subjects

Alterations in tear cytokine levels have been associated with various ocular disorders as compared to those in healthy subjects. However, age and sex are not always considered in these comparisons. In this study we aimed to establish age and sex reference intervals (RIs) for tear cytokine levels in healthy people. Tear samples were taken from 75 males and 82 females, aged 18–88 years, and tear cytokine levels were determined. Age- and sex-adjusted RIs for epidermal growth factor (EGF), fractalkine, interleukin (IL)-1 receptor antagonist (RA), IL-7, IL-8, interferon inducible protein (IP)-10, monocyte chemotactic protein (MCP)-1, and vascular endothelial growth factor (VEGF) tear cytokine levels in a healthy sample were established using generalized additive for location, scale and shape (GAMLSS) models. RIs were tested in two external samples: a validation sample of 40 individuals with normal results at four Dry Eye Disease (DED) clinical diagnostic tests (OSDI, T-BUT, corneal staining and Schirmer test); and a utility sample of 13 severe DED cases. IL-1RA, IL-8, IP-10, and MCP-1 levels showed a positive association with age, while EGF was negatively correlated. IL-7 concentration increased up to 40 years and again after 70 years, observing a quasi-linear decrease between them. For VEGF, higher levels were observed in the middle-aged range. Regarding sex-influence, fractalkine tear levels were higher in men, whereas those of IL-7, IL-8, and IP-10 were higher in women. Using the estimated age- and sex-adjusted RIs, more than 92% of the validation sample was correctly classified, and 100% of the severe DED patients in the utility sample had concentrations outside the RIs in at least two of the cytokines evaluated.


Introduction
Cytokine profiles have been proposed as biomarkers for numerous infectious and chronic diseases due to their pro-and anti-inflammatory effects [1]. Numerous studies have also reported important alterations in tear cytokines associated with various ocular disorders, such as primary open-angle glaucoma, ocular chronic graft versus host disease, uveitis, ocular allergies, and dry eye disease (DED), among others [2][3][4]. In particular, a kurtosis, respectively. Each of the parameters can be a parametric or a non-parametric smooth function of the explanatory variables. In this way, GAMLSS allow to overcome some of the generalized linear models' limitations. First, the distributional assumption for the response variable is relaxed and a general distribution family can be used. In addition, we can model as a function of explanatory variables, not only the location parameter but also other parameters of the distribution.
Different distributions were fitted to the observed distribution of each molecule concentration on log2 scale. Four distribution families were considered: normal; Box-Cox Cole and Green to explore the existence of skewness in the data; power exponential to explore the possibility of kurtosis; and Box-Cox power exponential distribution, to show if both skewness and kurtosis are present in the data. The age influence on parameters of the considered distributions was modeled as a constant, as a linear function or as a cubic spline. The sex was added to the model as qualitative variable. For a particular cytokine, Akaike Information Criterion (AIC) [24] was used in a stepwise algorithm to choose the best model. In the framework of each distribution family, the strategy for selecting relevant explanatory variables for each of the parameters that characterize the distribution is summarized in Figure 1. Finally, the best model was chosen by selecting the distribution whose final model was better fitted to the data.
(GAMLSS) [17]. GAMLSS are very general and flexible methods for modeling a response variable, in our case a cytokine level, as a function of explanatory variables: age and sex. Any parametric distribution, including highly skewed and kurtotic distributions, can be fitted to the response variable by modeling four parameters: a location parameter, denoted by µ; a scale parameter (σ); and two shape parameters: ν and τ modeling skewness and kurtosis, respectively. Each of the parameters can be a parametric or a non-parametric smooth function of the explanatory variables. In this way, GAMLSS allow to overcome some of the generalized linear models' limitations. First, the distributional assumption for the response variable is relaxed and a general distribution family can be used. In addition, we can model as a function of explanatory variables, not only the location parameter but also other parameters of the distribution.
Different distributions were fitted to the observed distribution of each molecule concentration on log2 scale. Four distribution families were considered: normal; Box-Cox Cole and Green to explore the existence of skewness in the data; power exponential to explore the possibility of kurtosis; and Box-Cox power exponential distribution, to show if both skewness and kurtosis are present in the data. The age influence on parameters of the considered distributions was modeled as a constant, as a linear function or as a cubic spline. The sex was added to the model as qualitative variable. For a particular cytokine, Akaike Information Criterion (AIC) [24] was used in a stepwise algorithm to choose the best model. In the framework of each distribution family, the strategy for selecting relevant explanatory variables for each of the parameters that characterize the distribution is summarized in Figure 1. Finally, the best model was chosen by selecting the distribution whose final model was better fitted to the data. Stepwise algorithm for selecting relevant explanatory variables given a distribution family. Forward selection stage starts with the simplest model, the model with all the parameters fitted as constants, and tests the addition of sex and age, as a constant, a linear function or a cubic spline, using the Akaike Information Criterion (AIC). Backward elimination stage starts with the best model fitted in the forward stage, and tests the deletion of each explanatory variable using the AIC. The distributional parameters are: µ location; σ scale; ν skewness; and τ kurtosis. Stepwise algorithm for selecting relevant explanatory variables given a distribution family. Forward selection stage starts with the simplest model, the model with all the parameters fitted as constants, and tests the addition of sex and age, as a constant, a linear function or a cubic spline, using the Akaike Information Criterion (AIC). Backward elimination stage starts with the best model fitted in the forward stage, and tests the deletion of each explanatory variable using the AIC. The distributional parameters are: µ location; σ scale; ν skewness; and τ kurtosis.

Model Diagnostics
To check the adequacy of GAMLSS fitted models we used normalized (randomized) quantile residuals [25]. When the model is correct, residuals will have a standard normal distribution, with mean 0, variance 1, and skewness and kurtosis 0 and 3, respectively. Each fitted model was checked by using summary statistics of the quantile residuals and four graphical tools: residuals against the fitted values of the µ parameter; residuals against each one of the included explanatory variables; a kernel density estimate of the residuals; and the QQ-normal plot of the residuals. In addition, worm plots [26] were used as a diagnostic tool to assess if adjustments for some of fitted model parameters were required (see Residual Analysis in Appendix A for complete information).

Reference Intervals Utility
In order to know the utility of the age-and sex-related RIs estimated with GAMLSS model, we used them to classify subjects from two external independent samples, hereinafter referred to as validation and utility samples, obtained from our own database. Individuals with tear cytokine concentration data along with ocular surface clinical parameters related to the DED diagnosis were selected.
The validation sample consisted of 40 healthy individuals with a mean age of 31.74 ± 15.85 years. There were 19 (47.5%) males and 21 (52.5%) females. Their respective age was 39.31 ± 17.92 years for males and 24.90 ± 9.88 years for females. This healthy sample had normal results at four DED diagnostic tests: the ocular surface index (OSDI) score ≤ 12 points, fluorescein tear break-up time (T-BUT) ≥ 7 s, corneal fluorescein staining (Oxford scale) = 0, and Schirmer test without topical anesthesia ≥ 5 mm in 5 min [27]. The subjects of this sample should be classified as "normal" subjects by the RIs.
By contrast, the utility sample was used to evaluate the capability of these RIs for the diagnosis of DED. Consequently, tear molecule data from 13 cases of severe DED were selected. There were 4 (30.8%) males and 9 (69.2%) females. Their age was 58 ± 7.39 years for males and 58.33 ± 9.72 years for females. The clinical selection criteria for these subjects were having an OSDI score ≥ 33 points, T-BUT < 3 s, corneal fluorescein staining ≥ 3, and Schirmer test < 3 mm in 5 min [27]. Tear sample collection and tear cytokine concentration analysis in those subjects had been performed following the same methodology as in this study.
It should be noted that information about the concentration of some cytokines was missing in some individuals' tear samples. Specifically, levels of IL-7 were missing in all the subjects.
The percentage of cytokine tear concentration values falling within the lower and upper RIs of each sample was determined.
Individual tear samples were analyzed. The minimum detectable concentrations (in pg/mL) and detection rates for the molecules analyzed are shown in Table 1.
Eight (EGF, Fractalkine, IL-1RA, IL-7, IL-8, IP-10, MCP-1, and VEGF) out of the 18 molecules analyzed were detected in at least the 82% of the participants. Table 2 summarizes the distribution of the levels of these eight cytokines by age and sex.
The rest of the molecules analyzed (IL-1β, IL-2, IL-6, IL-10, IL-12p70, IL-15, IL-17A,  IL-23, TNF-α, and IFN-γ) were not further considered for subsequent statistical analysis because their detection rates were below 50% (Table 1).  Table 3 shows the best GAMLSS models of each cytokine. Model selection involved the selection of the distribution for the corresponding cytokine and the selection of relevant explicative factors by the comparison of fitted models according to the smallest AIC. Four different distribution families were considered: Normal, Box-Cox Cole and Green, Power Exponential, and Box-Cox power exponential. For IL-8, IP-10, and VEGF, Power Exponential was chosen as the most appropriate distribution. Normal distribution was selected for Fractalkine, IL-7, and MCP-1. Finally, Box-Cox power exponential was the less selected distribution, chosen in two of the eight cytokines considered (EGF and IL-1RA).
The potential explanatory variables were age and sex. The age influence on the location parameter was common to EGF, IL-1RA, IL-7, IL-8, IP-10, MCP-1, and VEGF, but not for Fractalkine, whose mean parameter depends exclusively on sex. In addition, age was added to model variability of IP-10 and skewness of VEGF. Sex had a relevant influence on location parameter of Fractalkine, IL-7, IL-8, and IP-10 and on variability of IP-10 (Table 3). A residual analysis of each of the final fitted models was carried out. In this diagnostic stage, none of the cytokines showed inadequacies in the final fitted model (see Residual Analysis in Appendix A).
The age-and sex-related 90% reference (5%, 95%) and median (50%) curves calculated using the best-fitted GAMLSS model for each cytokine are shown in Figure 2; complete numerical values can be found in Appendix B. The levels of IL-1RA, IL-8, IP-10, and MCP-1 showed an age-related increase; in the case of IP-10 and MCP-1 this increase was approximately linear but for IL-1RA and IL-8, the increase was only relevant in older subjects (70 years old or more). In contrast, EGF tear concentration decreased linearly with age.
Age patterns were also observed for IL-7 and VEGF tear levels. IL-7 concentration increased up to approximately 40 years, then a mild decrease was observed up to 70 years and then increased again. For VEGF, higher levels were observed in the middle age range (30 to 60 years approximately). Regarding sex influence, Fractalkine levels were generally higher in men, whereas IL-7, IL-8, and IP-10 tear concentrations were higher in women ( Figure 2).
In order to validate the estimated RIs, we used them to evaluate the tear cytokine concentrations from an external independent sample from our laboratory database of 40 healthy individuals (validation sample). In addition, tear cytokine concentration values from 13 cases of severe DED (utility sample; also from our own database) were also used to show the utility of the estimated RIs for diagnosis of DED. RIs were evaluated for all cytokines except for IL-7 whose levels were missing in all subjects of both samples. Table 4 summarizes the proportion of values that fell within the appropriate age-and sex-specific RIs for all other cytokines in each sample. The complete demographic information and classification for each subject in the validation and utility samples is shown in Tables A11 and A12 (respectively) in Appendix C. More than 92% of validation sample values fell within the appropriate age-and sex-specific RIs for each cytokine. Concerning the utility sample, IL-1RA and IL-8 levels were above the upper RIs in 92.3% of those patients diagnosed with severe DED compared to 7.7% and 5% of those clinically classified as not having DED. Fractalkine and VEGF tear levels were also increased in the DED sample; in these cases, the percentages of values over the upper RIs were 72.7% and 75% versus 5.6% and 6.3% of the validation sample, respectively. On the other hand, EGF and IP-10 showed lower concentrations in the utility sample, finding that 23.1% and 38.5% of the DED patients had levels below the corresponding lower RI. These percentages are significantly higher than those found in the validation sample (0% and 2.9%, respectively). Interestingly, all the patients of severe DED sample had tear cytokine concentrations that fell outside the RIs in at least two of the cytokines evaluated (see Table A12 in Appendix C). and 6.3% of the validation sample, respectively. On the other hand, EGF and IP-10 showed lower concentrations in the utility sample, finding that 23.1% and 38.5% of the DED patients had levels below the corresponding lower RI. These percentages are significantly higher than those found in the validation sample (0% and 2.9%, respectively). Interestingly, all the patients of severe DED sample had tear cytokine concentrations that fell outside the RIs in at least two of the cytokines evaluated (see Table A12 in Appendix C).    Two external samples are evaluated: (1) validation sample, consisting in the cytokine levels of 40 healthy subjects, and (2) utility sample that involved 13 severe cases of dry eye disease (DED). The percentage of values falling both within and outside the lower and upper reference limits is showed. CI = Confidence interval. * For MCP-1, there was only a valid subject in utility sample and he was classified as Low. IL-7 RI was not validated as its concentration information was missing in all subjects of both samples.

Discussion
In this work, we have established the sex-and age-adjusted RIs of eight (EGF, Fractalkine, IL-1RA, IL-7, IL-8, IP-10, MCP-1 and VEGF) tear cytokine levels. RIs are used to describe the normal range of a parameter, and their important role in clinical medicine is beyond discussion as an initial indication for further investigation or treatment. They provide a basis for comparison of measured values of patients to a reference, and they could be very useful to diagnose the condition of a subject in epidemiological and clinical studies.
The most commonly used definition of a RI is the interval of values containing the central 90% or 95% of a healthy population [28]. It is possible to estimate population percentiles from an adequate sample by computing observed percentiles. As a result, single cut-off values will be obtained. However, this result may not be very useful in the case of cytokine levels, since the distribution of cytokine concentrations, as inflammatory marker measurements, could be heavily dependent on individual characteristics, such as age and sex among others [29]. In addition, immunological data very frequently have asymmetrical distributions [30] and statistical techniques that take into account this issue should be chosen. For these reasons, in this study, GAMLSS models [17] were applied to derive ageand sex-dependent centile curves. The main advantage of these type of models is that the shape of the distribution, rather than just the mean or the variance, for the response variable can vary according to more than one explanatory variable. In our case, for example, the skewness parameter for VEGF concentration varied according to the subject age and it would be wrongly modeled by other more classical statistical methodology.
Our results showed that tear levels of all studied cytokines, except Fractalkine, presented an age-dependent behavior in any of the distributional parameters. Age-associated changes in human tear film composition and proteome have already been shown in other studies [12,15,16] in agreement with that reported as well in other fluids, such as plasma or cerebrospinal fluid [31]. Indeed, the relationship between the aging process and immune response is not new. Even, the expression "inflamm-aging" has been coined to describe the progressive increase in pro-inflammatory status as age increases [32,33]. There is a lot of literature about age, the role of inflammation markers and its link with the development of disease (reviewed in [31]). The pro-and anti-inflammatory cytokine balance changes with increasing age, and low-grade inflammation is more common among elderly population, leading to a chronic low-grade inflammatory phenotype [34]. Accordingly, age-dependent behavior of several cytokines levels has widely been reported in different biological fluids, such as plasma, cerebrospinal fluid, and also in tears, although tear studies are still scarce. Nevertheless, cytokine levels in these different biological sources have shown to have similar evolution patterns and common to our results in tears. For instance, in our study MCP-1 and IP-10 tear levels showed a sustained age-dependent increase, similar to what has already been observed in serum levels [35]. In addition, IP-10 increased linearly with age in cerebrospinal fluid [36], and plasma MCP-1 levels were significantly higher in an elderly group (70 to 92 years) compared to adults (21 to 50 years) [37]. This latest study also showed that plasma levels of EGF were significantly lower in the elderly group, which is compatible with our results in tears. In our study, IL-1RA and IL-8 tear levels showed as well a very important increase but only in old ages (68 to 90 years). This same pattern has already been reported in tear [12] and in serum levels [38] of IL-8, and in plasma levels of IL-1RA [39]. Even, these two cytokines have been proposed as markers linked to longevity in very old people [38,40]. An increase of IL-8 levels with age has recently been described in cerebrospinal fluid as well, but, in this case, the evolution is practically linear [36].
The role of IL-7 in immunosenescence, the aging of the immune system, has been known for long [41]. In this work, we have found a cyclic evolution of IL-7 tear levels with higher levels around 30 years old and increasing again after 70 years old. This result disagrees with the findings by Micera et al. [12] which found an increase of IL-7 in middle age people (41-60 years old) with respect to young  and old (<60 years old). However, the observed behavior for IL-7 tear levels in our study could be consistent with the findings of Nasi et al. [42], that evaluated plasma IL-7 levels in three age groups: young (aged 20-45), middle-aged (aged [58][59][60][61][62], and centenarians (range 98-100 years). Despite the fact that they found no significant differences between the age groups, the evolution pattern of IL-7 levels was similar to ours. There is more confusion regarding VEGF evolution. In serum, for example, both an increase [43] and a decrease [44] of VEGF levels with aging have been found. In plasma samples of patients from 43 to 80 years old with stable coronary artery disease a significant negative correlation was observed [45]; in contrast, recently, Chakraborty et al. [46] detected a positive correlation in cerebrospinal fluid from subjects aged from 40 years. Other authors such as Rübenhagen et al. [47] found no association between the concentration of VEGF in synovial fluid and age. In our case, the relation between VEGF tear level and age can be described as an inverted-U shape curve, with maximum levels between 30-60 years old.
Our results regarding the age-dependence in tear cytokine levels are also in agreement with the findings by Nättinen et al. [16] for tear proteome in which the majority of the identified proteins in their study had the most notable increase/decrease in subjects aged 60 or more.
Having age-dependent RI of tear cytokine levels can be very useful for the diagnosis of ocular disorders. Several very common eye problems are related to the phenomenon of inflammaging. For example, an association between age and DED has been assumed for a long time, and recently, Di Zazzo et al. [48] have directly related inflammaging to a progressive impairment of ocular surface system function. An ocular disease in which age is an important risk factor is age-related macular degeneration, affecting people aged 65 years and older. Several inflammatory gene polymorphisms have been linked to this disorder indicating an important role of inflammation and immune-mediated processes, and, therefore, directly related to inflammaging term [49]. Finally, low-grade inflammation has been described as a causal factor in the pathogenesis of glaucoma [50]. Although glaucoma can affect people of all ages, is most common in adults in their 70s and 80s, so it is accepted that glaucoma risk increases with age [51].
Sex also has an influence on tear levels of some of the analyzed cytokines. Particularly, Fractalkine tear levels were higher in males compared to females and IL-7, IL-8, and IP-10 higher in females compared to males; whereas, EGF, IL-1RA, MCP-1, and VEGF levels were not affected by sex. Some other studies have also shown sex-dependence in cytokine levels in other fluids. For instance, Nasi et al. (2006) [42], found a sex-dependence for IL-7 serum levels, with higher levels in centenarian females than in males, in agreement with our results in tears. Moreover, the study by Nättinen et al. (2019) [16] revealed that protein-age correlations in tear film were, in most cases, more significant in males, although mammaglobin-A (SCGB2A2) protein had a higher correlation coefficient for females. Changes in sex hormones have been suggested as connected with these differences. In addition, very recently, stronger relationships between age and IL-6 and high-sensitivity C-reactive protein plasma levels have been found in females compared to males [52]. Although the cytokines analyzed in that work are not the same that those in this study, our results reflected the same behavior in relation to IL-7, IL-8, and IP-10 pro-inflammatory cytokines.
The association between tear concentration of cytokines and DED has been widely reported (reviewed in [3,4,53,54]). However, it is often difficult to use cytokines as diagnostic tools due to the lack of diagnostic cut-offs. In this study, we have shown the utility of our RIs to diagnose DED. All of the severe DED patients in our sample showed abnormal tear levels of at least two of the evaluated cytokines when compared to the established RIs. Fractalkine, IL-1RA, IL-8, and VEGF showed high levels in more than 70% of the DED patients. Moreover, EGF and IP-10 were found to be infra-expressed in 23.1% and 38.5% of patients respectively of our DED sample. Altered tear levels of all these molecules have already been reported in different clinical subgroups of patients with DED. So, for example, tear concentrations of Fractalkine, IL-1RA, IL-8, and VEGF have frequently been found significantly increased in DED patients as compared to normal controls [8,18,[55][56][57][58][59], whereas tear EGF level has been found significantly decreased in the more severe forms of DED [55,60,61]. IP-10 tear levels were also found decreased compared to normal in 38.5% of the DED patients' sample, in agreement with results in tears from other studies involving severe cases of DED associated to ocular graft versus host disease [62], Stevens-Johnson syndrome [63], and, as preliminary findings suggest, in primary Sjögren's syndrome subjects [64,65]. However, several authors have also found tear IP-10 concentration to be overexpressed in tears of patients with moderate evaporative DED [21] and with Sjögren's syndrome [66]. It has been suggested that the anti-angiogenic properties of tear IP-10 may contribute to normal ocular surface immune privilege [65] and the low levels in Stevens-Johnson syndrome have been suggested to lead to a severe smoldering inflammation and progressive subconjunctival fibrosis in these patients [63]. Further studies with a bigger DED sample are warranted to evaluate the role of IP-10 in DED.
A critical issue of establishing RIs is to recruit a valid group of reference individuals. Two fundamental aspects are important in this regard: the definition of "reference individual" and the sample size. In this case, a reference individual should be a subject without eye pathologies. However, the characterization of this type of subject is not always easy, so collecting medium/large samples from the reference population can be quite expensive. Certainly, although above 150 individuals, our sample could be considered not sufficiently large to properly represent the reference population. Nevertheless, it is very important to highlight that a single ophthalmologist evaluated all subjects included in our sample, so the eye health criterion was homogeneous. Moreover, sample-handling factors such as tear collection, storage, or processing protocol may be a strong influence on the value of cytokine level and they could help to invalidate the sample as a reference sample. In this study, the tear samples were collected and processed by a single person, so this source of variability was also controlled. It also should be considered that in this study environmental conditions during subjects' clinical evaluation and sample collection were not controlled; it is well known the fact that ocular surface parameters and tear mediator concentration are altered under adverse environmental conditions, particularly low humidity and/or high temperatures. Nevertheless, all subjects were evaluated and samples were obtained in the same air-conditioned rooms (in a normal range of temperature and humidity), assuring that all of them were evaluated under the same conditions. In general, aware that it will not always be possible to collect samples under controlled environmental conditions, our recommendation is to at least do it in the same office, under similar ambient conditions. Additionally, it is noted that our validation, based on an external sample of subjects without ocular pathology, was appropriate for all of studied cytokines. As per clinical and laboratory standards institute C28-A3 guidelines [28] validation can be considered successful, since more than 20 samples were evaluated and more than 90% of values fell within the RIs. In some way, this result reinforces the validity of our reference sample. Finally, although this study is related to the use of specific cytokine/chemokine X-MAP-Luminex kits for molecule determination, the results obtained about the dependence on sex and age in tear cytokine levels would remain the same regardless of the use of kits from other vendors or another molecule determination assay.

Conclusions
We have estimated adjusted RIs of eight tear cytokine levels by GAMLSS regression models. We were able to determine the effect of age and sex in the concentration of each of these tear molecules in subjects without ocular disorders and adjusted their cut-off values taking into account these dependencies. Finally, we have showed the effectiveness of our RIs as diagnostic tools for severe DED. These age-and sex-intervals for tear cytokine levels could therefore be used as reference values in other tear cytokine studies in patients, and can be an important tool to extend the use of tear cytokines as biomarkers of ocular diseases.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Residual Analysis
This appendix shows the diagnostic tools based on normalized quantile residuals [25] used for checking the adequacy of each one of the fitted generalized additive models for location, scale and shape (GAMLSS) models. For a particular model, the following diagnostic tools are used: If the residuals behave well, the density estimate of the residuals will be approximately normal, the normal Q-Q plot will be approximately linear (with intercept 0 and gradient 1) and the bottom two or three plots will show a random scatter around the horizontal line at 0, and similar distribution of residuals for both sexes. In these plots line on Y = 0 (solid line) and fitted smooth curve (dotted line) are also represented. The shape of cubic fitted curve to the points (solid curve) shows inadequacies in the model described in Table A1.
Multiple worm plots are also built in order to identify failures of the model within different ranges of the explanatory variables. Age is cut into four non-overlapping intervals with similar number of observations and sex splits the data into female and male groups.

Appendix B. Centile Curves Estimation
This appendix shows the estimations of centile curves for each evaluated cytokine based on the fitted generalized additive models for location, scale and shape (GAMLSS) models.

Appendix B. Centile Curves Estimation
This appendix shows the estimations of centile curves for each evaluated cytokine based on the fitted generalized additive models for location, scale and shape (GAMLSS) models.  Table A6. Percentiles of IL-7 levels estimated with the fitted GAMLSS model. Family distribution Power Exponential with age dependence modelled as a cubic spline for location and skewness parameters (µ and ν respectively) and as linear for scale parameter (σ). Sex is added as a qualitative variable for location and skewness parameters (µ and ν respectively).

Appendix C. Reference Intervals Classification
This appendix shows the complete classification of all individuals of validation and utility samples based on the cytokine reference limits estimated by the fitted generalized additive models for location, scale and shape (GAMLSS) models.