October 2022
Volume 11, Issue 10
Open Access
Artificial Intelligence  |   October 2022
Predictive Modeling of Long-Term Glaucoma Progression Based on Initial Ophthalmic Data and Optic Nerve Head Characteristics
Author Affiliations & Notes
  • Eun Ji Lee
    Department of Ophthalmology, Seoul National University College of Medicine, Seoul National University Bundang Hospital, Seongnam, Korea
  • Tae-Woo Kim
    Department of Ophthalmology, Seoul National University College of Medicine, Seoul National University Bundang Hospital, Seongnam, Korea
  • Jeong-Ah Kim
    Department of Ophthalmology, Kangwon National University School of Medicine, Chuncheon, South Korea
  • Seung Hyen Lee
    Department of Ophthalmology, Nowon Eulji Medical Center, Eulji University College of Medicine, Seoul, Korea
  • Hyunjoong Kim
    Department of Applied Statistics, Yonsei University, Seoul, Korea
  • Correspondence: Tae-Woo Kim, Department of Ophthalmology, Seoul National University Bundang Hospital, Seoul National University College of Medicine, 82, Gumi-ro, 173 Beon-gil, Bundang-gu, Seongnam, Gyeonggi-do 13620, Republic of Korea. e-mail: [email protected] 
Translational Vision Science & Technology October 2022, Vol.11, 24. doi:https://doi.org/10.1167/tvst.11.10.24
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to authenticated users only.
      Sign In or Create an Account ×
    • Get Citation

      Eun Ji Lee, Tae-Woo Kim, Jeong-Ah Kim, Seung Hyen Lee, Hyunjoong Kim; Predictive Modeling of Long-Term Glaucoma Progression Based on Initial Ophthalmic Data and Optic Nerve Head Characteristics. Trans. Vis. Sci. Tech. 2022;11(10):24. https://doi.org/10.1167/tvst.11.10.24.

      Download citation file:


      © ARVO (1962-2015); The Authors (2016-present)

      ×
  • Supplements
Abstract

Purpose: The purpose of this study was to develop a model, based on initial optic nerve head (ONH) characteristics, predictive of long-term rapid retinal nerve fiber layer (RNFL) thinning in patients with open-angle glaucoma (OAG).

Methods: This study evaluated 712 eyes with OAG that had been followed up for >5 years with annual evaluation of RNFL thickness. Baseline ophthalmic features were incorporated into the machine learning models for prediction of faster RNFL thinning. The model was trained and tested using a random forest (RF) method, and was interpreted using Shapley additive explanations. Factors associated with faster rate of RNFL thinning were statistically evaluated using a decision tree.

Results: The RF model showed that greater lamina cribrosa (LC) curvature, higher intraocular pressure (IOP), visual field mean deviation converging towards −5 dB, and thinner peripapillary choroid at baseline were the four most significant features predicting faster RNFL thinning. Partial interaction between the features showed that larger LC curvature was a strong factor for faster RNFL thinning when it exceeded approximately 12.0. When the LC curvature was ≤12, higher initial IOP and thinner peripapillary choroid played a role in the rapid RNFL thinning. Based on the decision tree, higher IOP (>26.5 mm Hg), greater laminar curvature (>13.95), and thinner peripapillary choroid (≤117.5 µm) were the 3 most important determinants affecting the rate of RNFL thinning.

Conclusions: Baseline ophthalmic data and ONH characteristics of patients with OAG were predictive of eyes at risk of faster progression. Combinations of important characteristics, such as IOP, LC curvature, and choroidal thickness, could stratify eyes into groups with different rates of RNFL thinning.

Translational Relevance: This work lays the foundations for developing prediction models to estimate glaucoma prognosis based on initial ONH characteristics.

