## Abstract

As the prevalence of diabetic retinopathy (DR) continues to rise, there is a need to develop computer-aided screening methods. The current study reports and validates an ordinary least squares (OLS) method to model optical coherence tomography angiography (OCTA) images and derive OLS parameters for classifying proliferative DR (PDR) and no/mild non-proliferative DR (NPDR) from non-diabetic subjects. OLS parameters were correlated with vessel metrics quantified from OCTA images and were used to determine predicted probabilities of PDR, no/mild NPDR, and non-diabetics. The classification rates of PDR and no/mild NPDR from non-diabetic subjects were 94% and 91%, respectively. The method had excellent predictive ability and was validated. With further development, the method may have potential clinical utility and contribute to image-based computer-aided screening and classification of stages of DR and other ocular and systemic diseases.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Diabetic retinopathy (DR) remains one of the leading causes of vision loss worldwide [1]. It is the most common microvascular complication of diabetes and is classified into non-proliferative DR (NPDR) and proliferative DR (PDR) [2,3]. It is projected the global number of people with DR will increase to 191 million by 2030 [4]. Identified by neovascularization and vitreous or preretinal hemorrhage, PDR and its complications form the most severe types of the disease, potentially leading to blindness [3]. Progression of DR may be asymptomatic, thus early screening of DR is crucial in reducing the risk of visual impairment [5].

Detection of retinal microvasculopathies through evaluation of fundus photographs and ophthalmoscopy are the most common methods to diagnose DR [6]. Furthermore, analysis of retinal microvasculature in digital fundus images facilitates prognosis of DR, as quantitative assessment of vessel morphology is feasible [7]. Increased retinal vessel tortuosity, decreased vessel density, and enlargement of the foveal avascular zone have been detected in DR optical coherence tomography angiography (OCTA) images [8–10].

Innovative methods involving machine learning and computational techniques have been developing to improve screening of DR. Deep learning algorithms have been designed to detect DR with high validity using retinal photographs [11–15]. Discrimination of DR images has been made possible by texture classification [16,17] and lesion detection in retinal images [18–20]. Furthermore, support vector machine classification methods have effectively detected DR using retinal images in recent studies [21–23]. More recently, features of OCTA images have been utilized in establishing computer-aided methods to classify mild NPDR [24,25].

The use of automated methods to screen DR will increase efficiency, reduce costs, and potentially improve patient outcomes by early detection of the disease [26]. As DR prevalence continues to rise [4], techniques to advance diagnosis must continue to be developed. Although vessel abnormalities are expected as DR progresses, development of methods to promptly identify pathology that is not visually recognizable is important to address early screening and diagnosis of DR through simple and objective classification.

We have previously established a fine structure analysis method to discriminate between normal and demented brain images [27], as well as DR stages using conjunctival microvasculature images [28] and retinal images [29]. The current study reports and validates a method of ordinary least squares (OLS) modeling of images for classification of DR from non-diabetic subjects. The method was validated in advanced DR due to established retinal vascular pathologies and applied for classification of early DR. The purposes of the current study were to test the hypotheses that estimated parameters derived from an OLS model applied to OCTA images are related to retinal vessel metrics and can classify advanced and early DR from non-diabetic subjects.

## 2. Materials and methods

#### 2.1 Subjects

The study was conducted at the University of Illinois at Chicago and University of Southern California and was approved by their corresponding Institutional Review Boards. After the study was explained to the subjects, informed consents were obtained in accordance to the tenets of Declaration of Helsinki. Subjects were stratified into non-diabetic (N=22), no DR (N=33), mild NPDR (N=26), and PDR (N=13) groups based on clinical retinal examination by expert retinal specialists. Diabetic subjects with no DR or mild NPDR were grouped as no/mild NPDR (N=59), as categorized in previous publications [30–32]. The ages of non-diabetic, no/mild NPDR, and PDR subjects were 55 ± 10 years (mean ± standard deviation), 56 ± 15 years, and 57 ± 11 years, respectively (P ≥ 0.62). The non-diabetic group consisted of 7 males and 15 females, the no/mild NPDR group included 26 males and 33 females, and there were 8 males and 5 females in the PDR group (P ≥ 0.09). Eleven of the 13 PDR subjects had a previous history of panretinal photocoagulation. None of the mild NPDR subjects had a history of treatment or diabetic macular edema, according to previously published criteria [33].

Prior to imaging, subjects’ eyes were dilated using 1% tropicamide (Alcon Laboratories, Inc., Fort Worth, TX) and 2.5% phenylephrine hydrochloride (Paragon BioTek Inc., Portland, OR). Images acquired from one eye of each subject were analyzed. One eye was selected at random, when data were available in both eyes.

#### 2.2 Image acquisition

Images were obtained using a commercially available Avanti OCTA system (Optovue, Inc., Fremont, CA). OCTA images of the superficial capillary plexus were acquired in a 6 × 6 mm macular region centered on the fovea. The superficial capillary plexus was defined by the Optovue software within the nerve fiber and ganglion cell layers. The superficial capillary plexus was selected for analysis to avoid potential projection artifacts on deeper plexus.

#### 2.3 Image analysis

### 2.3.1 Vessel tortuosity

