Title Referable Diabetic Retinopathy Prediction Algorithm Applied to a Population of 120389 Type 2 Diabetics over 11 Years Follow-Up.

1. Introduction

Diabetes mellitus (DM) is a chronic disease differentiated into two types: type 1 patients (T1DM) (young patients with an initial onset that requires treatment with insulin and possibly associated with the group of autoimmune diseases with a possible viral trigger) and type 2 patients (T2DM), which is the more frequent form and associated with insulin resistance problems secondary to obesity and metabolic syndrome. DM is affects 463 million people globally, according to 2019 data from the International Diabetes Federation (IDF). The IDF also estimates that T2DM will affect 700 million people by 2045 [1].

In Spain, a study conducted by di@bet.es has reported that 7.9% of the population are known to have DM, although it estimates true prevalence to be double that, at 13.8% [2]. T2DM affects various organs of the body in the form of vascular involvement, such as stroke and myocardial infarction (more frequent conditions in patients with DM, and in the kidneys and eyes, affecting the small vessels in the form of nephropathy or diabetic retinopathies. Currently, the retina is considered to develop a neurovascular lesion by affecting both the retinal vessels and its neuronal cells [3].

The scientific societies that focus on diabetes and ophthalmology recommend screening by fundus photographs (retinographs) every 1 to 2 years, depending on a patient's risk of progression to DR [3,4,5].

Most current screening programmes, however, are failing to provide annual fundus monitoring [7,8]. Given the real difficulties of applying these clinical guidelines, the same scientific societies recommend at least biannual control for those patients who have a good control of DM but there is no consensus on what good control actually means.

To address this, we have developed a diabetic retinopathy risk algorithm [9] that we are piloting in our clinical practice [10]. The present study evaluates the effectiveness of that algorithm in a population of 120,389 patients with T2DM whom we have followed over 11 years from 2010 to 2021, and for whom we have the necessary clinical data.

2. Materials and Methods

2.1. Setting. Spain

Database from Catalonia with 120,382 DM2 patients. We provided more information in a previous publication [11]

2.2. Sample Size

We included retrospective data from the Electronic Health Records (EHR) of 120,382 DM2 patients, provided by SIDIAP (“Sistema d'informació pel Desenvolupament de la Recerca a Atenció Primària”) [12]. The SIDIAP database contains a variety of patient data, such as visits to healthcare professionals, diagnostic codes, demographic information, clinical variables, laboratory tests results, prescriptions, referrals to specialists and hospitals, and medication dispensed in pharmacies. For this analysis, data was extracted from an 11-year period (2010-2020 inclusive).

2.3. Inclusion Criteria

• Patients with type 2 DM.
• Patients without DR or with mild-DR.

2.4. Exclusion Criteria

• Patients with type 1 DM.
• Patients included in Diabetes group III and other specific types (i.e. Diseases of the exocrine pancreas, Endocrinopathy, Genetic defects of ß-cell function, Genetic defects in insulin action).
• Patients included in Diabetes group IV and gestational diabetes mellitus (GDM)
• Patients who did not have a complete EHR.
• Patients with DR higher than mild

2.5. Construction of the Algorithm

The diabetic retinopathy prediction algorithm (DRPA) has been in development since 2010 in patients with T2DM. A total of 16 risk factors were considered: Current age, age at diagnosis of DM, sex, type of DM, body mass index, duration of DM, DM treatment, smoker, control of high blood pressure, diastolic pressure rate , systolic blood pressure rate, HbA1c%, creatinine, estimated glomerular filtration rate (eGFR) measured by CKD-EPI, total cholesterol, LDL cholesterol, HDL cholesterol, triglycerides and microalbuminuria [13,14].Of all those risk factors, the decision tree statistical system with fuzzy rules was used to generate the first two DRPAs (RETIPROGRAM 1.0 and 2.0). An improved version 3.0, which has been piloted at Hospital Universitario Sant Joan [10] includes the following factors:

• Sex.
• Body mass index.
• Duration of T2DM in units of one year.
• T2DM treatment, diet, oral antidiabetics, insulin, insulin analogues.
• Control of arterial hypertension normal value systolic BP <140 diastolic BP <90.
• HbA1c% in 1% fractions
• Estimated glomerular filtration rate, calculated from plasma creatinine using the chronic kidney disease epidemiology collaboration equation (CKD-EPI equation).
• Microalbuminuria value 30 mg/min up to 300 mg/min

We have recently added another variable: 'mild-DR (yes or no)' to enable us to assess which patients might develop more serious forms such as moderate-DR, severe-DR or diabetic macular edema. This has now developed into RETIPROGRAM 4.0, which is the latest version. The objective of the present study is to evaluate the validity of RETIPROGRAM 4.0.

The algorithm was based on a fuzzy random forest, which is a set of decision trees with fuzzy rules. An algorithm was built with 100 trees in the forest and three variables in each node. The statistical results of the test set obtained an AUC value of 0.807%, sensitivity of 80.67%, and specificity of 85.96% [15].

2.7. Statistical Methods. Data was analysed using SPSS, version 22.0 (IBM® Statistics, Chicago, IL, USA). Descriptive statistical analysis of quantitative data was made by the determination of mean and standard deviation. For qualitative data we used the analysis of frequency and percentage in each category. The univariate study was carried out using the two sample Student T-tests to compare quantitative variables, and we used the chi-squared table for qualitative data. A p value of less than 0.05 was considered statistically significant.

We measured the screening performance of the study using a confusion matrix/contingency 2x2 table. Given a classified data set, there were four basic combinations of actual and assigned: correct positive assignments, or true positives [TP], correct negative assignments, or true negatives [TN], incorrect positive assignments, or false positives [FP], and incorrect negative assignments, or false negatives [FN].

The statistical evaluation of the data set included: accuracy (A), sensitivity [S], specificity [SP], positive predictive value or precision [PPV], Harmonic means (F1 score), negative predictive value [NPV], positive false discovery rate or type 1 error [α], negative false discovery rate or type 2 error [β], and the area under the curve or diagnostic effectiveness expressed as a proportion of correctly classified subjects.

Sensitivity or recall = TP / (TP+FN). These are the proportion of the population with the condition for which the test is correct. Specificity = FP / (FP + TN). In clinical diagnosis, the values of sensitivity and specificity are considered good when they exceed 80%.

Predictive positive value or precision = TP / (TP + FP). These are the proportion of the population with a given test result for which the test is correct. Negative predictive value = FN / (FN + TN).

Harmonic means or F1 Score = 2xprecisionxrecall/ precision+recall, also known as the Sørensen–Dice coefficient or Dice similarity coefficient (DSC) in statistical analysis of binary classification, the F1 score (also F-score or F-measure) is a measure of a test's accuracy, calculated from the positive predictive value and sensitivity of the test. It is a measure of performance of the test that combines precision predictive positive value) and recall (sensitivity) into one number via a choice of weighing, simply equal weighing. The highest possible value of F1 is 1, indicating perfect precision and recall, and the lowest possible value is 0, if either the precision or the recall is zero.

Also, we determined the Accuracy or diagnostic effectiveness = TP+TN / (TP+FP+TN+FN), which is expressed as a proportion of correctly classified subjects among all of them. The accuracy is affected by prevalence, with the same sensitivity and specificity, diagnostic accuracy of a particular test increases as the disease prevalence decreases. On the other hand, the area under the curve (AUC) is not affected by the prevalence and allows the performance of the classifier to be represented as a single value, in addition to allowing two or more different classifiers to be compared. Values of accuracy and AUC of >0.8 is almost perfect correlation, 0.6 - 0.79 is significant correlation, 0.4 – 0.59 is moderate correlation, 0.2 – 0.39 is regular correlation, and 0-0.2 is poor correlation.

3. Results

3.1. Demographic Data

Table 1 gives the demographic values of the sample size. The sample in the present study was 120,389 patients, of whom 68578 were men (57%) and 51811 were women (43%).

Table 2 shows patient status of DR at the beginning and at the end of the study. A total of 111,172 patients with no DR or with only mild-DR were recruited. During 11 years of follow up, 18,694 of those patients had been diagnosed with DR (15.5%), of which 5775 patients (4.77%) had developed RDR (DR moderate or plus), and 1581 patients (1.31%) had developed sight-threatening DR (severe-DR, proliferative-DR and diabetic macular edema).

3.2. Statistical Analysis of Confusion Matrix/Contingency

We can see in Table 3 that the algorithm has both sensitivity and specificity within the limits considered necessary by the British Diabetes Association [28]; sensitivity is above 80% both in the case of predicting any type of DR and also for predicting RDR. Furthermore, the metrics that determine whether the predicted cases are correctly classified, either positively or negatively, we can see that the precision (PPV) is 78% for any-DR and 76% for RDR, which is a significant correlation. If we look at the harmonic mean (HM, or F1 score), which combines both precision and recall (sensitivity), the values are significant at 83% for any-DR and 79% for the case of RDR. Regarding accuracy, it is identical in both cases, the value being 97% for both any-DR and RDR. Accuracy is affected by the prevalence of the disease. In this case, prevalence of DR is 15.5% after 11 years of follow-up, which is a somewhat low figure. In these cases, accuracy tends to have high values, such as the are in the present study.

On the contrary, the AUC is not affected by prevalence and indicates whether the positive and negative predictions in our study are balanced. The result differs for both types of DR, so in order to predict any-DR the AUC gets a value of 93% higher than the prediction of RDR which is 90%

Note: For any-DR we have taken into account that 9477 patients developed some form of DR during the study, and that 5775 patients developed referable -DR during the study.

Table 3. Statistical confusion matrix of algorithm in our sample size.

	Any DR	RDR
True Positive	8387	4727
False Positive	2324	1466
True Negative	108588	113148
False Negative	1090	1048
Accuracy	0.97 (95% CI, 0.96–0.98)	0.97 (95% CI, 0.95–0.99)
AUC (area under de curve ROC)	0.93 (95% CI, 0.92–0.94)	0.90 (95% CI, 0.89–0.91)
Sensitivity or recall	0.88 (95% CI, 0.86–0.90)	0.82 (95% CI, 0.80–0.84)
Specificity	0.98 (95% CI, 0.96–0.99)	0.99 (95% CI, 0.95–0.994)
HM or F1 score	0.83 (95% CI, 0.81–0.84)	0.79 (95% CI, 0.78–0.80)
Precision or Positive predictive values	0.78 (95% CI, 0.75-0.80)	0.76 (95% CI, 0.74-0.80)
Negative predictive values	0.99 (95% CI, 0.98-0.999)	0.99 (95% CI, 0.97-0.997)

Note: For any‐DR we have taken into account that 9477 patients developed some form of DR during the study, and that 5775 patients developed referable ‐DR during the study.

Table 3. Statistical confusion matrix of algorithm in our sample size.

	Any DR	RDR
True Positive	8387	4727
False Positive	2324	1466
True Negative	108588	113148
False Negative	1090	1048
Accuracy	0.97 (95% CI, 0.96–0.98)	0.97 (95% CI, 0.95–0.99)
AUC (area under de curve ROC)	0.93 (95% CI, 0.92–0.94)	0.90 (95% CI, 0.89–0.91)
Sensitivity or recall	0.88 (95% CI, 0.86–0.90)	0.82 (95% CI, 0.80–0.84)
Specificity	0.98 (95% CI, 0.96–0.99)	0.99 (95% CI, 0.95–0.994)
HM or F1 score	0.83 (95% CI, 0.81–0.84)	0.79 (95% CI, 0.78–0.80)
Precision or Positive predictive values	0.78 (95% CI, 0.75-0.80)	0.76 (95% CI, 0.74-0.80)
Negative predictive values	0.99 (95% CI, 0.98-0.999)	0.99 (95% CI, 0.97-0.997)

Note: For any‐DR we have taken into account that 9477 patients developed some form of DR during the study, and that 5775 patients developed referable ‐DR during the study.

4. Discussion

The objective of the present study was to assess whether the algorithm that we had previously developed, and to which we had incorporated the presence of mild-DR as a variable, was useful for predicting the evolution of patients to more advanced forms of DR. In our case it was. RDR includes moderate-DR, severe-DR, proliferative-DR and diabetic macualr edema forms. All of these forms affect eyesight so predicting them will allow early intervention. In those patients with the mildest form (mild-DR) adjusting the clinical variables (HbA1c, kidney function, blood pressure) will slow down the progression to forms that might affect vision.

We had proposed two phases. The first was to identify those patients who are likely to develop DR but who did not have DR at the beginning of the study, which is why we left out the patients with mild-DR for this first phase. The results are interesting, having an AUC of 93% and a harmonic mean of 83%, with a precision of 78% and a recall (sensitivity) of 88%, figures that indicate that the algorithm is useful in these cases. More importantly, specificity is 98% and the negative predictive values are 99% meaning that the algorithm predicts with sufficient certainty those patients who, due to their clinical characteristics, are not going to develop DR. By linking both positive and negative predictors, it allows us to identify those patients at low risk of developing DR and who would be safe being scanned less frequently, thus freeing up resources for targeting the more urgent cases.

A second phase focused on those patients who had developed advanced forms of DR (RDR) during the study period. Results are somewhat lower than those previously mentioned but it is important to be able to draw conclusions. In this second phase of the study, the AUC values were 90% with a precision of 76% and a harmonic mean of 79%, with the recall being 82%. These values are lower than those for the detection of any -DR but they are significant enough to be useful for predicting the development of RDR based on the clinical data of the algorithm and the presence or absence of mild-DR as an added risk variable. In this case, however, it is very important that both the specificity and the negative predictive values are 99%; when the algorithm predicts that a patient is at low risk of developing RDR it is quite safe.

The two phases of the study, predicting the development of both any-DR and RDR, meet the criteria defined by the British Diabetic Association, which proposed that diabetic retinopathy screening programmes should satisfy a minimum of 80% sensitivity and specificity, values that the present study has achieved.

When comparing our results with previous studies, it is important to note that only three authors have published work on the development of prediction algorithms, Aspleund [16], Scanlon [17] and Broadbent [18]. Table 4 shows that the AUC of the present study is higher than those obtained by the other authors.

Table 4 shows that the most widely applied algorithm has been developed by Aspelund [16]. It was was constructed from the Diabetic Retinopathy database managed at the Ophthalmology Department of Aarhus University Hospital in Denmark. The database has been built from the clinical data of 5199 patients over 20 years, thus allowing the algorithm to be tested prospectively. The model provides a recommended interval for follow-up fundus screening for the presence of retinopathy affecting vision (STDR) of between 6 and 60 months. The algorithm has been further tested in other countries such as Spain, where Soto Pedre [19] studied 508 patients with T2DM. Results showed that 3.1% developed STDR before their following screening visit, with the value of the area under the curve (AUC) being 0.74. In the Netherlands, Van der Heijden [20] tested the algorithm on 76 patients for a mean of 53 months, and reported an AUC value of 0.83. Finally, in the UK, Lund [21], with a sample of 9690 DM patients followed for 2 years, reported that the algorithm predicted the onset of DR stages with AUC values of 0.833 for T2DM patients.

The second algorithm that has published results is the model developed by Scanlon et al in Gloucestershire, United Kingdom, [17]. This algorithm was validated with 15877 patients, obtaining an AUC of 0.77 for predicting the development of STDR.

Finally, the model developed by Broadbent [18], also known as the Liverpool Risk Calculation Engine, was built to detect the risk of developing STDR. The statistical analysis of this algorithm found an AUC value of 0.88 in the prediction of STDR, with sensitivity values of 0.61 and specificity of 0.93. We should highlight that it is the only report that provides sensitivity and specificity data. This algorithm was tested with a population of DM patients from Liverpool (UK) and has been published in the form of an action protocol under the name of the ISDR protocol [22].

Our algorithm is different in that it predicts RDR and any-DR rather than STDR in the other three. Also, Broadbent uses 7 variables, Aspelund uses 5 and Scanlon only 3, whereas our algorithm uses 10. Only our algorithm and that of Broadbent take into account whether there is DR at the beginning of the study or not as a variable for predicting STD, which suggests that using more variables might improve the robustness of the results, since Broadbent's algorithm has an AUC of 0.88 similar to ours of 0.90 for predicting STDR.

Table 5. Risk factors used by each algorithm.

	Aspelund	Scanlon	Broadbent	Authors
Current age	√	√		√
Age at diagnosis			√
Sex			√	√
DM duration	√		√	√
DM treatment				√
Systolic blood pressure	√		√	√
Diastolic blood pressure				√
Total-Cholesterol		√	√
HbA1c %	√	√	√	√
Microalbuminuria				√
Glomerular filtration rate measured by CKD EPI				√
Body mass index				√
DM Type	√
Diabetic retinopathy			√	√

Table 5. Risk factors used by each algorithm.

	Aspelund	Scanlon	Broadbent	Authors
Current age	√	√		√
Age at diagnosis			√
Sex			√	√
DM duration	√		√	√
DM treatment				√
Systolic blood pressure	√		√	√
Diastolic blood pressure				√
Total-Cholesterol		√	√
HbA1c %	√	√	√	√
Microalbuminuria				√
Glomerular filtration rate measured by CKD EPI				√
Body mass index				√
DM Type	√
Diabetic retinopathy			√	√

The limitation of the present study is that it is retrospective, meaning that patients who do not have all the variables collected are excluded, which can generate bias in data collection. Likewise, we have had to rely on the data collected in Electronic Health Records, meaning that we have to be certain that the data has been collected correctly. To avoid this bias, we have the reliability of SIDIAP in providing the veracity in the results collected.

5. Conclusions

We have developed an algorithm that reliably predicts the development of any-DR in type 2 diabetic patients and the development of referable forms of DR with a high degree of certainty. Additionally, we can identify those patients who are not going to develop DR by applying our algorithm, which will allow us to extend intervals for screening for diabetic retinopathy with greater security. Let us not forget that one of the particular difficulties of DR screening is lack of attendance, often because of annual frequency of DR screening. In these cases, our longer intervals will encourage the most reluctant of patients to attend.