Introduction
Glaucoma is a progressive optic neuropathy, with phenotypes that vary widely among patients. Some patients rapidly progress to significant visual impairment despite extensive treatment, whereas others do not progress or progress very slowly even without treatment. 
The pathogenic mechanism underlying glaucoma is not fully understood. Because lowering intraocular pressure (IOP) is the only proven way to slow the progression of glaucoma, current treatments mainly involve reducing the IOP. The IOP alone, however, cannot explain the mechanism of glaucomatous damage. Thus, patients may experience disease progression despite treatments that reduce IOP. 
Eyes that show rapid progression of glaucoma despite treatment require treatment adjustment, either intensification of IOP lowering medication or a switch to surgical methods. This treatment would be excessive, however, in patients with stable disease and may negatively affect their quality of life. The ability to predict glaucoma progression at baseline in individual patients may allow the determination of treatment goals, optimizing both treatment and the patient’s quality of life. In addition, identification of individual factors predictive of progressive disease may help understand the mechanism of glaucomatous optic neuropathy for individual eyes, enabling tailored treatment of each eye. Even if a practical method to reduce the rate of progression is unavailable, information on disease prognosis may help in advising individual patients about future progression and enable them to seek other strategies. 
Disease prediction models stratify patients based on their probability of achieving certain outcomes, thereby identifying patients at increased risk of an event. Using various methodologies, efforts have been made to develop glaucoma prediction models, with these models showing promising results in their ability to predict the occurrence of glaucoma1,2 and the need for glaucoma surgery.3 Relatively little is known, however, about the association of individual ophthalmic features with faster or slower glaucoma progression. In addition, despite increasing evidence about the pathogenic importance of peripapillary structures, such as the lamina cribrosa (LC) and choroid, in glaucomatous optic neuropathy, the influence of the characteristics of these peripapillary structures has not been investigated.49 
We hypothesized that machine learning models trained with baseline ophthalmic data may be prognostic in identifying patients at high risk of rapid glaucoma progression. The purpose of this study was to establish a prediction model based on initial ophthalmic characteristics for rapid, long-term glaucoma progression. Baseline ophthalmic characteristics, including optic nerve head (ONH) characteristics, were incorporated into machine learning models to develop a model predictive of rapid progressive retinal nerve fiber layer (RNFL) thinning over 5 years. 
Methods
The study included patients with primary open angle glaucoma (POAG) who had been enrolled in the Investigating Glaucoma Progression Study, an ongoing prospective study of patients with glaucoma at the Glaucoma Clinic of Seoul National University Bundang Hospital. All subjects provided written informed consent to participate. The study protocol was approved by the Institutional Review Board of Seoul National University Bundang Hospital and adhered to the Declaration of Helsinki. 
Inclusion and Exclusion Criteria
Patients were included if they had been diagnosed with POAG between September 17, 2009, and February 11, 2014, were followed up in the glaucoma clinic for >5 years, and had undergone annual spectral-domain optical coherence tomography (OCT) examinations to measure circumpapillary RNFL thickness (RNFLT). Only eyes with a global average RNFLT ≥50 µm at baseline were included in the present study, to avoid the floor effects of OCT measurements that can result in a falsely slow rate of RNFL thinning in eyes with extremely thin baseline RNFL.10,11 
POAG was diagnosed based on open iridocorneal angle on gonioscopy and signs of glaucomatous optic nerve damage (i.e. neuroretinal rim thinning, notching, or an RNFL defect) with consistent visual field (VF) defect. Glaucomatous VF defect was defined by one or more of the Anderson-Patella criteria: (1) glaucoma hemifield test results outside normal limits, (2) clusters of ≥3 non-edge points on a pattern deviation plot with probabilities <5%, with one having a probability <1%, and (3) <5% probability of pattern standard deviation confirmed on two consecutive reliable tests (fixation loss rate ≤20% and false-positive and false-negative error rates ≤25%). 
Eyes were excluded if they had a best-corrected visual acuity worse than 20/40; a spherical equivalent of <–8.0 diopters (D) or >+3.0 D; a cylinder correction of <–3.0 D or >+3.0 D; a history of intraocular surgery, except for uneventful cataract surgery; or a history of retinal (e.g. diabetic retinopathy, retinal vessel occlusion, or retinoschisis) or neurological (e.g. pituitary tumor) disease. When both eyes were eligible, one eye was randomly selected for this study. 
Baseline Ophthalmic Examination
All participants underwent comprehensive ophthalmic examinations, including assessments of best corrected visual acuity, Goldmann applanation tonometry, slit-lamp biomicroscopy, gonioscopy, dilated stereoscopic examination of the optic disc, disc photography, red-free fundus photography (EOS D60 digital camera; Canon, Utsunomiyashi, Tochigiken, Japan), spectral-domain OCT scanning of the circumpapillary RNFL and the ONH (Spectralis; Heidelberg Engineering, Heidelberg, Germany), and standard automated perimetry (Humphrey Field Analyzer II 750; 24-2 Swedish interactive threshold algorithm, Carl Zeiss Meditec). Other ophthalmic examinations included measurements of corneal curvature (KR-1800; Topcon, Tokyo, Japan), central corneal thickness (CCT; Orbscan II; Bausch & Lomb Surgical, Rochester, NY, USA), and axial length (AXL; IOLMaster version 5; Carl Zeiss Meditec, Dublin, CA, USA). 
Measurements of the Curvature of the Lamina Cribrosa and Peripapillary Choroidal Thickness
Features of the ONH were evaluated in two ways: by measuring the LC curvature and peripapillary choroidal thickness (CT), using the ONH volume images obtained by enhanced depth-imaging (EDI) OCT, and using the peripapillary circular scanning centered on the Bruch's membrane opening, respectively. The corneal curvature of each eye was entered into the Spectralis OCT system before performing the OCT scanning to compensate for potential magnification errors. 
The LC curvature is an index of the degree of posterior bowing or deformation of the LC, which is considered to represent the IOP-induced mechanical stress on the ONH.8,12,13 The degree of LC deformation has been shown to decrease after reducing IOP,1416 with a larger degree of LC deformation being associated with rapid glaucoma progression.7,8,17 
The LC curvature was evaluated using the LC curve index (LCCI), defined as the degree of inflection of the curve representing a section of the LC.13 The LCCI was determined by measuring the width of the anterior LC surface reference plane (W) and then by measuring the LC curve depth (LCCD) within the anterior LC surface plane in each B-scan, with LCCI calculated as (LCCD/W) × 100. Because the curvature was normalized to the LC width, the LCCI is an indicator of the shape of the LC independent of the actual size of the ONH. Only the LC within the W was considered because the LC was often not clearly visible outside this region. In eyes with LC defects, the LCCI was measured using a presumed anterior LC surface that best fit the curvature of the remaining part of the LC or by excluding the area of the LC defect. A manual caliper tool in the Heidelberg Eye Explorer (version 1.10.4.; Heidelberg, Germany) was used to measure the W and LCCD in seven selected B-scan images spaced equidistantly across the vertical optic disc diameter in each eye, with the mean LCCI of each eye calculated using the measurements made from these seven B-scans. 
The peripapillary choroid is considered an important vascular structure that perfuses the ONH. In addition to the choroidal and ONH vasculatures having the same origin (i.e. the short posterior ciliary artery), the peripapillary choroid has been shown to directly perfuse ONH tissue.5,18 A thinner peripapillary choroid has been regarded as indicative of decreased perfusion of the ONH.1921 
Peripapillary CT was measured on the 360 degree 3.5-mm diameter peripapillary circle scan centered on the Bruch's membrane opening (BMO) using the manual segmentation function built into the Heidelberg Eye Explorer software. The posterior edge of the retinal pigment epithelium and the sclerochoroidal interface, representing the inner and outer boundaries of the choroid, respectively, were manually delineated. Using the circumpapillary RNFLT measurement algorithm, the peripapillary CT in the global area and in the six sectors based on the foveal–BMO axis were automatically generated. 
LCCI and peripapillary CT were measured by two experienced observers (authors J.A.K. and S.H.L.) who were masked to the clinical information of participants. The measurements by the two observers were averaged for analyses. The interobserver agreements for measuring the LCCI and peripapillary CT were assessed by calculation of intraclass correlation coefficients (ICCs) and 95% confidence intervals (CIs). 
Input Variables
Data from the baseline examinations were incorporated into a machine learning modeling. Eleven features were selected as input variables: age, sex, highest IOP during the initial 6 months, glaucoma surgery during the initial 6 months, mean LCCI, global peripapillary CT, global RNFLT, VF mean deviation (MD), VF pattern standard deviation (PSD), AXL, and CCT. 
Not all participants had measurements of pretreatment IOP, thus the highest IOP recorded during the initial 6 months (highest IOP during the initial 6 months) was defined as the initial IOP, and used for the analysis. To exclude any potential effect of IOP lowering therapy on the morphology of the LC1416 and the CT,22,23 LCCI and CT were measured on OCT scans obtained >6 months after starting IOP-lowering treatment or 6 months after glaucoma surgery if surgery was performed during the initial 6 months. 
Output Variable
The rate of glaucoma progression was defined as the output variable and was assessed by calculating the rate of change of RNFL thickness relative to the baseline thickness per year, expressed as percent per year. The change of global average RNFL thickness over time was subjected to linear regression analysis to determine the rate of change for each subject. 
Machine Learning Prediction Models and Training Protocols
The machine learning model used for fitting and prediction was a random forest (RF) model.24 Training data were used to train the model, and test data for prediction. When fitting the training data, hyperparameters were tuned using the five-fold cross-validation method to prevent overfitting of the RF model. Because machine learning models are difficult to interpret due to their complexity, the results of the RF model in this study were interpreted using the Shapley additive explanations (SHAP) method,25 a type of explainable artificial intelligence method. 
Training and Test Sets
The data were randomly separated so that both eyes of the subjects belonged to only one of the test or training sets. The number of observations in the training and test data were 356 each. 
Random Forest
RF is an ensemble learning method for classification and regression that operates by constructing multiple decision trees in the learning phase. Although RF has been shown to outperform the predictive performance of other machine learning methods, it is difficult to interpret the predicted results of the RF method, a disadvantage similar to that of other machine learning methods. The number of decision trees for RF learning was set at 200, with 5-fold cross-validation showing that the optimal number of features for each node was 3. 
Shapley Additive Explanation
The SHAP method was developed to explain the output of any machine learning model. After constructing a model with several features, the SHAP value was obtained by determining the average change relative to the presence or absence of any individual feature. The SHAP value for each feature was an indicator of how strongly that feature affects the prediction of the model in a positive or negative direction, with a larger absolute SHAP value indicating that the feature had a greater impact on prediction by the model. Feature importance and partial dependence plots were obtained using the SHAP values. Interactions between features were visualized by partial interaction plots. 
Regression and Decision Tree Model
Univariate and multivariable regression analyses were performed to identify the factors associated with the rate of RNFL thinning, with the multivariable model including variables with P < 0.10 on univariate analyses. A decision tree model was used to stratify patients with faster or slower OCT RNFL thinning, based on factors influencing the rate of RNFL thinning. All possible splits of this type were considered, with the cutoff that best separated the data into groups being chosen. The end points of this tree were the rate of OCT RNFL thinning (micrometers per year). Factors found to be associated with the rate of OCT RNFL thinning on multivariable regression analyses were regarded as possible risk factors in the decision tree model. Insignificant splits were removed by 10-fold cross validation. 
Data Analysis
Except where stated otherwise, data are presented as mean ± standard deviation. All statistical and machine learning analyses were performed using Python version 3.8.5, with the scikit-learn package version 1.0 used for regression analysis and decision tree construction and the SHAP package version 0.39.0 used for SHAP analysis. Linear regression analyses were performed using Statistical Package for the Social Sciences software (version 22.0; SPSS, Chicago, IL), with P values < 0.05 considered statistically significant. Model performance was assessed by calculating mean absolute error (MAE) and mean squared error (MSE). 
Results
Table 1 shows the baseline characteristics of the 712 participants. Data from 357 and 355 eyes were designated as the training and test sets, respectively, to build the final model. The MAEs for the RF, regression, and decision tree models were 0.075, 0.115, and 0.128, respectively, and the MSEs were 0.023, 0.031, and 0.076, respectively. There were excellent interobserver agreements in measurements of LCCI and peripapillary CT, with ICCs of 0.996 (95% CI = 0.994–0.998), and 0.960 (95% CI = 0.940–0.974), respectively. 
Table 1.
 
