Original research

Prediction of refractive error and its progression: a machine learning-based algorithm

Abstract

Objective Myopia is the refractive error that shows the highest prevalence for younger ages in Southeast Asia and its projection over the next decades indicates that this situation will worsen. Nowadays, several management solutions are being applied to help fight its onset and development, nonetheless, the applications of these techniques depend on a clear and reliable assessment of risk to develop myopia.

Methods and analysis In this study, population-based data of Chinese children were used to develop a machine learning-based algorithm that enables the risk assessment of myopia’s onset and development. Cross-sectional data of 12 780 kids together with longitudinal data of 226 kids containing age, gender, biometry and refractive parameters were used for the development of the models.

Results A combination of support vector regression and Gaussian process regression resulted in the best performing algorithm. The Pearson correlation coefficient between prediction and measured data was 0.77, whereas the bias was −0.05 D and the limits of agreement was 0.85 D (95% CI: −0.91 to 0.80D).

Discussion The developed algorithm uses accessible inputs to provide an estimate of refractive development and may serve as guide for the eye care professional to help determine the individual best strategy for management of myopia.

What is already known on this topic

  • Percentile curves are currently used to statistically compute the further progression of (spherical equivalent) refractive errors and axial length and for this approach, a big data approach is needed to get a distribution that represents a population.

What this study adds

  • The conduction of large cross-sectional studies or the analysis of retrospective data of a given population is time and cost intensive. The study investigated how machine learning-based algorithms perform in the computation of refractive errors over time, based on different input variables.

How this study might affect research, practice or policy

  • The use of machine learning to research and develop algorithms for the prediction of refractive errors will help to generate algorithms for a high-quality prediction of refractive development with only a small subset of data needed compared with the situation before.

Introduction

Already in 2012, the WHO declared visual impairment as major health issue, while uncorrected refractive errors cover 43% of these impairments.1 Since then, uncorrected refractive errors still appear on the list of major causes for visual impairments.2–5 Myopia shows considerable variation in prevalence among children of different ethnic origins, geographic locations and ages.6 A recent review on the prevalence of myopia on children worldwide pointed out a myopia prevalence of 60% in Asiatic countries, with the eastern region being the most affected with prevalence values of 73%. In European countries, the study showed a 40% of myopia prevalence, while 42% in North American children. On the contrary, African and South American countries showed myopia prevalence under 10% in children.7 Myopia prevalence is also estimated to increase worldwide within the next decades.8 9

Some biometric components of the eye influence refractive errors. Among the most common are the corneal shape, the crystalline lens shape and thickness and the axial length of the eye. Regarding the axial length, myopia may be developed when the eye grows to an extent that overruns the focal length of the eye,10 in case of relaxed accommodation. In some cases, when the length of the eye is superior or equal to 26 mm, additional complications must be considered, such as glaucoma, retinal detachment or myopic maculopathy.11–13 The risk for these complications might be reduced by evaluating the risk of myopia onset and accordingly, applying an appropriate and prompt management. The age of myopia onset and its progression rate are of great importance for the evaluation of high-myopia risks and eventual complications.13 14 A wide number of variables that influence its onset and progression has been suggested, but only a hand full of them are frequently named. Among these parameters are the gender, age, ethnicity, behaviour and parental myopia.2 14–21 Regarding proper and prompt management, several options have been identified to delay myopia onset or slow myopia progression, such as more time spent outdoors, less near-work time, topical atropine treatments, optical interventions, among others.21–25 To better assess and predict the individual risk of myopia development, different strategies have been developed and explored, for example, percentile growth curves or artificial intelligence for growth curve estimation. Recently, growth curves of refractive power and axial length have helped to reference kids on different populations to make an estimate on the eventual development.13 26–30 The use of artificial intelligence in the prediction of refractive error prediction might give benefits when compared with the statistical approach by using percentile curves, but needs research and development as indicated by Foo et al,31 therefore, the present study investigates the application of artificial intelligence for the prediction of refractive errors and discusses its implications in the fields of optometry and ophthalmology.

Methods

Data and subjects

The study was a community-based cross-sectional data collection. Every child available was included in the study. Participants signed an informed consent before the data collection and all procedures performed were in accordance with the Declaration of Helsinki. Patients or the public were not involved in the design, or conduct, or reporting, or dissemination plans of our research.