Retinal vessel tortuosity was quantified by our previously published tortuosity index using 6 × 6 mm OCTA images of the retinal superficial capillary plexus [34]. A binary vessel map was produced from detection of retinal vessels using a k-means clustering algorithm. Distance transformation was used to extract centerlines between the bifurcation points by selection of vessel endpoints on the binary vessel map. The mathematical derivation of tortuosity index can be found elsewhere [34]. Tortuosity index was calculated for each of the extracted centerlines and averaged per eye. The minimum value for tortuosity index is theoretically zero, which corresponds to a straight line, and a theoretical maximum value for tortuosity index does not exist.

### 2.3.2 Vessel density and spacing

A previously validated local fractal dimension analysis method was applied to OCTA images to assess retinal vessel density [35–37]. To calculate the local fractal dimension of each pixel, a moving window size of 3 × 3 pixels was utilized, which varied with the distribution of vessels around it. The fractal dimension ratio (FDR) was derived from the ratio of local fractal dimension of each pixel to the maximum local fractal dimension. A graphic representation of the probability index of presence of vessel of a certain size at each pixel is provided by the FDR. FDR between 0.7 and 1 represented large and small vessels (vessel density), between 0.3 and 0.7 represented small vessel spacing, and less than 0.3 represented large vessel spacing, as previously reported [35]. Metric values of 0 and 1 indicate 0% and 100% of the total image, respectively.

#### 2.4 Image model assumptions

Our fundamental assumption is that a two-dimensional (2D) digital image with linearly independent columns and rows can be a solution to a 2D partial difference equation (PdE). Based on differential equation parameter estimation literature [38], a 2D PdE with estimable parameters can be derived by finite difference approximation, employing backward differences of all derivatives in the general linear stationary 2D partial differential equation. The PdE has the same autocorrelation function and parameters as the image. As described in our previous publications [27–29], Eq. (1) expresses the intensity (y) of the image at pixel locations (i,j) as a solution of the general linear, stationary, autoregressive PdE.

In Eq. (1), k + l>0 and u(i,j) is a random process error term which can be minimized in variance to estimate the β_{kl} parameters. In fact, it has been shown that a majority of images have an autoregressive 2D autocorrelation function [39]. We propose that the identified parameters in the autoregressive model are informationally sufficient to discriminate normal and diseased conditions. This is a classic binary discrimination exercise [40].

#### 2.5 Ordinary least squares model derivation and estimation

Derivation of the OLS model parameters has been described in depth by O’Neill and associates [27]. An image of m by n pixels, expressed as a matrix, can be transformed into a one dimension vector of length m x n by 1 by stacking the pixels into a single column. Specifically, if A is any m by n matrix [a_{1} a_{2 . . . }a_{n} ], then vec(A)= [a_{1}^{T} a_{2}^{T}_{ . . . }a_{n}^{T} ]^{T} is a mn by 1 vector. The following result transforms estimation of the autoregressive model parameters, β_{kl}, into a simple OLS estimation exercise [41].

If c_{t} and A_{t} are real scalars and matrices, respectively, then

_{o}= vec[y(i,j)], x

_{1}= vec[y(i,j-1)], x

_{2}= vec[y(i-1,j)], …, x

_{p+q+qp}= vec[y(i-p,j-q)], and u

_{o}= vec[u(i,j)]. Following Eq. (2), the expression in Eq. (1), with β as a column vector having p + q+qp elements, β

_{kl}, has the familiar OLS population regression form:

_{o}estimates the random residual vector u

_{o}. A necessary and sufficient condition that a vector b exists to minimize e

_{o}

^{T}e

_{o}is the columns of X are linearly independent; in which case b = (X

^{T}X)

^{-1}X

^{T}y

_{o}and e

_{o}= y

_{o}- Xb.

In our previous method for discriminating groups based on Fisher’s linear discriminant (FLD) method [27–29], 9 OLS parameters were used. Discrimination efficacy increases with the number of OLS parameters per image. However, there are validation constraints. In the FLD method, B_{1} and B_{2} denote the class matrices of parameter vectors of m_{1} and m_{2} from respective images. Each vector is of parameter size n (e.g., n = 3 for a spatial lag of 1), such that the class matrices are m_{1} by n and m_{2} by n. The estimated covariance matrices of the B matrices are expressed by G_{1} and G_{2} and the estimated covariance matrix of B_{p} = [B_{1} B_{2}]^{T} is G_{p}. The eigenvector v_{c} satisfying: (G_{1} – G_{2} – G_{p})v_{c}= λ(G_{1} + G_{2})v_{c}, for the unique eigenvalue λ ≠ 0, is the optimal projection vector for discriminating images from binary groups. A solution for v_{c} exists only if the number of OLS parameters per image is less than or equal to the total number of images minus 2 [42]. Also, the estimated covariance matrices must be positive definite for the FLD projections maximum separation to be valid [43]. Both of these constraints become tight as the number of OLS parameters increases. To validate the models, the order s of the assumed PDE model requires (s+1)^{2} OLS parameters to be determined within the above noted constraints. The algorithm to select s is quite simple. Increase s from one until any covariance matrix is not positive definite or the number of OLS parameters is greater than the number of images minus 2. For the image data in the current study, s=4, thus, (s+1)^{2} = 25 OLS parameters were required to estimate a 4^{th} order PDE model for each image. Note that this algorithm requires modeling all of the images for each s selected. Thus, to reach s=4 requires 4 passes through the data. The OLS and matrix positive definite were checked in MATLAB which typically takes 4.3 seconds to run through 100 images.