Baseline Characteristics or Participants (n = 712)
Table 1.
 
Baseline Characteristics or Participants (n = 712)
SHAP was used to determine the features that provide high predictive power. The four features with the highest mean absolute SHAP values were mean LCCI, initial IOP, VF MD, and peripapillary CT (Fig. 1A). The SHAP plot revealed that larger mean LCCI, higher initial IOP, VF damage of moderate degree (VF MD close to −3 dB), and thinner global CT were associated with faster RNFL thinning (Fig. 1B). 
Figure 1.
 
Interpretation of the final model based on baseline patient variables. (A) Feature importance plot based on mean SHAP values. (B) Interpretation of the importance of features using the SHAP plot. The red and blue colors indicate feature values of high and low levels, respectively. For example, a larger mean LCCI had a strong and negative impact on the rate of RNFL thinning (i.e. faster RNFL thinning). LCCI, lamina cribrosa curve index; IOP, intraocular pressure; VF, visual field; MD, mean deviation; CT, choroidal thickness; AXL, axial length; PSD, pattern standard deviation.
Figure 1.
 
Interpretation of the final model based on baseline patient variables. (A) Feature importance plot based on mean SHAP values. (B) Interpretation of the importance of features using the SHAP plot. The red and blue colors indicate feature values of high and low levels, respectively. For example, a larger mean LCCI had a strong and negative impact on the rate of RNFL thinning (i.e. faster RNFL thinning). LCCI, lamina cribrosa curve index; IOP, intraocular pressure; VF, visual field; MD, mean deviation; CT, choroidal thickness; AXL, axial length; PSD, pattern standard deviation.
Figure 2 illustrates partial dependence plots showing the influence of these four variables on SHAP values. There were approximate tipping points where SHAP values started to rapidly decrease for the initial IOP, mean LCCI, and global CT: when the highest IOP reached above 30 mm Hg, the mean LCCI reached above 13, and the global CT reached below 200 µm, the rate of RNFL thinning as determined by SHAP values appeared to be accelerated. On the other hand, the SHAP values rapidly decreased when the VF MD were converged to −5 dB. When the VF MD was under −5 dB, the rate of RNFL thinning became faster with increasing VF MD until it reached −5 dB, whereas it became slower when the VF MD was over −5 dB. 
Figure 2.
 
