Prediction of Subjective Refraction From Anterior Corneal Surface, Eye Lengths, and Age Using Machine Learning Algorithms

Abstract

Purpose: To develop a machine learning regression model of subjective refractive prescription from minimum ocular biometry and corneal topography features.

Methods: Anterior corneal surface parameters (Zernike coefficients and keratometry), axial length, anterior chamber depth, and age were posed as features to predict subjective refractions. Measurements from 355 eyes were split into training (75%) and test (25%) sets. Different machine learning regression algorithms were trained by 10-fold cross-validation, optimized, and tested. A neighborhood component analysis provided features’ normalized weights in predictions.

Results: Gaussian process regression algorithms provided the best models with mean absolute errors of around 1.00 diopters (D) in the spherical component and 0.15 D in the astigmatic components.

Conclusions: The normalized weights showed that subjective refraction can be predicted by only keratometry, age, and axial length. Increasing the topographic description detail of the anterior corneal surface implied by a high-order Zernike decomposition versus adjustment to a spherocylindrical surface is not reflected as improved subjective refraction prediction, which is poor, mainly in the spherical component. However, the highest achievable accuracy differs by only 0.75 D from that of other works with a more exhaustive eye refractive elements description. Although the chosen parameters may have not been the most efficient, applying machine learning and big data to predict subjective refraction can be risky and impractical when evaluating a particular subject at statistical extremes.

Translational Relevance: This work evaluates subjective refraction prediction by machine learning from the anterior corneal surface and ocular biometry. It shows the minimum biometric information required and the highest achievable accuracy.

RESUMEN:

Objetivo: El desarrollo de un modelo de regresión de aprendizaje automático prescripción refractiva subjetiva a partir de las características mínimas de la biometría ocular y la superficie corneal.

Métodos: Los parámetros de la superficie corneal anterior (coeficientes de Zernike y queratometría), además de longitudes axiales y de cámara anterior, edades y las refracciones subjetivas no ciclopléjicas de 355 ojos se dividieron en un conjunto de entrenamiento (75%) y otro de test (25%) y se entrenaron diferentes algoritmos de regresión de aprendizaje automático mediante validación cruzada 10 veces, se optimizaron y se probaron sobre el conjunto test.

Resultados: Los algoritmos de regresión del proceso gaussiano proporcionaron los mejores modelos con un error absoluto medio fue de alrededor de 1.00 D en el componente esférico y de 0.25 D en los componentes astigmáticos.

Conclusiones: Los pesos normalizados mostraron que la refracción subjetiva puede predecirse utilizando únicamente la queratometría, la edad y la longitud axial como características. El aumento del detalle de la descripción topográfica de la superficie corneal anterior que supone una descomposición de Zernike de alto orden frente al ajuste a una superficie esferocilíndrica realizado por queratometría no se refleja en una mejora de la predicción de la refracción subjetiva, que es pobre, en cualquier caso, principalmente en el componente esférico. Sin embargo, la máxima precisión alcanzada difiere en sólo 0,75 D de la de otros trabajos con una descripción más exhaustiva de los elementos refractivos del ojo. De todos modos, el aprendizaje automático y los datos masivos aplicados a la predicción de la refracción subjetiva pueden ser arriesgados y poco prácticos cuando se evalúa a un sujeto concreto en los extremos estadísticos, aunque los parámetros elegidos puedan no haber sido los más ineficaces.

Relevancia Traslativa: El trabajo evalúa la predicción de la refracción subjetiva mediante aprendizaje automático a partir de la superficie corneal anterior y la biometría ocular, mostrando la mínima información biométrica requerida y la máxima precisión alcanzable.

Introduction

Understanding how the human visual system works has become a recursive research task in recent years. Since the beginning of the past century, several optical models of the human eye have been developed to predict visual quality or subjective refraction at different complexity levels, and both comprise real measurements (in vivo and in vitro) and theoretical approximations.^1,2 The two fundamental elements are the cornea and crystalline lens. Corneal topography can be assessed by different techniques like interferometry, ultrasonography, profile photography, holography, Placido disk principle, and Scheimpflug photography. Corneal surfaces are usually described with the keratometry³ and/or Zernike coefficients,⁴ c_j, which result from fitting a series of Zernike polynomials to heights maps. The Zernike decomposition of a surface W can be expressed as polar coordinates⁵:  