In summary, the difference equation model of an image becomes a linear relation between the image, now a vector, and vectors representing spatial lags of the image, reducing the model to an OLS format. The vectorizing transformation of the image allows for estimation of the OLS model parameters, which weight a linear combination of lagged image vectors that minimizes the mean square error between the original image and the weighted lagged images. The resulting OLS parameters are estimates of the autoregressive model parameters, which under necessary conditions exist to minimize random residuals. The vectors of the OLS model are unstacked to produce a 2D OLS image model. Figure 1 illustrates OCTA and OLS images of non-diabetic and PDR subjects. R^{2} derived from the OLS regression for the NC and PDR OLS image models were 71% and 79%, respectively.

#### 2.6 Class probability estimates

The method of logistic regression can be used to estimate the probability a given image is in the diseased or normal class. Logistic regression is a nonlinear transformation of a linear discriminant for class membership. The estimated discriminating variable is the probability a given subject is in a certain class. The probability a particular subject is in class one (${\pi _i}\; $ = 1) is in agreement with the logistic model [44]. For a set of data to be classified into 2 classes, $\pi $ is either 0 or 1. If the predictor variables are known for all data, then the intercept coefficients can be estimated by iteratively reweighted least squares [45].

A natural choice of predictors when estimating the probability for classifying images are the 25 OLS parameters that are derived from each image. Even though iteratively reweighted least squares is a nonlinear algorithm, it has an overfitting problem (analogous to OLS estimation) which requires a careful choice of which OLS parameters to include in the logistic regression. The OLS parameters passed a Kolmogorov-Smirnov (KS) test for normality at level 0.05 or better for all images modeled [43]. The OLS regressions producing the parameters had 156,816 degrees of freedom (6272 samples per parameter estimated), thus the KS test outcomes were not unexpected. Student-t tests were made of all parameters after applying the White heteroscedastic transformation to the estimated parameter covariance matrices [46].

#### 2.7 Statistical analysis

All statistical analyses were conducted using SAS software (SAS, version 9.4; SAS Institute Inc., Cary, NC). A two-sided p-value less than 0.05 was considered to be statistically significant. From the 25 OLS parameters, the first coefficient was excluded, as it represents the intercept. Normality of data distribution of continuous variables and 24 OLS parameters was assessed using Shapiro-Wilk tests and graphical visualization of quantile-quantile plots. Assessment of Cook’s distance, difference in fits, and studentized residuals was performed to detect outliers, and no influential outliers were identified. Comparison of vessel metrics between non-diabetic and PDR groups was performed using unpaired t-tests or Wilcoxon Rank-Sum tests. Pearson or Spearman’s rank correlations were used to determine associations of OLS parameters with each vessel metric based on compiled data in non-diabetic and PDR groups.

Univariate logistic regressions were performed to examine the associations of individual OLS parameters with disease groups. OLS parameters that were significantly associated with presence of PDR or no/mild NPDR (significant Wald and likelihood ratio test) were then selected as predictors in binary multivariate logistic regression with disease group as outcome.

Performance of the binary multivariate logistic regression was assessed by likelihood ratio test, deviance, and Nagelkerke coefficient of determination. Classification rate at a cut-off point of 0.5 and area under receiver operating characteristic curve (AUROC) were determined to evaluate the predictive ability.

Since the sample size in PDR group was not adequate for external validation, the binary multivariate logistic regression model was validated by means of leave-one-out cross validation. Leave-one-out cross validation randomly divides the data set into N partitions, uses N – 1 partitions as the training set, and makes a prediction on the remaining partitions until there is a prediction for N partitions [47]. Cross-validated predicted probabilities were calculated for N response levels, corresponding to the maximum number of response levels across both groups. The cross-validated predicted probabilities were used to generate a receiver operating characteristic curve and to produce the sensitivity and specificity values used in estimating the AUROC. The equality of the AUROCs of the fitted models with and without cross validation was tested by a global chi-square test.

The OLS model was externally validated in the combined non-diabetic and no/mild NPDR groups, due to a larger sample size. Images in 22 non-diabetic and 59 no/mild NPDR subjects (N=81) were randomly split into training (70%; N=58) and validation (30%; N=23) sets. There were 16 non-diabetic and 42 no/mild NPDR subjects in the training set, and 6 non-diabetic and 17 no/mild NPDR subjects in the validation set. Classification rates and AUROCs derived from training and validation sets were compared to assess predictive ability. To determine the model’s performance using the external validation method, calibration, which is a measure of the agreement between observed and predicted outcomes, was assessed in training and validation sets using Hosmer-Lemeshow test and calibration plots [48].

## 3. Results

#### 3.1 Comparison of vessel metrics between non-diabetic and PDR groups

Figure 2 displays examples of OCTA and FDR images of non-diabetic and PDR subjects. The mean and standard deviation of each vessel metric in non-diabetic and PDR groups are presented in Table 1. Vessel density and small vessel spacing were significantly decreased (P < 0.0001) and large vessel spacing was increased (P < 0.0001) in the PDR compared to non-diabetic group. There were no significant differences in tortuosity index between groups (P = 0.34).

#### 3.2 Correlation of OLS parameters with vessel metrics