Partial dependence plots showing SHAP values versus feature values for the four most important variables. (A) Highest IOP during initial 6 months, (B) mean LCCI, (C) VF MD, and (D) global CT. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; VF, visual field; MD, mean deviation; CT, choroidal thickness.
Figure 2.
 
Partial dependence plots showing SHAP values versus feature values for the four most important variables. (A) Highest IOP during initial 6 months, (B) mean LCCI, (C) VF MD, and (D) global CT. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; VF, visual field; MD, mean deviation; CT, choroidal thickness.
Figure 3 shows partial interactions of initial IOP, global CT, and mean LCCI. VF MD was not included in the analysis of interactions, because its nonlinear influence on the rate of RNFL thinning resulted in a complicated interactions with other variables. Larger mean LCCI and larger initial IOP (see Fig. 3A) and larger mean LCCI and smaller global CT (see Fig. 3B) were found to contribute to faster rates of RNFL thinning. Large LCCI was a main contributor for faster RNFL thinning when it was larger than 12 (see Figs. 3A, 3B). When the LCCI was smaller than 12, the rate of RNFL thinning was affected more by higher initial IOP (see Fig. 3A), or smaller global CT (see Fig. 3B). 
Figure 3.
 
Partial interaction plots showing interactions between (A) initial IOP and mean LCCI and (B) global CT and mean LCCI. Overall, larger mean LCCI and larger initial IOP A, and larger mean LCCI and smaller global CT B together contributed to faster rates of RNFL thinning. However, in the region where mean LCCI was >12, the rate of RNFL thinning was mainly dependent on the mean LCCI, whereas higher IOP A and thinner global CT B were main contributors when the mean LCCI was <12. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; CT, choroidal thickness.
Figure 3.
 
Partial interaction plots showing interactions between (A) initial IOP and mean LCCI and (B) global CT and mean LCCI. Overall, larger mean LCCI and larger initial IOP A, and larger mean LCCI and smaller global CT B together contributed to faster rates of RNFL thinning. However, in the region where mean LCCI was >12, the rate of RNFL thinning was mainly dependent on the mean LCCI, whereas higher IOP A and thinner global CT B were main contributors when the mean LCCI was <12. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; CT, choroidal thickness.
Table 2 shows the results of linear regression analyses of factors significantly associated with the rate of RNFL thinning. Univariate analysis showed that larger initial IOP, glaucoma surgery during the initial 6 months, and larger mean LCCI were significantly associated with faster RNFL thinning (all P ≤ 0.001). Multivariable analysis revealed that larger initial IOP, larger mean LCCI, and smaller global CT were independently associated with rapid RNFL thinning. Scatter plots showing the relationship among the four variables and the rate of RNFL thinning are shown in the Supplementary Figure. Nonlinear correlations assessed using polynomial regression showed significant or marginally significant superiority than linear correlations for all four variables. 
Table 2.
 
Variables Associated With a Faster Rate of Retinal Nerve Fiber Layer Thinning
Table 2.
 