\begin{equation}W\left( {\rho ,\theta } \right) = \mathop \sum \limits_{j = 0}^{p - 1} {c_j}{Z_j}\left( {\rho ,\theta } \right);\end{equation}

(1)

 

\begin{eqnarray}j &=& 0.5\left( {n\left( {n + 2} \right) + m} \right);\nonumber\\ n &=& \left\lceil {0.5\left( { - 3 + \sqrt {9 + 8j} } \right)} \right\rceil ;\nonumber\\ m &=& 2j - n\left( {n + 2} \right);\end{eqnarray}

(2)

 

\begin{equation}{Z_{n,m}}\left( {\rho ,\theta } \right) = {N_{n,m}}{R_{n,m}}\left( \rho \right){M_m}\left( \theta \right)\end{equation}

(3)

 

\begin{eqnarray}{N_{n,m}} &=& \sqrt {\frac{{2\left( {n + 1} \right)}}{{1 + {\delta _{m0}}}}} ;\nonumber\\ {\rm{\;\;\;}}{R_{n,m}}\left( \rho \right) &=& \mathop \sum \limits_{s = 0}^{0.5\left( {n - \left| m \right|} \right)} \frac{{{{\left( { - 1} \right)}^s}\left( {n - s} \right)!}}{{s!\left[ {0.5\left( {n + \left| m \right|} \right) - s} \right]!\left[ {0.5\left( {n - \left| m \right|} \right) - s} \right]!}}{\rho ^{n - 2s}};\end{eqnarray}

(4.a)

 