Association of OLS parameters with vessel metrics was examined in compiled data in PDR and non-diabetic subjects. The strength of the relationships of OLS parameters that were correlated with at least 1 vessel metric are shown in Table 2. Vessel density was correlated with 4 OLS parameters (b6, b13, b17, b20) (P ≤ 0.04, N = 35). Small vessel spacing was correlated with 3 OLS parameters (b2, b6, b7) (P ≤ 0.0002, N = 35). Large vessel spacing was correlated with 5 OLS parameters (b6, b7, b17, b20, b24) (P ≤ 0.04, N = 35). Tortuosity index was correlated with 4 OLS parameters (b4, b5, b15, b23) (P ≤ 0.04, N = 35).

#### 3.3 Classification of non-diabetic and PDR groups

The statistical results of the associations of individual OLS parameters with disease group are reported in Table 3. Six OLS parameters (b2, b6, b11, b20, b23, b24) were significantly associated with presence of PDR (P ≤ 0.005). The individual OLS parameters correctly classified 71% to 77% of subjects and achieved AUROCs between 0.74 and 0.85.

Table 3 additionally presents the statistical results of associations of 6 significant OLS parameters with disease group. The binary multivariate logistic regression model had high performance as indicated by the likelihood ratio test statistic and correctly classified 94% of subjects, such that only 1 of 22 non-diabetic and 1 of 13 PDR subjects were misclassified. Furthermore, the AUROC was 0.99, implying excellent predictive ability.

#### 3.4 Validation of the binary multivariate logistic regression model

The AUROCs of the models with and without cross validation were 0.86 and 0.99, respectively. The difference between AUROCs was not statistically significant (X^{2} = 3.34, P = 0.07), implying validation of the method. However, further studies with a larger sample size are needed to dismiss the possibility of a type II statistical error in this result. The models with and without cross validation correctly classified 80% and 94% of subjects, respectively.

#### 3.5 Classification of non-diabetic and no/mild NPDR groups

In the training set, 5 OLS parameters (b14, b19, b20, b23, b24) were significantly associated with presence of no/mild NPDR (P ≤ 0.012). The individual OLS parameters correctly classified 71% to 79% of subjects and achieved AUROCs between 0.71 and 0.77. The binary multivariate logistic regression model correctly classified 86% of subjects, achieved an AUROC of 0.85, and was well-calibrated (X^{2 }= 6.26, P = 0.62). With the validation set and using the same 5 OLS parameters from the training set, the model correctly classified 91% of subjects, achieved an AUROC of 0.93, and was also well-calibrated (X^{2 }= 2.30, P = 0.93).

## 4. Discussion

The pressure of a continually growing population of individuals with vision impairment associated with DR has prompted ingenious techniques to aid in screening of the disease. The emerging field of image-based computer-aided screening has shown to discriminate between DR and healthy retinal images [16,17,22,23,25,28]. Automated screening has a promising future in alleviating the demand for qualified specialists, reducing economic burden, and most importantly, preventing irreversible blindness [26,49,50]. In the current study, an OLS modeling method was used to generate estimated parameters that were associated with vessel metrics and classified PDR and no/mild NPDR from non-diabetic subjects using OCTA images.

Although increased retinal vessel tortuosity has been described as a biomarker of NPDR, it has not been exhibited in subjects with PDR [8,9]. In agreement with Lee and associates [9], tortuosity index did not differ between PDR and non-diabetic groups in the current study. It is plausible that sclerotic alterations due to DR progression [51,52] or changes following panretinal photocoagulation [53] may decrease vessel tortuosity in PDR.

The findings of decreased vessel density and increased large vessel spacing in PDR are consistent with published literature [9,10,35]. Interestingly, contrary to the results of Bhanushali et al., in the present study small vessel spacing was decreased in PDR [35]. However, Bhanushali and associates additionally reported a decrease in small vessel spacing with increasing severity of DR [35]. Lei et al. demonstrated enlarged retinal vessel diameter in capillaries of the PDR compared to non-diabetic group using 3 × 3 mm OCTA images centered on the fovea [54]. Areas of small vessels may be exhibiting a decrease in spacing possibly due to vasodilation.

In the current study, all vessel metrics were correlated with multiple OLS parameters. These relationships suggest image characteristics such as vessel metrics are involved in statistical classification of disease groups by OLS models. Vessel metrics likely influence OLS parameters and their associations with presence of PDR. Correlations of retinal vessel metrics with OLS parameters strengthens the validity of the models.

The OLS parameters used as predictors in binary multivariate logistic regression correctly classified 94% of PDR and non-diabetic subjects and achieved a 0.99 AUROC. Different stages of DR have been discriminated using retinal and conjunctival images by applying a similar OLS method, but using 9 parameters and Fisher linear discriminant analysis [28,29]. However, the current study is the first report of an OLS 25-coefficient modeling method and the use of these parameters to classify and derive predicted probabilities of PDR and non-diabetics. Previous studies describing deep learning approaches to classify stages of DR using retinal images have similarly achieved AUROCs greater than 0.93 [11,12,14,15,55,56]. Despite demonstrating comparable predictive ability, the OLS modeling method is favorable because it does not require thousands of images to train a neural network. Furthermore, machine learning methods used for detection of DR and DR-related lesions predominantly use fundus photography [11,12,14,22,23,57]. Application of computational techniques to OCTA rather than fundus images is presumably more powerful due to improved resolution and detailed quantification of microvasculature [58–60].