Variables Associated With a Faster Rate of Retinal Nerve Fiber Layer Thinning
The three common significant variables (initial IOP, mean LCCI, and global CT) were included in a decision tree model, which identified four groups of eyes with different rates of RNFL thinning (Fig. 4). The strongest discriminating variable was initial IOP. Eyes with initial IOP >26.5 mm Hg showed the fastest rate of RNFL thinning (−1.76%/year, n = 52), with no further splitting being necessary. In eyes with initial IOP ≤26.5 mm Hg, the next discriminating variable was mean LCCI; eyes with a mean LCCI >13.95 had a faster rate of RNFL thinning (−1.41%/year, n = 53) than eyes with a mean LCCI ≤13.95 (−0.86%/year, n = 607). Global CT was the third splitting variable for eyes with lower initial IOP (≤26.5 mm Hg) and smaller mean LCCI (≤13.95), with eyes having a global CT ≤117.5 µm having a faster rate of RNFL thinning (−1.13%/year, n = 75) than those with a global CT >117.5 µm (−0.85%/year, n = 532). 
Figure 4.
 
Decision tree model stratifying groups with faster or slower glaucoma progression based on factors best explaining the rate of RNFL thinning. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; CT, choroidal thickness.
Figure 4.
 
Decision tree model stratifying groups with faster or slower glaucoma progression based on factors best explaining the rate of RNFL thinning. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; CT, choroidal thickness.
Discussion
The model described in this study identified ophthalmic factors that could predict the future rate of RNFL thinning for over 5 years in patients with OAG. The model incorporated key features of the ONH that have recently been shown to be representative of the two essential mechanisms of glaucoma. In addition to previously established variables, such as IOP, ONH features were found to be strong predictors of the rate of RNFL thinning. The strongest features in the prediction model were also identified as significant factors in the multivariable regression analysis, validating the accuracy of the prediction model. 
Conventional regression models assume that all relationships are linear, which is not true in real-world applications. The strength of machine learning models is that they consider all potential nonlinear relationships and interactions among features, enabling the development of more realistic prediction models. The RF is recognized as one of the powerful machine learning methods in the development of disease prediction models.3,2629 Strength of the RF over neural network models is that less hyperparameters are selected during training, reducing the chance of overfitting, which makes it easier to establish a model. The major limitation of the RF has been its high complexity, limiting its interpretation. However, the recent development of explainable artificial intelligence (i.e. SHAP) has facilitated the interpretation of the RF. It has been shown that SHAP performs well in explaining machine learning models even when the data is of a moderate size.2,30 Using SHAP, the present study could successfully explain the prediction model, along with the importance and dependence of each feature. The interactions between features could be readily illustrated by partial interaction plots. 
In the present study, a decision tree model was also used for the analysis. A decision tree model itself can illustrate groups having similar characteristic features that may have interactions together affecting the dependent variable. Using the decision tree, we could stratify patients having differed rates of RNFL thinning with different dominant features that may co-contribute to the rate of progressive RNFL thinning. 
Larger mean LCCI and higher IOP were the two most important predictors of rapid RNFL thinning. These two variables were previously shown to be risk factors for glaucoma progression.8,17,3335 Posterior deformation of the LC is the principal event that induces axonal damage in glaucoma.4 LCCI is regarded as a reliable indicator of LC deformation, as it can be used to estimate the degree of mechanical strain on the LC.13,17 The partial interaction between mean LCCI and IOP showed that an LCCI larger than a certain degree (>12) was significantly associated with faster RNFL thinning at any level of IOP, whereas influence of IOP was predominant when the LCCI was smaller than 12. 
VF MD was the third important feature, having nonlinear relationship with the rate of RNFL thinning. When the VF damage was within earlier stage (i.e. VF MD >−5dB), the rate of RNFL thinning became faster with worsening of VF MD, whereas the opposite was true when the VF was moderately or severely damaged (i.e. VF MD < −5dB). This may be attributable to the curvilinear relationship between structure and function in the glaucoma continuum,36 where the rate of structural damage accelerates until functional damage reaches a certain degree, and then decelerates with functional deterioration going further.3739 Because of this nonlinear relationship between structure and function, the influence of VF MD could only be revealed by an RF model, not by a conventional regression model. 
Smaller peripapillary global CT was the fourth predictor of faster RNFL thinning. Both partial interaction plots and the decision tree revealed that global CT was particularly important when the LCCI was smaller. The peripapillary choroid and the ONH have the same vascular origin, with the peripapillary choroid considered an important vascular structure for perfusion of the ONH.5,18 A thinner peripapillary choroid was found to be associated with localized dropout of the choroidal microvessels at the location of glaucomatous damage,20 which has been shown as a significant indicator of future glaucoma progression.8,9,40,41 These findings suggest that decreased CT is associated with decreased ONH perfusion.19,20 The thinner CT plus a smaller LCCI predicting rapid progression may indicate that ischemia plays a more important role than mechanical stress in the pathogenesis of glaucoma in certain eyes. 
The SHAP values for initial IOP, LCCI, and CT in our model rapidly decreased when these parameters were above or below certain levels. These findings indicate that cutoff levels for these parameters are associated with accelerated axonal damage. On the other hand, partial interaction plots showed some exceptions depending on the interactions between variables. Global CT appeared to be more influential than mean LCCI when the latter was below a certain level (approximately 12 in the present study). This finding may indicate the importance of a thin choroid in eyes with smaller LCCI, suggesting that the development and progression of glaucoma in eyes with smaller LCCI could be mainly attributed to a thin choroid (i.e. vascular incompetency). Among glaucomatous eyes with small LCCI, vascular incompetency may be worse in those with smaller LCCI, leading to a more rapid rate of progression. 
The features that can be influenced by a change in IOP, such as LC and CT, were evaluated after the initiation of medical treatment or after surgical treatment when surgery was performed during the initial 6 months. Measurements at these time points avoided the immediate effect of IOP lowering that might affect baseline features. Determining the baseline features in treatment naïve patients could allow the model to predict the natural course of glaucoma from the time before initiation of treatment. However, this may be less feasible in a real-world clinical setting. 
This model had several limitations. First, the prediction model was based on initial features. Glaucoma is a dynamic disease influenced by various factors. Any event during the follow-up, such as IOP fluctuation or change in treatment, could affect the rate of RNFL thinning. Therefore, it may not be practically possible to develop a highly accurate prediction model based only on initial features. Second, ONH parameters were measured manually by methods that are time-consuming and subject to variability. Development of a reliable algorithm that enables automated segmentation and analysis of imaging data could enhance the usefulness of the prediction model in patients diagnosed with glaucoma. There are numerous variables that can possibly change the rate of RNFL thinning during the follow-up. 
In conclusion, a machine learning approach to developing a prediction model enabled the identification and interpretation of important ophthalmic features suggestive of future glaucoma progression. These features also allowed patients with glaucoma to be classified by risk of disease progression. Predictive modeling may facilitate automated risk prediction, which may help in determining mechanisms underlying the pathogenesis of glaucoma and establishing long-term treatment strategies more efficiently in individual patients with glaucoma. 
Acknowledgments
Supported by the Patient-Centered Clinical Research Coordinating Center, funded by the Ministry of Health & Welfare, Republic of Korea (grant no. HI19C0481, HC19C0276), and by the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science, and Technology (No. 2016R1D1A1B02011696). 
Disclosure: E.J. Lee, None; T.-W. Kim, None; J.-A. Kim, None; S.H. Lee, None; H. Kim, None 
References
Mehta P, Petersen CA, Wen JC, et al. Automated Detection of Glaucoma With Interpretable Machine Learning Using Clinical Data and Multimodal Retinal Images. Am J Ophthalmol. 2021; 231: 154–169. [CrossRef] [PubMed]
Oh S, Park Y, Cho KJ, Kim SJ. Explainable Machine Learning Model for Glaucoma Diagnosis and Its Interpretation. Diagnostics (Basel). 2021; 11(3): 510. [CrossRef] [PubMed]
Baxter SL, Marks C, Kuo TT, Ohno-Machado L, Weinreb RN. Machine Learning-Based Predictive Modeling of Surgical Intervention in Glaucoma Using Systemic Data From Electronic Health Records. Am J Ophthalmol. 2019; 208: 30–40. [CrossRef] [PubMed]
Burgoyne CF, Downs JC, Bellezza AJ, Suh JK, Hart RT. The optic nerve head as a biomechanical structure: a new paradigm for understanding the role of IOP-related stress and strain in the pathophysiology of glaucomatous optic nerve head damage. Prog Retin Eye Res. 2005; 24: 39–73. [CrossRef] [PubMed]
Hayreh SS. The blood supply of the optic nerve head and the evaluation of it - myth and reality. Prog Retin Eye Res. 2001; 20: 563–593. [CrossRef] [PubMed]
Lee EJ, Kim TW. Lamina cribrosa reversal after trabeculectomy and the rate of progressive retinal nerve fiber layer thinning. Ophthalmology. 2015; 122: 2234–2242. [CrossRef] [PubMed]
Lee EJ, Kim TW, Kim M, Kim H. Influence of lamina cribrosa thickness and depth on the rate of progressive retinal nerve fiber layer thinning. Ophthalmology. 2015; 122: 721–729. [CrossRef] [PubMed]
Lee EJ, Kim TW, Kim JA, et al. Elucidation of the Strongest Factors Influencing Rapid Retinal Nerve Fiber Layer Thinning in Glaucoma. Invest Ophthalmol Vis Sci. 2019; 60: 3343–3351. [CrossRef] [PubMed]
Lee EJ, Kim JA, Kim TW. Influence of choroidal microvasculature dropout on the rate of glaucomatous progression, a prospective study. Ophthalmology Glaucoma. 2020; 3: 25–31. [CrossRef] [PubMed]
Bowd C, Zangwill LM, Weinreb RN, Medeiros FA, Belghith A. Estimating optical coherence tomography structural measurement floors to improve detection of progression in advanced glaucoma. Am J Ophthalmol. 2017; 175: 37–44. [CrossRef] [PubMed]
Mwanza JC, Kim HY, Budenz DL, et al. Residual and dynamic range of retinal nerve fiber layer thickness in glaucoma: comparison of three OCT platforms. Invest Ophthalmol Vis Sci. 2015; 56: 6344–6351. [CrossRef] [PubMed]
Lee EJ, Kim TW, Kim H, Lee SH, Girard MJ, Mari JM. Comparison between lamina cribrosa depth and curvature as a predictor of progressive retinal nerve fiber layer thinning in primary open-angle glaucoma. Ophthalmology Glaucoma. 2018; 1: 44–51. [CrossRef] [PubMed]
Lee SH, Kim TW, Lee EJ, Girard MJ, Mari JM. Diagnostic power of lamina cribrosa depth and curvature in glaucoma. Invest Ophthalmol Vis Sci. 2017; 58: 755–762. [CrossRef] [PubMed]
Lee EJ, Kim TW, Weinreb RN. Reversal of lamina cribrosa displacement and thickness after trabeculectomy in glaucoma. Ophthalmology. 2012; 119: 1359–1366. [CrossRef] [PubMed]
Lee EJ, Kim TW, Weinreb RN, Kim H. Reversal of lamina cribrosa displacement after intraocular pressure reduction in open-angle glaucoma. Ophthalmology. 2013; 120: 553–559. [CrossRef] [PubMed]
Lee SH, Yu DA, Kim TW, Lee EJ, Girard MJ, Mari JM. Reduction of the Lamina Cribrosa Curvature After Trabeculectomy in Glaucoma. Invest Ophthalmol Vis Sci. 2016; 57: 5006–5014. [CrossRef] [PubMed]
Lee EJ, Kim TW, Kim H, Lee SH, Girard MJA, Mari JM. Comparison between Lamina Cribrosa Depth and Curvature as a Predictor of Progressive Retinal Nerve Fiber Layer Thinning in Primary Open-Angle Glaucoma. Ophthalmol Glaucoma. 2018; 1: 44–51. [CrossRef] [PubMed]
Lee KM, Kim JM, Lee EJ, Kim TW. Anterior Optic Nerve Head Perfusion is Dependent on Adjacent Parapapillary Choroidal perfusion. Sci Rep. 2019; 9: 10999. [CrossRef] [PubMed]
Lee SH, Lee EJ, Kim TW. Topographic correlation between juxtapapillary choroidal thickness and microstructure of parapapillary atrophy. Ophthalmology. 2016; 123: 1965–1973. [CrossRef] [PubMed]
Lee SH, Lee EJ, Kim TW. Topographic correlation between juxtapapillary choroidal thickness and parapapillary deep-layer microvasculature dropout in primary open-angle glaucoma. Br J Ophthalmol. 2018; 102: 1134–1140. [CrossRef] [PubMed]
Sullivan-Mee M, Patel NB, Pensyl D, Qualls C. Relationship Between Juxtapapillary Choroidal Volume and Beta-Zone Parapapillary Atrophy in Eyes With and Without Primary Open-Angle Glaucoma. Am J Ophthalmol. 2015; 160: 637–647.e631. [CrossRef] [PubMed]
Kara N, Baz O, Altan C, Satana B, Kurt T, Demirok A. Changes in choroidal thickness, axial length, and ocular perfusion pressure accompanying successful glaucoma filtration surgery. Eye (Lond). 2013; 27: 940–945. [CrossRef] [PubMed]
Kadziauskiene A, Kuoliene K, Asoklis R, Lesinskas E, Schmetterer L. Changes in choroidal thickness after intraocular pressure reduction following trabeculectomy. Acta Ophthalmol. 2016; 94: 586–591. [CrossRef] [PubMed]
Breiman L. Random forests. Machine learning. 2001; 45: 5–32. [CrossRef]
Lundberg SM, Lee S-I. A unified approach to interpreting model predictions. Proceedings of the 31st International Conference on Neural Information Processing Systems. 2017: 4768–4777.
Sugimoto K, Murata H, Hirasawa H, Aihara M, Mayama C, Asaoka R. Cross-sectional study: Does combining optical coherence tomography measurements using the “Random Forest” decision tree classifier improve the prediction of the presence of perimetric deterioration in glaucoma suspects? BMJ Open. 2013; 3: e003114. [CrossRef] [PubMed]
Nezu N, Usui Y, Saito A, et al. Machine Learning Approach for Intraocular Disease Prediction Based on Aqueous Humor Immune Mediator Profiles. Ophthalmology. 2021; 128: 1197–1208. [CrossRef] [PubMed]
Horn FK, Lammer R, Mardin CY, et al. Combined evaluation of frequency doubling technology perimetry and scanning laser ophthalmoscopy for glaucoma detection using automated classification. J Glaucoma. 2012; 21: 27–34. [CrossRef] [PubMed]
Muhammad H, Fuchs TJ, De Cuir N, et al. Hybrid Deep Learning on Single Wide-field Optical Coherence tomography Scans Accurately Classifies Glaucoma Suspects. J Glaucoma. 2017; 26: 1086–1094. [CrossRef] [PubMed]
Qiao N, Ma Y, Chen X, et al. Machine Learning Prediction of Visual Outcome after Surgical Decompression of Sellar Region Tumors. J Pers Med. 2022; 12(2): 152. [CrossRef] [PubMed]
Kim YW, Lee EJ, Kim TW, Kim M, Kim H. Microstructure of beta-zone parapapillary atrophy and rate of retinal nerve fiber layer thinning in primary open-angle glaucoma. Ophthalmology. 2014; 121: 1341–1349. [CrossRef] [PubMed]
Leung CK, Cheung CY, Weinreb RN, et al. Evaluation of retinal nerve fiber layer progression in glaucoma: a study on optical coherence tomography guided progression analysis. Invest Ophthalmol Vis Sci. 2010; 51: 217–222. [CrossRef] [PubMed]
The effectiveness of intraocular pressure reduction in the treatment of normal-tension glaucoma. Collaborative Normal-Tension Glaucoma Study Group. Am J Ophthalmol. 1998; 126: 498–505. [CrossRef] [PubMed]
[No authors listed] . The Advanced Glaucoma Intervention Study (AGIS): 7. The relationship between control of intraocular pressure and visual field deterioration. The AGIS Investigators. Am J Ophthalmol. 2000; 130: 429–440. [PubMed]
Leske MC, Heijl A, Hussein M, et al. Factors for glaucoma progression and the effect of treatment: the early manifest glaucoma trial. Arch Ophthalmol. 2003; 121: 48–56. [CrossRef] [PubMed]
Weinreb RN, Friedman DS, Fechtner RD, et al. Risk assessment in the management of patients with ocular hypertension. Am J Ophthalmol. 2004; 138: 458–467. [CrossRef] [PubMed]
Leung CK, Medeiros FA, Zangwill LM, et al. American Chinese glaucoma imaging study: a comparison of the optic disc and retinal nerve fiber layer in detecting glaucomatous damage. Invest Ophthalmol Vis Sci. 2007; 48: 2644–2652. [CrossRef] [PubMed]
Schlottmann PG, De Cilla S, Greenfield DS, Caprioli J, DF Garway-Heath. Relationship between visual field sensitivity and retinal nerve fiber layer thickness as measured by scanning laser polarimetry. Invest Ophthalmol Vis Sci. 2004; 45: 1823–1829. [CrossRef] [PubMed]
Leung CK, Chong KK, Chan WM, et al. Comparative study of retinal nerve fiber layer measurement by StratusOCT and GDx VCC, II: structure/function regression analysis in glaucoma. Invest Ophthalmol Vis Sci. 2005; 46: 3702–3711. [CrossRef] [PubMed]
Jo YH, Shin JW, Song MK, Won HJ, Kook MS. Baseline Choroidal Microvasculature Dropout as a Predictor of Subsequent Visual Field Progression in Open-angle Glaucoma. J Glaucoma. 2021; 30: 672–681. [CrossRef] [PubMed]
Park HL, Kim JW, Park CK. Choroidal microvasculature dropout is associated with progressive retinal nerve fiber layer thinning in glaucoma with disc hemorrhage. Ophthalmology. 2018; 125: 1003–1013. [CrossRef] [PubMed]
Figure 1.
 