\begin{equation}{M_m}\left( \theta \right) = \left\{ \begin{array}{@{}*{1}{c}@{}} {\cos \left( {m\theta } \right)\;\;\;\;\;\;\;\;\;m \ge 0}\\ {\sin \left( {m\theta } \right)\;\;\;\;\;\;\;\;\;m < 0} \end{array}\right.;\end{equation}

(4.b)

where n is the radial order, m is the azimuthal frequency, j is the single index for the Zernike polynomial, p is the number of terms in the expansion, c_j are the Zernike coefficients associated with their Zernike polynomial, δ_m0 is the Kronecker delta function, and ⌈ · ⌉ denotes the ceiling (round-up) operator. Regarding the crystalline lens, on the one hand, phakometry,⁶ Scheimpflug imaging,^7–10 magnetic resonance imaging,¹¹ and optical coherence tomography^12,13 have been used to assess the in vivo optical properties of lens shape and lens thickness. On the other hand, in vitro techniques^14–18 have also been followed to evaluate lens shape and power. 

Subjective refraction is responsible for determining the most adequate optical power needed to compensate a patient’s refractive errors. It relies on the patient's ability to discern and communicate possible improvements or distortions caused by corrective adjustments made by the optometrist, and thus, its determination with automated machines is still not accurate. This procedure takes time and, consequently, slows down the diagnosis process. To reduce the time needed for this process, machine learning has been recently used to predict refractive prescription from physical eye data obtained by wavefront aberrometers,^19,20 photorefraction images,²¹ retinal fundus images,²² other ophthalmologic devices,²³ and intraocular lenses characteristics.²⁴ The aim of this work is to develop a machine learning regression model that predicts patients’ subjective refractive prescription from the anterior corneal surface and ocular biometry. The first step consists of choosing physiologic descriptors as the model's features. It is a key point because using too many features can degrade prediction performance, even if all features are relevant.²⁵ The strong correlation between the spherical power components and astigmatic components²⁶ of the anterior and posterior corneal surfaces, as well as the fact that the latter contributes only about one-eighth of the eye's refractive power,²⁷ leads to the hypothesis that the anterior surface suffices to assess the whole corneal refractive effect on an eye model. Apart from anterior corneal topography, other physiologic parameters that are expected to be involved in patients’ subjective refraction were added to the model: axial length (AL) and then constructing a similar approach to Emsley schematic eye; the patient's age, on which refractive indices and crystalline lens morphology depend; and, finally, anterior chamber depth (ACD), which can be related to lens location. Compared to previous works,²⁰ this proposal offers the advantages of using simpler measurement devices, requiring fewer descriptive data and, therefore, faster performance. However, disregarding the crystalline lens effect will probably provide poor results. In fact, the distribution of some aberrations between the cornea and lens appears to be auto-compensated.²⁸ Next, the selection and tuning of machine learning models must be performed with training and test population subsets. The Results section shows an analysis of the main indispensable characteristics and the evaluation of the final predictive accuracy of the selected models. 

Methods

The electronic medical records of 229 patients were retrieved. This study followed the Declaration of Helsinki principles, and the participants gave their written informed consent. Any patients with any ocular abnormalities, except ametropia and ever-performed cataract surgery or artificial lens implants, were excluded. The resulting selected measurements consisted of the anterior corneal surface parameters, AL, ACD, age, and noncycloplegic subjective refractions conducted by optometrists with 201 patients (355 eyes), 154 patients’ bilateral eye records and 47 patients’ unilateral eye records. Subjective refraction (S, C, α) was transformed through equation (5) to power vector notation²⁹ (M, J₀, J₄₅), which consists of components that are independent of one another. Figure 1 shows the histograms of the age and power vector components of the population sample herein used.  

\begin{eqnarray}M &=& S + \frac{C}{2};\;\;\;\;\;\;\;\;\;{J_0} = - \frac{C}{2}\cos \left( {2\alpha } \right);\nonumber\\ {J_{45}} &=& - \frac{C}{2}\sin \left( {2\alpha } \right)\end{eqnarray}

(5)

 

Figure 1.

 

Histograms of the (a) age and (b–d) power vector components (M, J₀, J₄₅) obtained from the subjective refraction performed with the sample population.

Both ACD and AL were measured by an optical biometer IOLMaster 500 (Carl Zeiss Meditec, Inc., Dublin, CA, USA). Corneal surface measurements were taken with a Sirius Topographer (Sirius, CSO, Firenze, Italy). The Sirius is a Scheimpflug combined with a Placido disk imaging system that allows measurements of 256 meridians of corneas in 30 radial distances to be taken. It provides several corneal physiologic descriptors susceptible to be selected as features in machine learning (Zernike coefficients from both surfaces with different diameters, keratometry, etc.). The Zernike coefficients describing the anterior corneal surface at a pupil diameter of 4.5 mm and anterior corneal surface astigmatism (C_k,α_k) were selected from the measurements taken by the Sirius Topographer as corneal physiologic descriptors. As the device does not provide defocus coefficient c₄, we used the keratometry equivalent sphere (M_k) instead, which was computed from the keratometry for the same pupil diameter as  

\begin{equation}{M_k} = 0.5\left( {{n_k} - 1} \right)\left( {\frac{1}{{{R_f}}} + \frac{1}{{{R_s}}}} \right)\end{equation}

(6)

where n_k is the keratometric index and R_f and R_s are, respectively, the flattest and the steepest anterior corneal curvature radius. The conventional keratometric index (1.3375) was used in this work. However, its value is not relevant for machine learning algorithms because it is constant. Keratometry data (C_k,α_k) were also transformed into standard power vector notation²⁹ according to equation (5). 

The database was split into a training subset and a test subset with 75% and 25% of the patients, respectively, to avoid overfitting. It comprised the ocular measurements of both unilateral and bilateral patients. Splitting was performed by preventing the eyes of the same patient from appearing in both subsets. The training subset consisted of 36 unilateral eyes, plus both eyes of 115 bilateral patients. The test subset included 11 unilateral eyes, plus both eyes of 39 bilateral patients. Table 1 presents the descriptive statistics that summarize the refractive data and ages of the different sets. 

Table 1.

 

Descriptive Statistics Summarizing the Subjective Refractive Data and Ages of the Different Sets

 

Data set size determines the maximum number of features that can be postulated as significant characteristics in the machine learning algorithm. In this work, the training data included 266 eyes. However, a 10-fold cross-validation technique was applied to the training data to assess how the results of a statistical analysis would generalize to an independent data set. This process consists of randomly partitioning the training sample is into 10 equally sized subsamples. Of these 10 sets, a single subsample is retained as the validation data to test the model. The remaining nine sets are used as data to train the model. Therefore, the effective size of the training data is 266 × 9 / 10 ≈ 239. A rule of thumb establishes that the maximum number of parameters for a good-performing model is limited to one-tenth of the amount of training data. This left the maximum number of features at 23. 

Sirius provides Zernike coefficients up to the seventh order except piston (c₀) and defocus (c₄), that is, 34 coefficients from the anterior corneal surface decomposition that are above the maximum limit of 23 features. Bearing in mind that a patient’s age, M_k, ACD, and AL are all features that characterize the eye apart from cornea, a thus-detailed anterior corneal surface description can be foregone to meet that limit. Therefore, those coefficients, and excluding tilts (c₁ and c₂) for their little relevance as they can be naturally compensated by eye movements, were selected in ascending order without exceeding that limit. The fifth order includes 17 Zernike coefficients that, together with age, M_k, ACD, and AL, make up a set of 21 selected features as a first approach, to characterize each eye and to train the models. As part of the preprocessing data step, a filter-type feature selection algorithm³⁰ that used a diagonal adaptation of the neighborhood component analysis³¹ was applied to determine the features’ normalized weights (NWs). This reports about the importance of each feature in the regression models of the power vector components and allows different feature selections to train and tune the new models to be run. An extra parameter with a random value (rv) was added to check the significance of the 21 selected features in the model. The NW of the extra random parameter is expected to be zero. If not, any feature with an NW that equals or goes below that can be discarded because its significance would be equal or worse than that of a random variable. 

Regression models were trained with the regression learner app³² of MATLAB (version R2021b; MathWorks, Inc., Natick, MA, USA), which provides linear regression models,³³ regression trees,³⁴ support vector machines (SVMs),³⁵ Gaussian process regression (GPR) models,³⁶ ensembles of trees,^37–39 and neural networks.⁴⁰ Algorithms were trained to predict each refraction vector component separately because of their mathematical independence. The models that obtained the best root mean squared errors (RMSEs) were later hyperparameter-tuned by Bayesian optimization. Finally, the best models were tested with the test subset data. Otherwise, a model might only perform well with the training data but may fail to predict anything useful in yet unseen data. 

Results

The neighborhood component analysis feature selection with regularization was applied to determine feature importance to predict the power vector components of the subjective refraction. Figure 2 is a bar graph showing the computed NWs of all 21 features, plus the random one to predict the three power vector components. 

Figure 2.

 

Normalized weights of the features for predicting the three power vector components. c_j stands for the jth Zernike coefficient in line with (1).

The next step consisted of determining the machine learning algorithm that worked the best. Different regression models were trained using the regression learner app³² with the 266 eyes from the training subset characterized by the 21 features and by following the above-described procedure. Table 2 shows the RMSE and the coefficient of determination (R²) values obtained for each power vector component with the trained models, which gave better results. R² compares the trained model to the model with a constant response, and it equalled the training response mean. If the model is worse than this constant model, then R² is negative and the model is discarded. 

Table 2.

 

RMSE and Coefficient of Determination (R²) Values Obtained for Each Power Vector for the Training Subset With All the Trained Models

 

According to the RMSE, the GPR models provided the best results in the first approaches for the three power vector components, although other models approximately matched those values for some components. Therefore, the GPR, Ensembles, and SVM models were subsequently optimized by hyperparameter tuning to obtain the best results in the GPR models. The results appear in the last row of Table 2. In Figures 3a and 3b, the responses obtained with the optimized GPR models for the power vector components of the training set are represented versus the true data. 

Figure 3.

 

Predicted responses obtained with the optimized GPR models versus the true data of the training set for the three power vector components.

The performance of the optimized models in the test subset (89 eyes) was also evaluated. Figures 4a and 4b show the predicted responses versus the true response for the three power vector components. Table 3 includes the RMSE, mean absolute error (MAE), and R² obtained from the test subset for each power vector through the optimized GPR model trained for 21 features. 

Figure 4.

 

Predicted responses obtained with the optimized GPR models versus the true data of the test set for the three power vector components.

Table 3.

 

RMSE, MAE, and Coefficient of Determination (R²) Obtained for Each Power Vector for the Test Subset With the Optimized GPR Trained Model and the Proposed Features

 

Discussion and Conclusions

Machine learning regression algorithms were been trained to predict the subjective refraction prescription from the anterior corneal shape parameters, eye lengths, and the patients’ age. The Gaussian process regression algorithms provided the best models with different accuracies for each subjective refraction power vector component. Algorithms were applied to a small data set compared to other studies found in the bibliography. However, data adhered to the rule of 10 in the worst case, namely, the amount of training data needed for a good-performing model being 10 times the number of parameters in the model. Therefore, the performance of the obtained models was reliable. 

A selection of the features with weight that appears in Figure 2 would provide equal or better results when performing the regression models for each power vector component than those from the 21 features. The effect of neglecting some of those features can be assessed by establishing different weight thresholds to make the selections. Two different thresholds were set. The first comprised features with NWs above zero. The second contained those with NWs above 20%, which depreciated the ACD effect, with higher-order aberrations c₆ (vertical trefoil), c₇ (vertical coma), c₁₂ (primary spherical), c₁₃ (vertical secondary astigmatism), and c₁₉ (oblique secondary trefoil). Hence, the M component features were consecutively in NW importance terms AL, the keratometry equivalent sphere, and age. Both J₀ and J₄₅ were respectively modeled using the primary astigmatism Zernike coefficient (c₃ and c₅). The independency of the three components was evidenced because they showed different NWs of features. Hence, a different subset of features could be selected for each one to perform the regression models. However, common selections of features have been sought for three components for simplicity's sake and to obtain only a set of features for each threshold. Table 4 shows the results, for the test subset, of the trained and optimized GPR models by considering features according to these two different thresholds for the three power vector components. 

Table 4.

 

RMSE, MAE, and Coefficient of Determination (R²) Obtained for Each Power Vector for the Test Subset With the Optimized GPR Trained Model With Different Eye Features

 

The fact that the regression performed with the features of NW >0 slightly improved RMSE is the first point to be noted. As expected, using too many features degrades prediction performance. Indeed, the R² of component J₄₅ improved. The results of the optimized GPR model with all features and of that with those of NW >20% (five features) were similar. Consequently, a regression model with only these five parameters can achieve maximum accuracy. Disregarding high-order aberrations and ACD did not significantly make the results worse. 

Based on the above result, good results can be expected by considering age, AL, and keratometry, which describes the corneal defocus and astigmatism. Table 4 also present the obtained errors and R² from the test subset for the models trained with keratometry, age, and AL. One again, the results were not different from those obtained using the models with all the features. 

The analysis of the computed normalized weights of features allowed the number of required parameters to be minimized. This work revealed the highest accuracy that can be obtained using the available information and confirmed some preestablished facts. The models’ predictive ability for components J₀ and J₄₅ was good (around 0.25 diopters [D] of RMSE), but the coefficient of determination for J₄₅ was poor. For spherical component M, although R² was good, the error (MAE around 1.00 D) seemed to be too high to employ it as a technique to predict subjective refraction. Notwithstanding, the achieved accuracies were not much worse than those reported in previous works,²⁰ although the presented proposal used simpler measurement devices (a keratometer would work) and required fewer descriptive data. 

Our results are probably not as accurate as the models contained no data about posterior corneal surface and crystalline lens aberrations. Machine learning algorithms employed age, ACD, and AL to approach the eye's inner optical part effect. On the one hand, it proved that the geometric parameters of the human crystalline lens and refractive indices were age dependent and, therefore, so was lens refractive power. The analysis of the NW of the features showed that the NW of the age was above 20% when modeling the M component. This fact confirmed the marked dependence of the refractive spherical component on this feature. The refractive power of the missing elements in the eye model was slightly described by the age feature. A patient’s age also appeared in the NW features analysis for predicting J₀ and J₄₅. This could indicate the age dependence of astigmatism.^42,43 The ACD effect, which was hypothesized to be related to the lens location in the eye, and the higher-order Zernike coefficients effect were assessed from the results obtained with the parameters of NW >0 in the M component model. The result obtained considering these effects was not significantly better (only about 0.1 D of MAE) than those obtained in the model with the features of NW >20%, which excluded these parameters in relation to the previous ones. The comparison made of the astigmatic power vector components showed that errors were no higher when those coefficients were not considered. Hence, we conclude that their influence was slight and, therefore, machine learning algorithms can dispense these features without making the predictions worse. 

On the other hand, as AL determines the location of the retina, it is fundamental to model the spherical component. Indeed, as this feature had the highest NW for the M component prediction, not taking it into account would provide poor results. It was also present for the J₀ component prediction but with an NW below 20% and, hence, its poor significance. 

Astigmatic components are expected to be solely determined by the anterior corneal surface shape. However, as we have just cited, the NW analysis of features to predict J₄₅ and J₀ showed other features apart from c₃ and c₅ (oblique and vertical astigmatism Zernike coefficients, respectively), albeit with lower NWs. The models that optimized GPR and only employed the astigmatic components were also trained and tested. The results in Table 5 are practically no different from those obtained using all the parameters with NWs other than zero, which implies that anterior corneal surface topography determines astigmatism. These results, together with the fact that the errors of models for astigmatic components were lower than the error of the spherical ones, confirmed the strong correlation between the astigmatic components²⁶ of the anterior and the posterior corneal surfaces. 

Table 5.

 

RMSE, MAE, and Coefficient of Determination (R²) Obtained for the Astigmatic Power Vector for the Test Subset With the Optimized GPR Trained Model With c₃ and c₅ Coefficients

 

The weight feature analysis also allowed us to conclude that no precise corneal surface description is needed for subjective refractive prescriptions. The errors of the astigmatic components for all the evaluated models remained unchanged, regardless of whether keratometry or a high-order anterior corneal Zernike decomposition was considered or not. The obtained MAE values were of the same order as those obtained by Rampat et al.,²⁰ who contemplated whole-eye aberrations. Therefore, the contribution of the posterior corneal surface and the crystalline lens to the eye's astigmatism falls within those errors. For the spherical component, the difference with that work²⁰ was below 0.75 D. Therefore, by assuming the strong correlation between the spherical power²⁶ of corneal surfaces, disregarding the lens would be partially, but not completely, compensated by considering age. The difference in patients’ demographics between both works lay in them having a refractive spherical equivalent between −6.75 D and 6.13 D, which, in this work, is between −11 D and 6.25 D. This could also be a cause of the worse results herein obtained for the widest range of ametropies. 

Machine learning allowed us to confirm that the low-order aberrations obtained from anterior corneal surfaces, together with age and AL, suffice to approximately predict subjective refraction. This is an advantage in diagnosis speed terms over the subjective traditional technique. The roles of the proposed features in modeling each power vector component were assessed. 

The traditional subjective refraction is nowadays still considered the universal gold standard in the evaluation of refractive error. However, it is time-consuming, taking about 5 minutes per eye⁴⁴ by a well-trained eye care professional, and its subjective nature implies deviations in the spherical component around ±0.25 D^45,46 and in the astigmatic component above 0.75 D⁴⁷ in both intra- and interoptometrist variability. 

The aim of this work was to develop of a machine learning regression model that predicts patients’ subjective refractive prescription from the anterior corneal surface and ocular biometry. It is an objective technique that can be used to measure the spherical and cylindrical refractive errors of the human eye. Objective refraction is not only useful but often essential, for example, when examining young children and patients with mental or language difficulties. Cycloplegic refraction, for example, requires little or no cooperation from the patient. This is an advantage over other techniques, such as the one proposed in this study, which requires biometry and corneal surface measurements. However, the administration of cycloplegics may cause undesirable ocular and/or systemic side effects. Moreover, although objective refraction can provide good visual outcomes, neural processing must also be considered. 

Our proposal, compared to the work of Rampat et al.,²⁰ offers the advantages of using simpler measurement devices, requiring fewer descriptive data, and, therefore, having faster performance. It fits in with the unsupervised methods’ trend of minimizing misunderstandings between the clinician and the patient and the measurement variability and time. However, we must consider that the proposed approaches statistically predict subjective refraction within a tolerance range for a population, and, individually, the errors that they may make can be intolerable. So although the choice of input parameters in this work might not be the most predictive efficient, a general model should clearly show its extreme errors and not merely its average performance. Otherwise, establishing a standard predictive method based on machine learning algorithms and big data can be risky and impractical when evaluating a particular patient. 

Acknowledgments

The authors thank Vissum Grupo Miranza (Alicante) for providing the electronic medical records. 

Disclosure: J. Espinosa, None; J. Pérez, None; A. Villanueva, None