LMS parameters, percentile, and Z-score growth curves for axial length in Chinese schoolchildren in Wuhan

Understanding the ocular structural changes are fundamental to defining strategies for myopia prevention and management. This study aimed to establish age-gender specific normative LMS parameters for axial length to generate percentile and Z-score growth curves in a population of Chinese schoolchildren. A total of 14,760 individuals aged 6 to 15 years from Wuhan, central China, contributed to this study. The LMS method was used for the calculation of LMS parameters and the generation of percentile and Z-score growth curves for axial length. Growth curves derived from the LMS parameters were compared with those originally calculated. Axial elongation was age- and percentile-dependent. The highest elongation rate occurred at the 98th percentile in the range 6 to 9 years, being up to 1.46 mm in boys and 1.42 mm in girls. The largest differences between original and newly generated growth curves were detected at the 98th percentile at age 15; 0.78 mm (females) and 0.63 mm (males). Multinomial logistic regression and receiver operating characteristic analyses revealed Z-scores as a good predictor for estimating high myopia development. The axial length growth curves presented in this study provide a technically solid instrument that depicts the best description of physiological eye growth for Chinese schoolchildren aged 6 to 15 years.

and later extended by Cole and Green 18 . The LMS method is based on the use of Box-Cox transformations 21 to normality through the calculation of a skewness parameter 22 . Given its efficiency and flexibility it has been applied as a standard method in many benchmark studies 1,[23][24][25] . Particularly, LMS technique has not been previously used for the generation of growth curves for ocular variables, such as axial length. Therefore, in this paper, we aim at demonstrating the application of the LMS method for population-based ocular axial length data. We report the basis for the calculation of LMS parameters and the subsequent estimation of axial length growth curves in a Chinese school-age population. In addition, the new percentile growth curves derived from the LMS parameters are compared with the original percentile curves previously calculated and published from the same study population.