Interpretation of the final model based on baseline patient variables. (A) Feature importance plot based on mean SHAP values. (B) Interpretation of the importance of features using the SHAP plot. The red and blue colors indicate feature values of high and low levels, respectively. For example, a larger mean LCCI had a strong and negative impact on the rate of RNFL thinning (i.e. faster RNFL thinning). LCCI, lamina cribrosa curve index; IOP, intraocular pressure; VF, visual field; MD, mean deviation; CT, choroidal thickness; AXL, axial length; PSD, pattern standard deviation.
Figure 1.
 
Interpretation of the final model based on baseline patient variables. (A) Feature importance plot based on mean SHAP values. (B) Interpretation of the importance of features using the SHAP plot. The red and blue colors indicate feature values of high and low levels, respectively. For example, a larger mean LCCI had a strong and negative impact on the rate of RNFL thinning (i.e. faster RNFL thinning). LCCI, lamina cribrosa curve index; IOP, intraocular pressure; VF, visual field; MD, mean deviation; CT, choroidal thickness; AXL, axial length; PSD, pattern standard deviation.
Figure 2.
 
Partial dependence plots showing SHAP values versus feature values for the four most important variables. (A) Highest IOP during initial 6 months, (B) mean LCCI, (C) VF MD, and (D) global CT. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; VF, visual field; MD, mean deviation; CT, choroidal thickness.
Figure 2.
 