The data was collected by the Wuhan Center for Adolescent Poor Vision Prevention and Control in Wuhan, China. Data from 12 780 kids that consisted of eye information on spherical power (spherical component of the refractive error) and biometrical parameters, together with age and gender was used. The biometrical parameters included the axial length, corneal radius (mean of the flattest and steepest radii), lens thickness and anterior chamber depth. Spherical refraction was collected using the Topcon CV-3000 autorefractor (Topcon, Tokyo, Japan). All autorefractive measurements were performed while subjects were under cyclopaedia (cyclopentolate at 0.5%, applied four times each for 5 min). Biometric data was collected using the Lenstar LS 900 SN 1914, V.1.1.0 (Haag-Streit AG, Koeniz, Switzerland). The data was separated into two, a cross-sectional (12 554 kids) and a longitudinal (226 kids), data sets. The cross-sectional data set was divided into 48.2% of girls and 51.8% of boys and their age ranged from 5 to 16 years. Their spherical refraction was −0.93±1.85 D (SD) for the girls and −0.88±1.83 D (SD) for the boys, while their age was 9.99±2.47 years (SD) and 9.90±2.48 years (SD), respectively. The longitudinal data set was divided into 40.7% girls and 59.3% boys and measurements were taken in three separate visits. Each dataset was used for the training of separate regression models. For the validation, a separated longitudinal data set of 81 kids that consisted of the same named parameters was used.

Development of algorithm and modelling

The development of the complete prediction algorithm and the modelling was performed using MATLAB, V.R2018a (MathWorks, Natick, Massachusetts, USA). In the specific case of modelling, supervised machine learning from the Regression Learner Application found at the Statistics and Machine Learning Toolbox V.11.3 was used.

The algorithm was defined as the set of code lines and models that in combination allowed the spherical refraction prediction. The intention was to take as inputs only the age, gender and spherical power and deliver as output the spherical power as a function of the age. This definition ensured the development of an algorithm that may be deployed in several applications without the need to measure many eye variables, that is, biometric parameters. The same applied for the models trained and used on the algorithm. To keep the models with the less possible independent variables, only variables that were significant and of common use for eye care professionals were taken. The initial independent variables evaluated to be fitted on the models were the age, gender, axial length of the eye, corneal radius, anterior chamber depth and lens thickness. However, only the age, gender, spherical power and axial length to corneal radius ratio (Axl/Cr) were defined as commonly used variables.

Statistical analysis and validation

The statistical evaluations were performed using MATLAB, V.R2018a (MathWorks, Natick, Massachusetts, USA). Stepwise multiple regressions were performed to identify whether the independent variables were of significance for the model. The performance of the prediction model was evaluated using Bland-Altman plots to evaluate the bias and limits of agreement. Moreover, scatter plots and the Pearson correlation coefficients (r2) between prediction and measured data were also evaluated.32

In addition, a validation was performed that compared results of the developed algorithm with the prediction of two already available solutions for the prediction of myopia, namely the Myopia Calculator (Brien Holden Vision Institute, Sydney, Australia) and the MyAppia (MyopiaCare by Eyetific, Renens, Switzerland). The statistical test used to evaluate for normal distribution and to compare the distributions between the solutions were performed using the Kolmogorov-Smirnov test and the Wilcoxon signed-rank test from Matlab. Differences were considered statistically significant when the associated p values were lower than 0.05.

Results

Selection of variables

A forward stepwise regression was performed to evaluate the predictors eligible for the trained models. The evaluation added to the trained model variables that would fulfil a threshold p value of 0.06. A separated evaluation was performed using the cross-sectional and the longitudinal data sets. For the cross-sectional data set, two dependent variables were tested, namely the sphere (D) and the Axl/Cr. When the sphere was the dependent variable, the independent variables inputted were the age, gender, Axl/Cr, anterior chamber depth and lens thickness. In case the Axl/Cr was the dependent variable, the independent variables inputted were the age, gender, sphere, anterior chamber depth and lens thickness. For the longitudinal data set, the procedure was the same as when using the cross-sectional data, the only difference was that the output dependent variable was also added as the input independent variable measured at baseline. For the forward stepwise regression using the cross-sectional data, the evaluation showed that all inputted parameters were of significance for both tested models. For the forward stepwise regression using the longitudinal data, the evaluation showed that when the output was the Axl/Cr, the sphere, Axl/Cr and anterior chamber depth were of significance for the tested model. When the output was the sphere, the Axl/Cr, the sphere and the lens thickness were of significance for the tested model. Specific values for the F-Stat and p values of all evaluations are provided in the online supplemental material S1.