Classification of early stage of DR is essential for early detection of disease. In the current study, in addition to mild NPDR, subjects with no DR were included because retinal damage may be present in these subjects without any clinical signs [61]. In the externally validated set, 91% of no/mild NPDR and non-diabetic subjects were correctly classified with an AUROC of 0.93. Discrimination of non-diabetic and early DR subjects has been reported by other methods. Using a support vector machine classifier, Alam and associates achieved an AUROC of 0.92 when classifying OCTA images of non-diabetic and mild NPDR subjects [25]. Additionally, Khansari et al. used a fine structure analysis method applied to retinal fundus images that correctly classified 88% of non-diabetic and no DR subjects [29].

The most quantitative method to date that is directly related to our study is that by Kar and Maity [19]. This study employs a mutual information maximization of complementary optimal filters to effect classification of data sets of 40 to 1200 total subjects with an average of 98% accuracy. Each image required 13 minutes to model. Therefore, their method of image modeling applied to 100 images to determine the optimal PdE order would require 3.61 days compared to our method that takes 20 seconds. Further, image pre-processing and mutual information calculation preclude any possibility of parameter significance tests.

Deep learning methods (DLM), various forms of convolutional neural networks (CNN), dominate other retinal image classification literature [11,12,14,15]. In [11], sets of 9963 and 1748 images are subject to a CNN classifier resulting in accuracies of 87.7 to 98.5%. Computer times are not cited, but CNN is a notorious computer resource consumer [62]. CNN also precludes statistical significance statistics except for confidence intervals of estimated sensitivity and specificity. The authors admit unease at the lack of meaning of image features estimated. The image grading system in [12] mimics that in [11] but with a larger, 71043 images -12,239 DR positive, database. A 98.5-92.5% sensitivity – specificity is achieved as is a “black box” dilemma attempting to interpret image properties estimated. The DLM employed in [14] achieved virtually identical results as those in [12] with a slightly larger data base of 75,137 images. It is not clear how the outcomes could possibly be differentiated. A DLM is also the classifier of choice in [15] but a salient contribution is achieved by use of multiethnic subject classes. For each class the shortcomings of using a DLM are obvious and separating classes requires human intervention so it is not clear how “automatic” the process is.

For specific classification outcomes in specificity, sensitivity, and area under the receiver operating characteristic, our method is comparable to those reviewed for a given sample size. For smaller number of subjects, DLMs are not feasible while OLS models exist for any size of image subjects having linearly independent pixel columns. This condition also implies an image has a statistically significant 2D autocorrelation function as opposed to a degenerate one. OLS models may be preferable since parameters are: 1) estimated over millisecond time frames using a modest number of images, 2) normally distributed and Student-t testable, and 3) interpretable in terms of biological features presented on images.

In the current study, only 45% of the OLS parameters were associated with presence of PDR or correlated with vessel metrics. This finding suggests that other unknown image factors besides vessel metrics contribute more than 50% to determination of OLS parameters. The study was limited to classifying based on vascular perfusion in the parafoveal region, though the method may be applicable to images obtained from other modalities and retinal regions. Additionally, the method was established for classifying binary groups and further studies are needed to extend the application of the method for classifying multiple groups. Further development is needed to gain knowledge on the basis of disease classification using OLS parameters. Although a limited sample size prohibited external validation of the model for classifying PDR and non-diabetic subjects, leave-one-out cross validation which is an established and reliable method of validation was employed [47,63]. Nevertheless, despite the small sample size, the number of available images exceeded the minimum number of images needed for OLS modeling. Future studies with larger cohorts are needed to perform external validation of the models.

In conclusion, an OLS modeling method estimated parameters that correlated with retinal vessel metrics and classified PDR and no/mild NPDR from non-diabetic subjects with excellent predictive ability. Compared to other methods of classification, the OLS modeling method has advantages of requiring smaller data sets and shorter processing time, as well as generating results that are amenable to statistical testing and biological interpretation. With further development, the method may have potential clinical utility and contribute to image-based computer-aided screening and classification of stages of DR and other ocular and systemic diseases.

## Funding

National Institute of Diabetes and Digestive and Kidney Diseases (DK104393); National Eye Institute (EY029220, EY030115).

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **W. H. Organization, “Blindness and vision impairment,” https://www.who.int/news-room/fact-sheets/detail/blindness-and-visual-impairment.

**2. **N. Cheung, P. Mitchell, and T. Y. Wong, “Diabetic retinopathy,” Lancet **376**(9735), 124–136 (2010). [CrossRef]

**3. **C. P. Wilkinson, F. L. Ferris 3rd, R. E. Klein, P. P. Lee, C. D. Agardh, M. Davis, D. Dills, A. Kampik, R. Pararajasegaram, and J. T. Verdaguer, “Proposed international clinical diabetic retinopathy and diabetic macular edema disease severity scales,” Ophthalmology **110**(9), 1677–1682 (2003). [CrossRef]

**4. **Y. Zheng, M. He, and N. Congdon, “The worldwide epidemic of diabetic retinopathy,” Indian J. Ophthalmol. **60**(5), 428–431 (2012). [CrossRef]

**5. **R. A. Gangwani, J. X. Lian, S. M. McGhee, D. Wong, and K. K. Li, “Diabetic retinopathy screening: global and local perspective,” Hong Kong Med. J. **22**, 486–495 (2016). [CrossRef]

**6. **M. M. Nentwich and M. W. Ulbig, “Diabetic retinopathy - ocular complications of diabetes mellitus,” World J. Dairy Food Sci. **6**(3), 489–499 (2015). [CrossRef]