Partial dependence plots showing SHAP values versus feature values for the four most important variables. (A) Highest IOP during initial 6 months, (B) mean LCCI, (C) VF MD, and (D) global CT. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; VF, visual field; MD, mean deviation; CT, choroidal thickness.
Figure 3.
 
Partial interaction plots showing interactions between (A) initial IOP and mean LCCI and (B) global CT and mean LCCI. Overall, larger mean LCCI and larger initial IOP A, and larger mean LCCI and smaller global CT B together contributed to faster rates of RNFL thinning. However, in the region where mean LCCI was >12, the rate of RNFL thinning was mainly dependent on the mean LCCI, whereas higher IOP A and thinner global CT B were main contributors when the mean LCCI was <12. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; CT, choroidal thickness.
Figure 3.
 
Partial interaction plots showing interactions between (A) initial IOP and mean LCCI and (B) global CT and mean LCCI. Overall, larger mean LCCI and larger initial IOP A, and larger mean LCCI and smaller global CT B together contributed to faster rates of RNFL thinning. However, in the region where mean LCCI was >12, the rate of RNFL thinning was mainly dependent on the mean LCCI, whereas higher IOP A and thinner global CT B were main contributors when the mean LCCI was <12. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; CT, choroidal thickness.
Figure 4.
 
Decision tree model stratifying groups with faster or slower glaucoma progression based on factors best explaining the rate of RNFL thinning. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; CT, choroidal thickness.
Figure 4.
 
Decision tree model stratifying groups with faster or slower glaucoma progression based on factors best explaining the rate of RNFL thinning. IOP, intraocular pressure; LCCI, lamina cribrosa curve index; CT, choroidal thickness.
Table 1.
 
Baseline Characteristics or Participants (n = 712)
Table 1.
 
Baseline Characteristics or Participants (n = 712)
Table 2.
 
Variables Associated With a Faster Rate of Retinal Nerve Fiber Layer Thinning
Table 2.
 
Variables Associated With a Faster Rate of Retinal Nerve Fiber Layer Thinning
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×