Classification of advanced and early stages of diabetic retinopathy from non-diabetic subjects by an ordinary least squares modeling method applied to OCTA images

: 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 quantiﬁed from OCTA images and were used to determine predicted probabilities of PDR, no/mild NPDR, and non-diabetics. The classiﬁcation 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 classiﬁcation of stages of DR and other ocular and systemic diseases. America under the terms of the OSA Open Access Publishing Agreement


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][9][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][12][13][14][15]. Discrimination of DR images has been made possible by texture classification [16,17] and lesion detection in retinal images [18][19][20].
Furthermore, support vector machine classification methods have effectively detected DR using retinal images in recent studies [21][22][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.

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][31][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.

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.

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.

Vessel density and spacing
A previously validated local fractal dimension analysis method was applied to OCTA images to assess retinal vessel density [35][36][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.

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][28][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].

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 The finite sample form of Eq. (3) employed to estimate β is: In Eq. (4) In our previous method for discriminating groups based on Fisher's linear discriminant (FLD) method [27][28][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 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.

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 (π i = 1) is in agreement with the logistic model [44]. For a set of data to be classified into 2 classes, π 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].

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].  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).

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.

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.  a OLS = ordinary least squares; b PDR = proliferative diabetic retinopathy; c LR = Likelihood Ratio; d R 2 = Nagelkerke coefficient of determination (ranging from 0 to 1); e AUROC = area under receiver operating characteristic curve. Results of individual and multiple OLS parameters are derived from univariate and binary multivariate logistic regression models, respectively. Individual OLS parameters without significant Wald X 2 and LR X 2 are not displayed. f P < 0.05; g P ≤ 0.005; h P > 0.20.

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).

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 nondiabetics. 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][59][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).