The model was developed to allow eye care professionals to evaluate the risk of children regarding their onset and/or progression of myopia in children feasible input variables. This led to a final selection of age, gender, refractive error and Axl/Cr.

Prediction algorithm and validation of performance

After definition of the significant variables, a training of all 19 regressions available in Matlab V.R2018a (MathWorks, Natick, Massachusetts, USA) was performed for each differently defined combination of inputs and outputs. This was done for evaluation of the default performance delivered by Matlab for each regression and consequently, the best models were selected for each data set and input–output combination. Afterwards, from the selected regressions, a hyperparameter optimisation was performed including a separate performance evaluation using a separated dataset to avoid overfitting. The final model selection was done when the validation using the separate dataset showed values for the coefficient of determination (R2) and for the root mean squared error (RMSE) that were comparable to the values delivered after training. The details on the hyperparameter optimisation can be found in the online supplemental material, S2.

For a prediction algorithm which takes as input the age, gender and sphere and outputs the sphere as a function of the age, different combinations of regressions and data sets were tested. Tables 1 and 2 show the best model that resulted from the training of cross-sectional and longitudinal data, using different inputs and outputs. The column that shows the best regression on tables 1 and 2 specifies the name of the model together with the coefficient of determination (R2) and RMSE delivered after training.

Table 1
|
Best regression model obtained for the use of cross-sectional data and specified input and output
Table 2
|
Best regression model obtained for the use of longitudinal data and specified input and output. The table shows Support Vector Machines SVM as the best model.

A total of five models resulted from the training of cross-sectional and longitudinal data. Whereas for the models trained using cross-sectional data, the output was obtained for the same input age, for the models trained with longitudinal data the output was obtained for a future age (age+n). From the five models, only model 5 fulfilled the definition of input and output (input: age, gender and sphere and output: sphere at age+n). Thus, the models from 1 to 4 were combined in such a way that this definition was fulfilled.

The performance of the combined models and of model 5 was evaluated using the longitudinal data of 81 subjects measured two times and that were not included in the training. Figure 1 shows the input, algorithms and output that resulted from the combination of the trained models. Additionally, the right column shows the validation values obtained using the longitudinal data.

Figure 1
Figure 1

Input, algorithms and output resulting from the combination of the trained regression models. Additionally, the right column shows the validation values r2 for Pearson correlation coefficient, Bias and limits of agreement, obtained for each algorithm. Axl, axial length; Cr, corneal radius; R2, coefficient of determination; RMSE, root mean squared error; SVR-linear, linear support vector regression; SVM-quadratic, quadratic support vector machine; GPR- Sq exponential, square exponential gaussian process regression; GPR- Matern, matern gaussian process regression;

The algorithms are listed from top to bottom, starting with the algorithm that showed the lower Pearson correlation coefficient (r2) on the validation and ending with the algorithm that showed the highest r2. Additionally, information on the algorithm bias and limits of agreement are shown. Regarding the bias, all values improved as the r2 also improved, excepting from the algorithm at the bottom that showed the second highest bias. This could be explained due to the use of only one model and the relatively low number of data used on the training of this model (226 kids). This effect was a decision criterion for leaving this algorithm out of the options.

From the three remaining algorithms, the combination of linear support vector regression and Matern Gaussian process regression resulted with the lowest bias and best limits of agreement. More details on this performance evaluation can be seen on figure 2. The figure shows on the left side the prediction as a function of the measured true data for the two available data points of each kid (n:162). The Pearson correlation value between prediction and measurement was r2: 0.77, whereas the total sample size of 162. On the right side, the Bland-Altman plot shows the difference between prediction and measurement as a function of their mean. The obtained bias value of −0.05 D (p:0.14) is shown and the limits of agreement for all ages pooled of 0.85 D (95% CI: −0.91 to 0.80D) is also shown.