**7. **A. W. Stitt, T. M. Curtis, M. Chen, R. J. Medina, G. J. McKay, A. Jenkins, T. A. Gardiner, T. J. Lyons, H. P. Hammes, R. Simo, and N. Lois, “The progress in understanding and treatment of diabetic retinopathy,” Prog. Retinal Eye Res. **51**, 156–186 (2016). [CrossRef]

**8. **M. B. Sasongko, T. Y. Wong, T. T. Nguyen, C. Y. Cheung, J. E. Shaw, and J. J. Wang, “Retinal vascular tortuosity in persons with diabetes and diabetic retinopathy,” Diabetologia **54**(9), 2409–2416 (2011). [CrossRef]

**9. **H. Lee, M. Lee, H. Chung, and H. C. Kim, “Quantification of Retinal Vessel Tortuosity In Diabetic Retinopathy Using Optical Coherence Tomography Angiography,” Retina **38**(5), 976–985 (2018). [CrossRef]

**10. **C. C. Hsiao, H. M. Hsu, C. M. Yang, and C. H. Yang, “Correlation of retinal vascular perfusion density with dark adaptation in diabetic retinopathy,” Graefe's Arch. Clin. Exp. Ophthalmol. **257**(7), 1401–1410 (2019). [CrossRef]

**11. **V. Gulshan, L. Peng, M. Coram, M. C. Stumpe, D. Wu, A. Narayanaswamy, S. Venugopalan, K. Widner, T. Madams, J. Cuadros, R. Kim, R. Raman, P. C. Nelson, J. L. Mega, and D. R. Webster, “Development and Validation of a Deep Learning Algorithm for Detection of Diabetic Retinopathy in Retinal Fundus Photographs,” JAMA **316**(22), 2402–2410 (2016). [CrossRef]

**12. **Z. Li, S. Keel, C. Liu, Y. He, W. Meng, J. Scheetz, P. Y. Lee, J. Shaw, D. Ting, T. Y. Wong, H. Taylor, R. Chang, and M. He, “An Automated Grading System for Detection of Vision-Threatening Referable Diabetic Retinopathy on the Basis of Color Fundus Photographs,” Diabetes Care **41**(12), 2509–2516 (2018). [CrossRef]

**13. **Y. Kanagasingam, D. Xiao, J. Vignarajan, A. Preetham, M. L. Tay-Kearney, and A. Mehrotra, “Evaluation of Artificial Intelligence-Based Grading of Diabetic Retinopathy in Primary Care,” JAMA Netw. Open **1**(5), e182665 (2018). [CrossRef]

**14. **R. Gargeya and T. Leng, “Automated Identification of Diabetic Retinopathy Using Deep Learning,” Ophthalmology **124**(7), 962–969 (2017). [CrossRef]

**15. **D. S. W. Ting, C. Y. Cheung, G. Lim, G. S. W. Tan, N. D. Quang, A. Gan, H. Hamzah, R. Garcia-Franco, I. Y. San Yeo, S. Y. Lee, E. Y. M. Wong, C. Sabanayagam, M. Baskaran, F. Ibrahim, N. C. Tan, E. A. Finkelstein, E. L. Lamoureux, I. Y. Wong, N. M. Bressler, S. Sivaprasad, R. Varma, J. B. Jonas, M. G. He, C. Y. Cheng, G. C. M. Cheung, T. Aung, W. Hsu, M. L. Lee, and T. Y. Wong, “Development and Validation of a Deep Learning System for Diabetic Retinopathy and Related Eye Diseases Using Retinal Images From Multiethnic Populations With Diabetes,” JAMA **318**(22), 2211–2223 (2017). [CrossRef]

**16. **P. Porwal, S. Pachade, M. Kokare, L. Giancardo, and F. Meriaudeau, “Retinal image analysis for disease screening through local tetra patterns,” Comput. Biol. Med. **102**, 200–210 (2018). [CrossRef]

**17. **T. Nazir, A. Irtaza, Z. Shabbir, A. Javed, U. Akram, and M. T. Mahmood, “Diabetic retinopathy detection through novel tetragonal local octa patterns and extreme learning machines,” Artif. Intell. Med. **99**, 101695 (2019). [CrossRef]

**18. **D. Sidibe, I. Sadek, and F. Meriaudeau, “Discrimination of retinal images containing bright lesions using sparse coded features and SVM,” Comput. Biol. Med. **62**, 175–184 (2015). [CrossRef]

**19. **S. S. Kar and S. P. Maity, “Automatic Detection of Retinal Lesions for Screening of Diabetic Retinopathy,” IEEE Trans. Biomed. Eng. **65**(3), 608–618 (2018). [CrossRef]

**20. **P. Chudzik, S. Majumdar, F. Caliva, B. Al-Diri, and A. Hunter, “Microaneurysm detection using fully convolutional neural networks,” Comput. Meth. Prog. Bio. **158**, 185–192 (2018). [CrossRef]

**21. **M. Alam, D. Le, J. I. Lim, R. V. Chan, and X. Yao, “Supervised Machine Learning Based Multi-Task Artificial Intelligence Classification of Retinopathies,” J. Clin. Med. **8**(6), 872 (2019). [CrossRef]

**22. **S. Long, X. Huang, Z. Chen, S. Pardhan, and D. Zheng, “Automatic Detection of Hard Exudates in Color Retinal Images Using Dynamic Threshold and SVM Classification: Algorithm Development and Evaluation,” BioMed Res. Int. **2019**, 3926930 (2019). [CrossRef]