Results
LMS parameters. Table 1 shows the different combinations of edf values chosen for the modelling and optimization process of the new calculated percentile curves. The edf values were increased and reduced one at a time until the smallest SBC value was found. First, we set the initial edf values for L to 0 and for M and S to 1. Considering edf (M) > edf (S) > edf (L), edf values for M were increased first, followed by edf values for S. While higher edf values for M and S generated smaller SBC values, greater edf values for L produced higher SBC values (Table 1). For both genders, the best model optimization was achieved with edf values of L = 0, M = 4, and S = 2, as they provided the highest reduction in SBC and a reasonable smoothing. As seen in Table 1, a subsequent increase in edf values did not lead to a reduction in SBC values. Tables 2 and 3 show the LMS parameters for axial length. In both females and males, value of 1 was used to represent the L curve at all ages. M curve increased with age in both genders. On average, males presented higher but not significantly different than females (average difference = 0.54 ± 0.04 mm; p = 0.10, t-value = − 1.75, df = 18; two-sample t-test). S curve increased with age in both genders, reaching the highest values at age 15. Both groups showed similar S curves (p = 0.81, t-value = 0.24, df = 18; two-sample t-test).
Percentile growth curves derived from LMS parameters. Figure 1a and Table 2 exhibit the new percentile growth curves for axial length generated using the LMS parameters. In both genders, all percentiles trended upward with age. The higher the percentile rank, the greater the increase in axial length. From 6 to 15 years, the 2nd percentile had a total axial elongation of 1.59 mm (females) and 1.71 mm (males), while the 98th increased by 2.35 mm (females) and 2.46 mm (males). Concretely, in the female group, the 10th percentile increased by 1.74 mm (21.47  Percentile axial elongation rates trended positively but unevenly with age. From 6 to 9 years of age, in females, the elongation rate between the 2nd and 98th percentiles ranged from 0.97 to 1.42 mm (4.67-5.85%), while in males ranged from 1.04 to 1.46 mm (4.87-5.93%). In this age range, those percentiles below and above the median showed a significant increase in axial length (two-sample t-test: p = 0.03, t-value = − 2.89, 95% CI [− 1.92  Table 3. LMS parameters and axial length Z-scores in millimeters, as a function of age and gender.  www.nature.com/scientificreports/ www.nature.com/scientificreports/ to − 0.16], for females; and p = 0.02, t-value = − 3.00, 95% CI [− 2.00 to − 0.20], for males). From 9 years onwards, the older the child, the lower the axial elongation rate. In the range 9-12 years, the elongation rate between the 2nd and 98th percentiles ranged from 0.40 to 0.60 mm (1.81-2.35%) in females. In males, ranged from 0.45 to 0.65 mm (2.00-2.50%). Compared to the 6-9-year range, this meant an average reduction in elongation rate of 61% in females and 58% in males at all percentiles. In the range 12-15 years, the elongation rate between the 2nd and 98th percentiles ranged from 0. 22  Comparison with original growth curves (Sanz Diez et al. 10 ). Figure 1b shows the comparison between the originally published axial length growth curves 10 and those derived from the LMS parameters. In females, between 6 and 15 years, the new percentile cut-off values ranged from 20.97 to 27.32 mm, while the originals ranged from 20.84 to 26.55 mm. The maximum discrepancies between the original percentiles and those newly generated appeared at the 95th and 98th percentiles at 15 years, where the differences were equal to 0.63 mm and 0.78 mm, respectively. Along the 50th and 90th percentiles, the highest differences between axial length curves were also found at 15 years of age, being 0.13 and 0.15 mm respectively. Along the 2nd, 5th and 10th percentiles, the maximum differences were on average 0.13 mm within the age range of 6 to 7 years. The minimum differences occurred along the 25th, 50th and 75th percentiles (mean difference: 0.05 ± 0.01 mm). In females, no significant differences between original and newly generated percentile curves were found (p > 0.77, t-value = [− 0.12, 0.29], df = 18; all superimposed percentile curves, two-sample t-test).
In males, between 6 and 15 years of age, the new percentile cut-off values ranged from 21.27 to 27.16 mm, while the originals ranged from 21.44 to 27.78 mm. The greatest discrepancies also appeared at 15 years at the 95th and 98th percentiles. Axial length differences were 0.54 mm for the 95th and 0.63 mm for the 98th percentile. Along the 2nd, 5th and 90th percentiles, the largest differences were also found at the age of 15 years (mean difference: 0.23 ± 0.01 mm). The smallest differences were also observed at the 25th, 50th and 75th percentiles (mean difference: 0.05 ± 0.03 mm). In males, no statistical differences were found between the original curves and those derived from the LMS parameters (p > 0.62, t-value = [− 0.02, 0.50], df = 18; all superimposed percentile curves, two-sample t-test).
Z-scores. Z-score curves are given for axial length based on age and gender (Fig. 1c, Table 3). From Eq. (2), the age-gender-specific LMS parameters provided in Tables 2 and 3 allow to calculate the Z-score corresponding to individual child's axial length data. For example: a 7-year-old boy with an axial length of 24.51 mm would have a Z-score of 1.25. A 10-year-old girl with an axial length of 22.75 mm would have a Z-score of − 1.26. Accordingly, using age, gender, and axial length data from the second cross-sectional dataset, the corresponding Z-scores were calculated and graphically superimposed over the reference Z-scores curves for both gender groups (Fig. 2a). As depicted in Fig. 2a,b, there is a clear distribution of Z-scores according to the refractive status and age. From 6 to 15 years, in females, 82.35% of hyperopes and 64.71% of emmetropes were below the zero Z-score (0.0 SD or median) while 60% of myopes and 100% of high myopes were above it. Particularly, in females, the refractive values of those myopes located below (− 1.49 ± 0.88D) and above (− 2.49 ± 1.36D) the median differed statistically (p < 0.001, Mann-Whitney test). In males, 80.56% of hyperopes and 68.75% of emmetropes were below the 0.0 SD curve, whereas 64.14% of myopes and 90.91% of high myopes were above it. In males, refractive errors were also significantly different between those myopes below (− 1.60 ± 0.93D) and above (− 2.37 ± 1.29D) the 0.0 SD curve (p < 0.001, Mann-Whitney test).

Z-scores as a predictor variable.
To assess the predictive role of Z-scores on axial and refractive development, the longitudinal dataset (226 children) was employed. Age, gender, and axial length data at the first visit were used to determine the Z-scores. A linear regression analysis between Z-scores (first visit) and axial length values (third visit) revealed a significant direct relationship in females (F(1,85)  www.nature.com/scientificreports/ A multinomial logistic regression was performed to predict the refractive condition at third visit (hyperopia, emmetropia, myopia or high myopia) from the first visit Z-scores. The final model significantly predicted the dependent variable better than the intercept-only model alone in females (χ 2 (3) = 18.47, p < 0.001) and males (χ 2 (3) = 14.42, p = 0.002). Particularly, Z-scores were able to significantly discriminate high myopia from the other three types of refractive errors developed at the third visit (χ 2 (3) = 26.96, p < 0.001). For each one unit increase in Z-score, the log-odds of an individual falling into the hyperopia, emmetropia or myopia categories (relative to the high myopia category) is predicted to decrease by 3.12, 1.94 and 1.22 units, respectively. Moreover, as reported in the odds ratio (OR) column, with increasing Z-score values, the odds of falling into the hyperopia, emmetropia and myopia categories as changing by a factor of 0.044, 0.144 and 0.297 (Table 4). In both genders, those individuals classified as high myopes at the third visit were more likely to have higher Z-scores at the first visit (female group: β = 4.02, Std. error = 2.03, χ 2 (1) = 3.92, p < 0.05; male group: β = 1.46, Std. error = 0.65, χ 2 (1) = 5.07, p = 0.02).
The ROC analysis revealed Z-scores as a good predictive factor for high myopia, with an area under the curve of 0.84 for females (standard error: 0.04; 95% CI 0.76-0.92, p < 0.001) and 0.76 for males (standard error: 0.12; 95% CI 0.52-0.99, p = 0.03). As mentioned earlier for both genders, high myopia developed in most of those children whose Z-score (calculated from axial length) at the first visit was near to or above + 1.0 SD. In females, the Z-score cut-off value of + 0.88 SD showed sensitivity and specificity of 100.00% and 83.70%, respectively. In males, the Z-score cut-off value of + 0.91 SD exhibited sensitivity of 71.40% and specificity of 89.40%.

Discussion
In this paper we have illustrated the use of the LMS method for the development of axial length growth curves in a Chinese school-age population. The LMS method has been used to calculate the LMS parameters, from which the percentile and Z-score growth curves for axial length have been generated. Additionally, a comparison was made between the growth curves derived from the LMS parameters and the original curves calculated from the same study population 10 .
Interestingly, we have observed a close agreement between the percentiles derived from the LMS parameters and those originally developed 10 . Both methodologies showed comparable percentile trends with larger discrepancies at ages 14 and 15. Percentile curves have exhibited an increase in axial length with age. Particularly, the increase in axial length was percentile-dependent and inconsistent across different ages, with the youngest individuals increasing the fastest. These findings agree with the trends reported by other percentile growth studies. However, the rate of axial elongation differs considerably from those reported in European population studies 11,13,15 . Interestingly, the largest discrepancies between cohorts occur in the 6 to 9-year age range. In fact, this is consistent with the annual changes observed within this range, being on average 0.19 mm/year in the European cohort, while 0.41 mm/year in the Asian. Similar axial elongation rates have been reported by other authors, who have showed variances based on the refractive condition 8,9 . This emphasizes the importance of a close monitoring of the ocular components during the school years to carefully supervise possible structural and refractive changes 26,27 .
Given the LMS parameters calculated from the reference population, Z-score of any child can be calculated from their age and axial length. Z-scores showed a defined distribution pattern according to individuals' refractive condition and age. Considering the median as a reference, all hyperopic and emmetropic subjects were below it. In contrast, myopic and high myopic subjects were above it. Z-scores were found to be good predictors for high myopia. Multinomial analysis revealed Z-scores to be able to discriminate high myopia from the other refractive conditions. Moreover, ROC analysis revealed good performance values. These findings suggest the feasibility of Z-scores in the study of ocular component growth patterns, since are based on a reference population and can be studied as a continuous variable 16 .
Despite the reported findings, the main limitation of the study is the failure to consider myopia risk factors in the generation of the percentile and Z-score growth curves. Applying the effect of risk factors could result in improved accuracy. Another limitation is the extrapolation of results to other populations. The development and application of axial length growth curves should consider the geographic area and ethnicity of the study population. Axial length percentile growth curves for different geographic regions are needed to define global strategies for myopia prevention.
In summary, we provide a demonstration of the LMS method on ocular axial length data. While the LMS method has been widely applied in the anthropometric assessment for children and adolescents, its Table 4. Multinomial logistic regression model used to predict high myopia at third visit given Z-score from first visit. The reference category is high myopia. www.nature.com/scientificreports/ implementation has not been extended to the analysis of ocular axial length growth patterns. In this research we have established for the first time the LMS parameters for axial length, from which we have generated axial length growth curves comparable to those originally published. Additionally, we have observed the practicality and helpfulness of the LMS methodology for the calculation of axial length Z-score curves, which enable a more precise assessment of eye growth during childhood. We believe, the comprehensive study of ocular growth patterns in school children necessitates the use of valid methodologies to generate percentile growth curves and Z-scores that can shed light on the field of vision science. These findings may assist organizations and governments in evaluating and designing appropriate myopia prevention programs for children and adolescents.

Methods
Study population and data collection. The current study is based on a secondary analysis of the population from a previously published prospective cross-sectional study of school-aged children from Wuhan, China 10 . Eye data was collected by the Wuhan Center for Adolescent Poor Vision Prevention and Control. A total of 14,760 individuals (7133 girls and 7627 boys) were included in the study. Three datasets were used: two cross-sectional datasets and one longitudinal dataset. The first cross-sectional dataset was composed of 6054 girls (9.99 ± 2.47 years) and 6500 boys (9.90 ± 2.48 years) and was used to generate the growth curves for axial length, which served as a reference. The second cross-sectional set of data included 987 girls (9.36 ± 2.39 years) and 993 boys (9.32 ± 2.45 years) and assisted in showing the methodology implementation and in observing the distribution of generated growth curves as a function of refractive errors and axial length data. The longitudinal dataset consisted of 226 children (92 girls and 134 boys) with a total of three appointments over time and a mean time difference between the first and the third appointment of 2.67 ± 2.96 years for girls and 2.57 ± 2.77 years for boys. The longitudinal dataset was used to explore the efficiency of the newly calculated curves in predicting high myopic values.
Cycloplegic spherical refractive error data was obtained using the Topcon CV-3000 autorefractometer (Topcon, Tokyo, Japan). Axial length data was measured using a non-contact optical biometer, Lenstar LS900 (Haag-Streit AG, Koeniz, Switzerland). For more information on the study population and data acquisition, see Sanz Diez et al. 10 .
The study and data acquisition were approved by the Ethical Committee of the Wuhan Center for Adolescent Poor Vision Prevention and Control. Written and oral information was given, after which written informed consent was obtained from all participants or legal representatives. The study was conducted in accordance with the Declaration of Helsinki.
Percentile curves based on LMS method. Considering the purposes of the current study and the percentile range chosen in groundbreaking studies on axial length percentiles 10,11 , nine percentile curves (2nd, 5th, 10th, 25th, 50th, 75th, 90th, 95th and 98th) were generated. Here, percentiles were constructed through the LMS method, described by Cole [18][19][20] . The LMS method is a mathematical model that allows to fit longitudinal and cross-sectional anthropometric data, ocular axial length here, to obtain normalized percentile curves 20,21 . LMS method summarizes the distribution of the variable of interest according to age, based on three parameters or curves: L (λ), M (µ), and S (σ). These three parameters indicate the power in the Box-Cox transformation for the skewness adjustment (L), the median (M), and the generalized coefficient of variation (S) for each annual measurement. Assuming the residuals follow a normal distribution and by using the three estimated parameters (L(x), M(x) and S(x)) for each age (x), axial length values can be converted into percentiles and Z-scores. Percentiles at age (x) can be calculated from Eq. (1), where C 100α (x) is the 100 α -th percentile rank and Z α is the desired percentile in standard deviation units.
The corresponding Z-scores at age (x) can be obtained from Eq. (2), where y indicates the individual axial length measurement.
Based on the maximum penalized likelihood methodology, LMS curves can be fitted as cubic splines by nonlinear regression 18 . The degree of smoothness is defined by three smoothing parameters or equivalent degrees of freedom (edf). Each edf indicates the complexity of each L, M and S curve. A more detailed description of the LMS method can be found elsewhere 18-20 . Curve modelling. LMSchartmaker Pro software (version 2.54, Medical Research Council, UK) was used to calculate the LMS parameters and subsequently construct the axial length growth curves based on the LMS www.nature.com/scientificreports/ method. LMSchartmaker is a software developed by Pan and Cole to fit smooth percentile curves to reference data using the LMS method 28 . Under the software protocol, curve modelling involves selecting the appropriate age scale and choosing the proper edf values to optimize the L, M and S curves. Deviance measure is the benchmark for model fitting and curves optimization: the smaller the deviance measure, the better the optimization of the L, M and S curves. As a deviance measure, this study followed the authors' recommendations, and therefore Schwarz Bayesian Criterion (SBC) was used as a deviance measure for model fitting 28 . A more complete information of the software modelling statistics can be seen elsewhere 28 . Following Sanz Diez et al., axial length growth curves were computed considering both gender and age 10 . Statistical analysis. Statistical analyses were performed using the MATLAB R2020a statistical toolbox (MathWorks, USA) and SPSS statistical software package, v-27.0 (SPSS, Chicago, Illinois, USA). Data distribution was inspected visually (frequency distribution and quantile-quantile plot) and statistically (Kolmogorov-Smirnov test). Statistical analyses were done using the appropriate tests depending on the data distribution. Differences were considered statistically significant when p < 0.05. Results are provided as mean ± standard deviations. Multinomial logistic regression was performed, with refractive state condition as dependent outcome variable. Cox and Snell's, Nagelkerke's and McFadden's goodness-of-fit tests were used to evaluate the fit of the model. Likelihood ratio tests were used to assess the contribution of the independent variable to the model. Receiver operating characteristic (ROC) analysis was conducted to assess the diagnostic performance of Z-scores as a prediction model. ROC area of 1.0 indicates an ideal test, while an area of 0.5 describes an inaccurate test.
Comparisons of percentile curves of axial length data between original growth curves 10 and growth curves derived from the LMS method were performed in two ways. First, a comparison of the original and the newly calculated percentiles was prepared by means of a graphical superimposition of each of them. Second, differences between the percentile cut-off values of both growth curves were assessed by the Student's t-test at all ages.

Data availability
The datasets generated during the current study are available from the corresponding author on reasonable request.