Figure 2
Figure 2

Validation diagrams for the selected algorithm that combines an linear support vector regression model with a Matern Gaussian process regression model. The left diagram shows the prediction as a function of the measured data for the total sample size of 162. Pearson correlation value r2 and the sample size of 162 data points are also shown. The right diagram shows the difference between prediction and measurement as a function of their mean. Additionally, a bias of −0.05 D (p:0.14) and limits of agreement of 0.85 D (95% CI: −0.91 to 0.80) are shown.

A slight refraction-dependent bias can be observed on the Bland-Altman plot. For mean refractive errors approaching emmetropic values, the difference between prediction and measurement seems to be higher than for myopic values ranging from −2.00 to −4.00 D. For an eventual application of the algorithm, this kind of bias can be compensated by fitting a linear regression on the data set and applying the slope value on the predictions.

Comparison to available solutions

The selected algorithm was compared with the Myopia Calculator (Brien Holden Vision Institute, Sydney, Australia) and the MyAppia (MyopiaCare by Eyetific, Renens, Switzerland). The longitudinal data set of the 81 kids used for the performance evaluation was also used to generate predictions of refractive power using both solutions. The absolute difference between the predicted and measured true data was calculated. Figure 3 shows boxplots for the absolute differences obtained for each solution. On each box, the central line indicates the median absolute error, and the box indicates the 25th and 75th percentile. The error bars are extensions to the most extreme data points that are not considered outliers, whereas the outliers are marked by the ‘+’ symbol.

Figure 3
Figure 3

Absolute difference on prediction obtained for the Myopia Calculator, MyAppia and the support vector regression and Gaussian process regression (SVR&GPR) algorithm. The central line indicates the median value, whereas the boxes indicate the 25th and 75th percentiles.

To further compare the algorithm with the two available solutions and given that under absolute values the distributions were not normal (Kolmogorov-Smirnov test, p<0.001 for all solutions), a Wilcoxon signed-rank test was performed. The test showed significantly different distributions between the algorithm and MyAppia (p=0.89), whereas there were no significantly different distributions between the algorithm and the Myopia Calculator (p=0.03).

Discussion

The purpose of this study was to develop an algorithm for the prediction of spherical power progression, using the less possible number of inputs without hindering an acceptable performance. Although all the available data resulted of significance for the prediction of refractive development, the training variables selected were the age, sphere, gender and Axl/Cr. Moreover, the resulting algorithm with the best performance allowed the use of the age, sphere and gender as inputs. This arrangement may lead to a flexible application, without the need to collect all available biometrical parameters of the eye.

Resulting algorithm

A stepwise multiple regression evaluation delivered the variables of significance for the model. Similarly, Magome et al33 also used stepwise multiple regression to evaluate the significance of biometrical parameters on the prediction of spherical power and found that axial length, anterior chamber depth, lens thickness, corneal refractive power and corneal astigmatism acted as significant parameters. In this study, when using the longitudinal data set of 226 children, anterior chamber depth and lens thickness were not of significance for the tested model and these differences in significance might be related to the differences in the data sets. As the authors evaluated eyes of 617 children from which the SD of the provided ocular biometry was <0.5%, the number of subjects and stability of measurements might be responsible for these differences.33

A combination of linear support vector regression and Matern Gaussian process regression delivered the best results in terms of bias (−0.05 D) and limits of agreement (±0.85 D). The bias obtained can be considered of no significance for the optometric practice, since eye care professionals assess the refractive powers in ±0.25 D steps; however, the limits of agreement of ±0.85 D should be improved as more data is available.

Lin et al34 used refraction data of 129 242 individuals and trained a random forest model for the prediction of high myopia onset. The validation of the model suggested a clinically acceptable prediction of high myopia at a future point in time, namely, at age 18 years. However, for earlier ages like 8 years, 95% of the predicted values differed from the true value within 0.5–0.8 D. The difference to our combination of models might be explained by the difference in the amount of data available. Whereas the authors had 129 242 individuals and a total of 687 063 multiple visits records,34 our data consisted of 12 780 cross-sectional data points and data points of 226 kids measured either two or four times. Additionally, the large number of data points allowed the training of a single model.34 This was not applicable for the current, as the data availability was a limitation that forced the training of two separated models for the estimation of refractive error in a future point in time. The combination of two models instead of a single model might explain the differences in prediction performance between these two solutions. Future attempts should try to simplify the model into one, to avoid risks of overfitting or contradictory predictions between models.