**23. **K. Balasubramanian and N. P. Ananthamoorthy, “Analysis of hybrid statistical textural and intensity features to discriminate retinal abnormalities through classifiers,” Proc. Inst. Mech. Eng., Part H **233**(5), 506–514 (2019). [CrossRef]

**24. **H. S. Sandhu, N. Eladawi, M. Elmogy, R. Keynton, O. Helmy, S. Schaal, and A. El-Baz, “Automated diabetic retinopathy detection using optical coherence tomography angiography: a pilot study,” Br. J. Ophthalmol. **102**(11), 1564–1569 (2018). [CrossRef]

**25. **M. Alam, Y. Zhang, J. I. Lim, R. V. P. Chan, M. Yang, and X. Yao, “Quantitative optical coherence tomography angiography features for objective classification and staging of diabetic retinopathy,” Retina **40**, 1 (2018). [CrossRef]

**26. **C. Y. Cheung, F. Tang, D. S. W. Ting, G. S. W. Tan, and T. Y. Wong, “Artificial Intelligence in Diabetic Eye Disease Screening,” Asia Pac J Ophthalmol (Phila) (2019).

**27. **W. O’Neill, R. Penn, M. Werner, and J. Thomas, “A theory of fine structure image models with an application to detection and classification of dementia,” Quant. Imaging Med. Surg. **5**, 356–367 (2015). [CrossRef]

**28. **M. M. Khansari, W. O’Neill, R. Penn, F. Chau, N. P. Blair, and M. Shahidi, “Automated fine structure image analysis method for discrimination of diabetic retinopathy stage using conjunctival microvasculature images,” Biomed. Opt. Express **7**(7), 2597–2606 (2016). [CrossRef]

**29. **M. M. Khansari, W. D. O’Neill, R. D. Penn, N. P. Blair, and M. Shahidi, “Detection of Subclinical Diabetic Retinopathy by Fine Structure Analysis of Retinal Images,” J. Ophthalmol. **2019**, 1–6 (2019). [CrossRef]

**30. **R. Rajalakshmi, R. Subashini, R. M. Anjana, and V. Mohan, “Automated diabetic retinopathy detection in smartphone-based fundus photography using artificial intelligence,” Eye **32**(6), 1138–1144 (2018). [CrossRef]

**31. **F. Li, Z. Liu, H. Chen, M. Jiang, X. Zhang, and Z. Wu, “Automatic Detection of Diabetic Retinopathy in Retinal Fundus Photographs Based on Deep Learning Algorithm,” Trans. Vis. Sci. Tech. **8**(6), 4 (2019). [CrossRef]

**32. **M. D. Abràmoff, Y. Lou, A. Erginay, W. Clarida, R. Amelon, J. C. Folk, and M. Niemeijer, “Improved Automated Detection of Diabetic Retinopathy on a Publicly Available Dataset Through Integration of Deep Learning,” Invest. Ophthalmol. Visual Sci. **57**(13), 5200–5206 (2016). [CrossRef]

**33. **G. Panozzo, B. Parolini, E. Gusson, A. Mercanti, S. Pinackatt, G. Bertoldo, and S. Pignatto, “Diabetic macular edema: an OCT-based classification,” Semin. Ophthalmol. **19**(1-2), 13–20 (2004). [CrossRef]

**34. **M. M. Khansari, W. O’Neill, J. Lim, and M. Shahidi, “Method for quantitative assessment of retinal vessel tortuosity in optical coherence tomography angiography applied to sickle cell retinopathy,” Biomed. Opt. Express **8**(8), 3796–3806 (2017). [CrossRef]

**35. **D. Bhanushali, N. Anegondi, S. G. Gadde, P. Srinivasan, L. Chidambara, N. K. Yadav, and A. Sinha Roy, “Linking Retinal Microvasculature Features With Severity of Diabetic Retinopathy Using Optical Coherence Tomography Angiography,” Invest. Ophthalmol. Visual Sci. **57**(9), OCT519 (2016). [CrossRef]

**36. **J. Cano, S. Farzad, M. M. Khansari, O. Tan, D. Huang, J. I. Lim, and M. Shahidi, “Relating retinal blood flow and vessel morphology in sickle cell retinopathy,” Eye (Lond) (2019).

**37. **S. G. Gadde, N. Anegondi, D. Bhanushali, L. Chidambara, N. K. Yadav, A. Khurana, and A. Sinha Roy, “Quantification of Vessel Density in Retinal Optical Coherence Tomography Angiography Images Using Local Fractal Dimension,” Invest. Ophthalmol. Visual Sci. **57**(1), 246–252 (2016). [CrossRef]

**38. **A. K. Jain, “Partial differential equations and finite-difference methods in image processing, part 1: Image representation,” J. Optim. Theory Appl. **23**(1), 65–91 (1977). [CrossRef]

**39. **J. J. Koenderink, “The structure of images,” Biol. Cybern. **50**(5), 363–370 (1984). [CrossRef]

**40. **D. H. Moore, “Evaluation of Five Discrimination Procedures for Binary Variables,” J. Am. Stat. Assoc. **68**(342), 399–404 (1973). [CrossRef]

**41. **A. Ben-Israel and T. N. E. Greville, * Generalized Inverses: Theory and Applications* (Springer-Verlag, 2003).

**42. **W. J. Krzanowski, P. Jonathan, W. V. McCarthy, and M. R. Thomas, “Discriminant Analysis with Singular Covariance Matrices: Methods and Applications to Spectroscopic Data,” Bull. Pure Appl. Sci. **44**(1), 101–115 (1995). [CrossRef]