Other studies, such as Rampat et al,34 investigated the prediction of subjective refraction using aberrometry eye data from 3.729 individuals and machine learning. The results suggested an acceptable prediction of subjective refraction from polynomial wavefront data. However, unlike in the current study and by Lin et al34 and Rampat et al,35 authors studied refractive development always for a present point in time and not for a projected time. To our best knowledge, the presented study and the study by Lin et al34 study showed for the first time that big data and machine learning can be used to support prediction of myopia prognosis of Asiatic children on a future point in time.

Comparison with available solutions

The performance of the algorithm was compared with the Myopia Calculator (Brien Holden Vision Institute, Sydney, Australia) and the MyAppia (MyopiaCare by Eyetific, Renens, Switzerland). A slightly better performance was obtained when the algorithm was compared with the MyAppia and similar performance when the algorithm was compared with the Myopia Calculator. Saunders et al36 evaluated the prediction of the Myopia Calculator using available data from 80 European kids. The results showed that for kids with age between 9 and 10 years, 58% of the predictions were overrated as myopic, whereas 26% were categorised erratically. For older kids aged between 12 and 13 years, the progress of their refractive power was more likely to be in accordance with the prediction. They concluded that overestimations were more likely for younger kids with lower myopic error and suggested to handle predictors with care. In this study, longitudinal data of 81 Asiatic kids were used for the comparison of the three solutions; however, the dependency of prediction on age could not be evaluated due to the small number of data sample used for the comparison and their inhomogeneous age distribution. Another difference worthwhile to mention is that whereas the MyAppia and Myopia Calculator can only take as input myopic refractive data, our solution can also consider data from hyperopic kids. This option may be of help for the eventual prognosis of myopia onset.

Limitations

The validation of this study was performed using longitudinal data of 81 kids. The age and refractive power distribution of the 81 kids was not homogeneously distributed, leading to the only alternative to evaluate the error on prediction by pooling for all ages available. Future developments should validate for each age, separately. In addition, the training was performed with a combination of cross-sectional and longitudinal data. For a better performance on the prediction of spherical power development, big data on longitudinal measurements might be used.

An additional limitation was the fact that refractive power progression depends on several more parameters, for example, parental myopia, near-work time and time spent outdoors. Thus, algorithms that do not consider these variables must be taken with specific care. Furthermore, the pace of refractive development and/or progression is still not fully understood and may vary among kids. Hence, although growth curves and prediction algorithms might be useful to reference the kid within a population, they must be used only as a reference and not as a unique tool for risk assessment.

Moreover, although artificial intelligence may be gaining strength in healthcare fields such as ophthalmology,31 the by-products of this technique must be handled with care. As time passes, different circumstances influence the lifestyle of people, and this might lead to changes in visual health. Wang et al37 showed how the COVID-19 confinement affected the myopia development and its onset, accounting for a significant myopic shift on Chinese children aged between 6 and 8 years. Additionally, technological changes and always more common use of devices on the eye’s proximity might also lead to impacts on visual health. McCrann et al38 studied the use of smartphones as a possible risk factor of myopia and found that myopic participants used almost two times more smartphone data than non-myopes. Social changes affecting visual health raise the need to gradually collect data and constantly update the prediction models.

Conclusion

The present study showed that using machine learning and large data of Chinese kids, an algorithm for the prediction of spherical power as a function of the age could be developed. The algorithm that showed an acceptable performance consisted of a support vector regression and a Gaussian process regression. The performance evaluation covered a correlation between prediction and measured true data of 0.77, a bias of −0.05 D and a limit of agreement of ±0.85 D. Moreover, the inputs selected for the algorithm, that is, age, gender and spherical power, are inputs which can normally be accessed by eye care professionals or parents. This may allow a flexible application of the algorithm and might support eye care professionals to evaluate the prognosis of refractive errors on kids.