**43. **S. S. Wilks, * Mathematical Statistics* (John Wiley & Sons, Ltd, 1963).

**44. **D. Kleinbaum and M. Klein, * Logistic Regression: A Self-Learning Text* (Springer, New Yorl, 2010). [CrossRef]

**45. **K. P. Murphy, * Machine Learning: A Probabilistic Perspective* (MIT Press, 2012).

**46. **H. White, “A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity,” Econometrica **48**(4), 817–838 (1980). [CrossRef]

**47. **C. Rushing, A. Bulusu, H. Hurwitz, A. B. Nixon, and H. Pang, “A leave-one-out cross validation SAS macro for the identification of markers associated with survival,” Comput. Biol. Med. **57**, 123–129 (2015). [CrossRef]

**48. **E. W. Steyerberg, A. J. Vickers, N. R. Cook, T. Gerds, M. Gonen, N. Obuchowski, M. J. Pencina, and M. W. Kattan, “Assessing the Performance of Prediction Models: A Framework for Traditional and Novel Measures,” Epidemiology **21**(1), 128 (2010). [CrossRef]

**49. **M. F. Norgaard and J. Grauslund, “Automated Screening for Diabetic Retinopathy - A Systematic Review,” Ophthalmic Res. **60**(1), 9–17 (2018). [CrossRef]

**50. **N. Salamat, M. M. S. Missen, and A. Rashid, “Diabetic retinopathy techniques in retinal images: A review,” Artif. Intell. Med. **97**, 168–188 (2019). [CrossRef]

**51. **Y. Aso, T. Inukai, K. Tayama, and Y. Takemura, “Serum concentrations of advanced glycation endproducts are associated with the development of atherosclerosis as well as diabetic microangiopathy in patients with type 2 diabetes,” Acta Diabetol. **37**(2), 87–92 (2000). [CrossRef]

**52. **R. Klein, A. R. Sharrett, B. E. Klein, S. E. Moss, A. R. Folsom, T. Y. Wong, F. L. Brancati, L. D. Hubbard, and D. Couper, “The association of atherosclerosis, vascular risk factors, and retinopathy in adults with diabetes : the atherosclerosis risk in communities study,” Ophthalmology **109**(7), 1225–1234 (2002). [CrossRef]

**53. **T. L. Torp, R. Kawasaki, T. Y. Wong, T. Peto, and J. Grauslund, “Temporal changes in retinal vascular parameters associated with successful panretinal photocoagulation in proliferative diabetic retinopathy: A prospective clinical interventional study,” Acta Ophthalmol. **96**(4), 405–410 (2018). [CrossRef]

**54. **J. Lei, E. Yi, Y. Suo, C. Chen, X. Xu, W. Ding, N. S. Abdelfattah, X. Fan, and H. Lu, “Distinctive Analysis of Macular Superficial Capillaries and Large Vessels Using Optical Coherence Tomographic Angiography in Healthy and Diabetic Eyes,” Invest. Ophthalmol. Visual Sci. **59**(5), 1937–1943 (2018). [CrossRef]

**55. **J. Sahlsten, J. Jaskari, J. Kivinen, L. Turunen, E. Jaanio, K. Hietala, and K. Kaski, “Deep Learning Fundus Image Analysis for Diabetic Retinopathy and Macular Edema Grading,” Sci. Rep. **9**(1), 10750 (2019). [CrossRef]

**56. **J. Krause, V. Gulshan, E. Rahimy, P. Karth, K. Widner, G. S. Corrado, L. Peng, and D. R. Webster, “Grader Variability and the Importance of Reference Standards for Evaluating Machine Learning Models for Diabetic Retinopathy,” Ophthalmology **125**(8), 1264–1272 (2018). [CrossRef]

**57. **U. R. Acharya, C. M. Lim, E. Y. Ng, C. Chee, and T. Tamura, “Computer-based detection of diabetes retinopathy stages using digital fundus images,” Proc. Inst. Mech. Eng., Part H **223**(5), 545–553 (2009). [CrossRef]

**58. **E. D. Cole, E. A. Novais, R. N. Louzada, and N. K. Waheed, “Contemporary retinal imaging techniques in diabetic retinopathy: a review,” Clin. Exp. Ophthalmol. **44**(4), 289–299 (2016). [CrossRef]

**59. **T. E. de Carlo, A. Romano, N. K. Waheed, and J. S. Duker, “A review of optical coherence tomography angiography (OCTA),” Int. J. Retin. Vitr. **1**(1), 5 (2015). [CrossRef]

**60. **J. Lee and R. Rosen, “Optical Coherence Tomography Angiography in Diabetes,” Curr. Diabetes Rep. **16**(12), 123 (2016). [CrossRef]

**61. **T. Y. Wong, C. M. G. Cheung, M. Larsen, S. Sharma, and R. Simó, “Diabetic retinopathy,” Nat. Rev. Dis. Primers **2**(1), 16012 (2016). [CrossRef]

**62. **E. J. Topol, “High-performance medicine: the convergence of human and artificial intelligence,” Nat. Med. **25**(1), 44–56 (2019). [CrossRef]

**63. **G. C. Cawley and N. L. Talbot, “Fast exact leave-one-out cross-validation of sparse least-squares support vector machines,” Neural Netw **17**(10), 1467–1475 (2004). [CrossRef]