# Automated multimodal severity assessment of diabetic retinopathy using ultra-widefield color fundus photography and clinical tabular data

Alireza Rezaei, Sarah Matta, Rachid Zeghlache, Pierre-Henri Conze, Capucine Lepicard, Pierre Deman, Laurent Borderie, Deborah Cosette, Sophie Bonnin, Aude Couturier, et al.

► To cite this version:

Alireza Rezaei, Sarah Matta, Rachid Zeghlache, Pierre-Henri Conze, Capucine Lepicard, et al. Automated multimodal severity assessment of diabetic retinopathy using ultra-widefield color fundus photography and clinical tabular data. Biomedical Signal Processing and Control, 2026, 112 (Part C), pp.108673. 10.1016/j.bspc.2025.108673. hal-05266495

HAL Id: hal-05266495

https://imt-atlantique.hal.science/hal-05266495v1

Submitted on 28 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Automated multimodal severity assessment of diabetic retinopathy using ultra-widefield color fundus photography and clinical tabular data

Alireza Rezaei<sup>a,b,c</sup>,*, Sarah Matta<sup>a,b</sup>, Rachid Zeghlache<sup>a,b</sup>, Pierre-Henri Conze<sup>a,c</sup>, Capucine Lepicard<sup>e</sup>, Pierre Deman<sup>g,h</sup>, Laurent Borderie<sup>h</sup>, Deborah Cosette<sup>i</sup>, Sophie Bonnin<sup>e,f</sup>, Aude Couturier<sup>e,f</sup>, Béatrice Cochener<sup>a,b,d</sup>, Mathieu Lamard<sup>a,b</sup>, Mostafa El Habib Daho<sup>a,b</sup>,<sup>i</sup>, Gwenolé Quellec<sup>a,1</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Diabetic retinopathy classification

Multimodal information fusion

Deep learning

UWF-CFP

Tabular data

## ABSTRACT

This study introduces an automatic deep-learning-based approach to diabetic retinopathy (DR) severity assessment by integrating two modalities: Ultra-Widefield Color Fundus Photography (UWF-CFP) from the CLARUS 500 device (Carl Zeiss Meditec Inc., Dublin, CA, USA) and a comprehensive set of clinical data from the EVIRED project. We propose a framework that combines the information from 2D UWF-CFP images and a set of 76 tabular features, including demographic, biochemical, and clinical parameters, to enhance the classification accuracy of DR stages. Our model uses advanced machine learning techniques to address the complexities of synthesizing heterogeneous data types, providing a holistic view of patient health status. Results indicate that this fusion outperforms traditional methods that rely solely on imaging or clinical data, suggesting a robust model which can provide practitioners with a supportive second opinion on DR severity, particularly useful in screening workflows. We measured a multiclass accuracy of 63.4% and kappa of 0.807 for our fusion model which is 2.1% higher in accuracy and 0.022 higher in kappa compared to the image unimodal classifier. Several interpretation methods are used to provide practitioners with an inside view of the workings of classification methods and allow them to discover the most important clinical features.

## 1. Introduction

Diabetic Retinopathy (DR), a prevalent ocular complication of diabetes, is a leading cause of blindness worldwide [1]. This disease is characterized by progressive retinal damage due to prolonged hyperglycemia, affecting roughly one-third of individuals with diabetes. DR severity is defined by the International Clinical Diabetic Retinopathy Severity Scale (ICDR), and includes the following grades: no DR, mild, moderate or severe non-proliferative DR (NPDR), and proliferative DR (PDR). The timely and accurate diagnosis of DR is essential for effective management and treatment, yet the subtle and complex nature of early-stage DR poses significant diagnostic challenges.

In the realm of medical image analysis, deep learning models have shown remarkable success. Convolutional Neural Networks (CNNs) are particularly effective for retinal image segmentation and classification. A CNN backbone was used by Maninis et al. [2] to perform segmentation for blood vessels and the optic disk in the same network. Similarly, El-Hag et al. [3] uses CNN to detect exudates in the retina and classify images into normal and abnormal cases of DR. Disease screening frameworks use CNNs to classify patients' conditions into different grades. Ting et al. [4] compared the performance of CNNs with expert graders in detecting cases of Moderate DR or worse and

<footnote><sup>a</sup> LaTIM UMR 1101, Inserm, Brest, France</footnote>

<footnote><sup>b</sup> University of Western Brittany, Brest, France</footnote>

<footnote><sup>c</sup> IMT Atlantique, Brest, France</footnote>

<footnote><sup>d</sup> Ophtalmology Department, CHRU Brest, Brest, France</footnote>

<footnote><sup>e</sup> Lariboisière Hospital, AP-HP, Paris, France</footnote>

<footnote><sup>f</sup> Hôpital Fondation Adolphe de Rothschild, Paris, France</footnote>

<footnote><sup>g</sup> ADCIS, Saint-Contest, France</footnote>

<footnote><sup>h</sup> Evolucare Technologies, Le Pecq, France</footnote>

<footnote><sup>i</sup> Carl Zeiss Meditec, Dublin, CA, USA</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: alireza.rezaei@univ-brest.fr (A. Rezaei).</footnote>

<footnote n="1">¹ Contributed equally as senior authors to this work.</footnote>
reported an average AUC score of 0.936; however they did not attempt classification of samples into specific DR severity grades. Ho et al. [5] examined five different CNN architectures for the screening process of 28 different ocular pathologies and reported an AUC score of 0.98 for the detection of DR; however, they did not perform DR severity grading.

These advances in medical image processing, along with innovations in imaging techniques, have greatly improved the detection and classification of DR. Ultra-Widefield Color Fundus Photography (UWF-CFP) provides a comprehensive panoramic view of the retina, enhancing the assessment of retinal pathology [6,7]. While the traditional color fundus images can cover around a third of the retina, an ultra-widefield image can show up to 82% of the retina and capture peripheral lesions [7].

Although these imaging techniques are effective, they offer limited information on the systemic factors that influence DR progression (diabetes types, diabetes duration, etc.). Integration of imaging data with various clinical features, recorded as tabular data, offers a promising avenue to improve diagnostic accuracy [8]. Traditionally, tabular data classification has been dominated by classical ML models such as Decision Trees [9]. However, recent years have seen the emergence of novel approaches that apply deep learning (DL) to tabular data. For instance, Gorishniy et al. [10] introduced FTTransformer, an adaptation of the transformer architecture for tabular data, which demonstrated improvements over existing algorithms. Another example is GANDALF [11], which features a representation learning block based on a gating mechanism [12]. To address the substantial data requirements of neural network models, Hollmann et al. [13,14] proposed TabPFN, a prior-data fitted transformer capable of solving classification problems with small datasets in a single forward pass.

The integration of imaging and tabular data enables the analysis of heterogeneous data types to provide a more holistic view of patient health [15]. Previous studies in medical imaging have utilized a limited number of tabular features to improve diagnosis. For example, patient's age has been used by Artzi et al. [16] in conjunction with MRI images in a CNN architecture to classify fossa tumors in children. They report an increase of 10% overall accuracy compared to an image only classifier. In a similar way, Wolf et al. [17] integrates a set of 9 tabular features, including biomarkers and demographic variables, with 3D images using a ResNet CNN structure to perform various tasks (diagnosis and progression) related to Alzheimer's Disease. They report an increase of 7% and 12% accuracy when compared to a tabular-only and image-only classifier respectively. In ophthalmology, Hsu et al. [18] used a set of 22 clinical and demographic tabular features alongside fundus color images to detect cases with *moderate or worse DR*. They reported an increase of 2.82% accuracy when compared to an image-only model. Although they demonstrated the potential of tabular data in detecting worse cases of DR, they did not perform multi-class grading for DR severity. In conclusion, to the best of our knowledge, the current work is the first study that jointly analyzes UWF-CFP and a large set of contextual data for multi-class grading of DR.

In this study, we provide a multimodal integration for the severity assessment of diabetic retinopathy. We leverage high-resolution UWF-CFP images alongside an extensive set of clinical features from the EVIRED<fnref n="2" /> project, encompassing 76 distinct clinical and biochemical markers. The novelty of our approach lies in multi-class classification of DR severity with the use of ultra-widefield imagery as well as a more expansive set of clinical measurements compared to previous studies. In addition, we explore both classic ML and deep learning methods for the process of tabular data and investigate multiple fusion combinations between image and tabular classifiers. Our aim is to capitalize on the complementary strengths of imaging and clinical data to improve understanding of DR severity and ultimately enhance diagnostic accuracy and patient outcomes. Finally, We replicate our experiments on two new external multimodal datasets (BRset and mBRset) and we show that the results are consistent with our findings on the EVIRED dataset.

## 2. Methods

In this section, we present our approach to DR classification using tabular data, imaging data, and a combination of both. For each modality, we explore classification models as well as the necessary data preprocessing steps.

### 2.1. Tabular classification

We perform tabular data classification using both deep learning and classic ML models. A model sweep is performed on the following models: Random Forest [19], FTTransformer [10], GANDALF [11], and TabPFN [13]. In addition, we include two simple architectures as described in Gorishniy et al. [10]:

* **A Multi-Layered Perceptron** with a single linear, activation, and dropout block:

$$MLP(x) = Linear(Dropout(ReLU(Linear(x)))) \quad (1)$$

* **A ResNet** with one residual block:

$$\begin{aligned} ResNet(x) &= ResBlock(Linear(x)) \\ ResBlock(x) &= x + Dropout(Linear( \\ & \quad Dropout(ReLU(Linear( \\ & \quad BatchNorm(x))))) \end{aligned} \quad (2)$$

In our experiments, we set the probability of all dropout layers to 0.4.

Most tabular models require the input data to follow certain conditions to work properly. Missing data, class imbalance, data scaling, and categorical feature encoding are some of the challenges that need to be addressed. To solve the problem of missing data, we follow a few different complementary approaches:

* **Filling with patient's other data:** If the patient has previous recorded values for the missing feature, the last recorded feature value for that patient is used (also known as LOCF: last observation carried forward).

* **Filling with global average:** If a patient does not have a previously recorded value for a numerical feature, the missing value is filled with a global median across the dataset. We have used median and not mean to ignore outlier values. Examples include blood pressure, height, weight, etc.

* **Filling with zero:** for binary features which are mostly zero in the dataset. Examples include whether the patient has undergone dialysis or received a kidney transplant.

* **Encoding as a category:** for categorical features, missing values are encoded as a new category. Examples include smoking habits or a particular treatment.

* **Dropping the feature:** If a feature has more than 90% missing values, it is dropped entirely from the dataset. In our dataset, the recorded visual acuity in ETDRS units was dropped for this reason.

Class imbalance is dealt with using weighted random sampling during the training process. Sampling weight $w_{c_i}$ for class $c_i$ is defined as

$$w_{c_i} = N_{train}/N_{c_i}, i \in \{0, 1, 2, 3, 4, 5\} \quad (3)$$

where $N_{train}$ is the total number of samples in the training set, and $N_{c_i}$ is the number of training samples in class $c_i$. Samples with a higher $w_{c_i}$ are more likely to be picked by the sampler during training. In general, other approaches such as random over-sampling, random under-sampling, and SMOTE [20] may also prove useful to combat class imbalance. However, for our dataset, initial experiments demonstrated superior performance on the validation set using the weighted random sampling strategy. Pure undersampling methods resulted in loss of

<footnote n="2">² https://evired.org.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. General fusion model Architecture.</figcaption></figure>

information in most classes and lead to lower accuracy. On the other hand, simple oversampling methods did not remove model's bias on the majority class and caused complete failure of the training process.

The features are divided into two groups: categorical and numerical. Numerical features are scaled to unit variance with a mean of zero: sample x is scaled to x':

$$x' = (x - u)/s, \qquad (4)$$

where u and s are the mean and standard deviation of the samples, respectively.

Categorical features (which include binary features) are encoded using ordinal encoding, resulting in integer feature values from 0 to K − 1, where K is the number of categories for that feature. Missing values are encoded as − 1.

## 2.2. Image classification

In the realm of medical image analysis, deep learning models, particularly CNNs, have demonstrated remarkable success. They are adept at learning complex hierarchical representations from images, making them a natural choice for integrating multimodal data [21–24]. In this work, each image I is randomly cropped with different scales and aspect ratios and then resized to 448 × 448 pixels. Random horizontal flips and color jitter for saturation, brightness, and contrast are used for augmentation. Finally, the image is normalized to be used in the classification model.

Similar to tabular classification, weighted random sampling (Eq. (3)) is used to alleviate the effects of class imbalance.

## 2.3. Multimodal fusion and classification

Multimodal deep learning approaches that combine imaging data with other data types, such as electronic health records (EHR), have shown significant potential in improving diagnostic performance [15].

Fusion can be applied at various stages: input fusion, intermediate fusion, or output fusion. In input fusion, data from different modalities are simply concatenated, though this approach lacks geometric relevance in this context. Output fusion, which combines independent class predictions from each modality, is generally less efficient than intermediate fusion, where predictions are merged earlier in the decision-making process. In this work, we thus combine imaging and tabular modalities through an intermediate fusion of feature vectors. Different tabular models produce feature vectors of varying sizes M. FTTransformer provide a feature vector of size M = 32; for GANDALD M = 139; and for TabPFN M = 192. Similarly, for MLP and Residual models we set M = 64. Finally, for Random Forest, since there are no intermediate features, we use its output of size 6. To perform the fusion, we simply reduce the size of the image feature vector to match M.

The overall architecture is depicted in Fig. 1. The model consists of two main branches: one for the image features extraction and the other for the tabular data. The outputs of these branches are concatenated and further processed to produce the final classification. Both feature vectors go through fully connected layers with ReLU activation before concatenation and the classification layer. Dropout layers are used to avoid over-fitting.

For image feature extraction, we used different pretrained ResNet models (following the results of El Habib Daho et al. [8]), denoted by F(). The obtained features, f$_{image}$ ∈ R$^{2048}$, are passed through a Fully Connected (FC) layer with ReLU activation and Dropout regularization:

$$f_{\text{image}} = F(I),$$

$$h_{\text{image}} = \text{Dropout}(\text{ReLU}(W_1 f_{\text{image}} + b_1)),$$

where W$_{1}$ ∈ R$^{M×2048}$ and b$_{1}$ ∈ R$^{M}$. In our experiments, M = 76 and dropout layers have a probability of 0.4.

# 3. Experiments

## 3.1. Data and training details

In this section, we describe the datasets used for our experiments and give details on the training process, hyperparameters, and the hardware.

### 3.1.1. The EVIRED dataset

The data utilized in this research were sourced from the ongoing EVIRED project, a collaborative effort involving 14 hospitals across France. This rich dataset integrates two fundamental elements for the study: high-resolution UWF-CFP images and an extensive set of 76 clinical features.

The UWF-CFP component includes images captured using the CLARUS 500 system from Zeiss, renowned for its high-resolution and wide-field imaging capabilities. The dataset contains around 3200 high-quality images from 1031 patients diagnosed with various stages of DR, ranging from no apparent retinopathy to proliferative DR. The images offer a panoramic view of the retina, significantly exceeding the standard field of view provided by traditional fundus photography, thus allowing for a more comprehensive assessment of peripheral retinal lesions critical for DR diagnosis. Each image varies in size from 3900 × 3900 to 7900 × 4900 pixels, reflecting the image stitching process employed to create a complete retinal montage.

The clinical data comprises 76 detailed tabular features for each patient, capturing demographic information (e.g., age, gender), diabetes duration, treatment specifics, a self-reported health questionnaire, and a comprehensive panel of biochemical markers, including glycated hemoglobin (HbA1c), fasting blood glucose, lipid profiles, and kidney
function tests. These features were meticulously collected and curated, reflecting a broad spectrum of systemic health parameters that could influence DR progression. Table 1 reports different types of measured features and some examples of each category. A complete list is available in Table A.12 in Appendix. Integrating these clinical features with imaging data aims to create a holistic model that not only assesses retinal changes but also considers systemic health variables indicative of microvascular complications in diabetes.

The dataset entries were labeled based on the International Clinical Diabetic Retinopathy Severity Scale (ICDR), comprising the following grades: no DR, mild, moderate or severe non-proliferative DR (NPDR), proliferative DR (PDR), to which we add an additional class indicating that the patient received panretinal photocoagulation treatment (PRP). Each image was independently graded by two orthoptists and then confirmed by an experienced ophthalmologist, ensuring high-quality, accurate annotations.

### 3.1.2. External datasets

To assess the generalizability of our approach, we employed two external datasets: the Brazilian Multilabel Ophthalmological Dataset (BRSET) [25,26] and the mBRSET [27,28].

The Brazilian Multilabel Ophthalmological Dataset (BRSET) [25,26] is a multi-labeled ophthalmological dataset designed to reduce the under-representation of diverse populations within ophthalmological data and to facilitate the development and validation of machine learning models in ophthalmology. It comprises 16,266 images from 8524 diabetic and non-diabetic Brazilian patients, acquired in three Brazilian ophthalmological centers in São Paulo. These images were captured using Nikon NF505 and Canon CR-2 cameras, focusing on macula-centered views. In addition to fundus images, the dataset includes the following tabular features: patient's age and sex, insuline treatment, and diabetes duration. All Images were labeled by a single ophthalmologist.

mBRSET [27,28] is an extension of BRSET. It comprises 5164 images from 1291 patients of diverse backgrounds. These images were captured using handheld retinal cameras in real-life, high-burden scenarios. They were graded by two ophthalmologists for quality control and DR severity. This dataset provides the following clinical and demographic tabular features: patient's age and sex, educational Level, diabetes duration, insulin treatment, systemic hypertension, alcohol consumption, smoking, obesity, vascular disease, acute myocardial infarction, nephropathy, neuropathy, and diabetic foot.

Both datasets include labels for DR severity according to the ICDR grading scale.

### 3.1.3. Training process

For training, the batch size was set to 32. A maximum learning rate of 0.0001 was used in conjunction with a cosine annealing scheduler and a warmup period of 5 epochs. Finally, AdamW [29] was used as the optimizer. In all networks, we use the pretrained weights (on ImageNet) for the image feature extraction module. The training process fine-tunes the network in an end-to-end fashion using our data.

Patients in our dataset were divided into a fixed test set (207 patients, 669 images), and multiple training and validation sets through 5-fold cross-validation to ensure thorough evaluation and unbiased model performance estimation. This way, one patient's data will never be divided between two sets guarantying that there is no data leakage between the training and the test set. The best model across the 5-fold cross-validation is used on the test set to report the metrics in each experiment.

The training and evaluation of all models were done using a system with an Nvidia Tesla V100 GPU, an AMD EPYC 7402 processor, and 32 GB of memory.

## 3.2. Tabular classification on the EVIRED dataset

We perform a model sweep using the algorithms presented in Section 2.1. For evaluation, we use accuracy, quadratic weighted kappa,

<table><caption>**Table 1**

Description of tabular features gathered and measured for each patient.</caption><thead><tr><td>Type</td><td>No.</td><td>Examples</td></tr></thead><tbody><tr><td>Demographic data</td><td>3</td><td>Age; Sex; Ethnicity.</td></tr><tr><td>Medical history</td><td>30</td><td>Diabetes type and duration; High blood pressure; etc.</td></tr><tr><td>Ophthalmological history</td><td>13</td><td>Cataract; IOP lowering drops; Corticoid injections; etc.</td></tr><tr><td>General medical exam</td><td>12</td><td>Weight; Height; BMI; SAP; DAP; HbA1c; etc.</td></tr><tr><td>Ophthalmological exam</td><td>12</td><td>Visual Acuity; Refraction Sphere; Intraocular pressure; etc.</td></tr><tr><td>Health questionnaire</td><td>6</td><td>Ability to perform tasks: mobility, self-care, usual activities; Amount of pain, anxiety; General health status.</td></tr></tbody></table>

weighted F1 score, and the Area Under the ROC Curve (AUC) as metrics. AUC is computed separately for each class versus the rest. The weighted F1 score is the average F1 score for each class c versus the rest, weighted by the fraction of test samples belonging to c. Table 2 shows the measured metrics for each algorithm.

The confusion matrix for each of these algorithms shows different strengths and weaknesses. For example, Fig. 2 compares the results of TabPFN and Random Forest. While Random Forest is able to classify NoDR, Moderate and PDR cases better, TabPFN performs better on Mild and Severe cases. The low accuracy of PDR cases for both algorithms is related to the dataset containing very few PDR samples.

## 3.3. Multimodal classification on the EVIRED dataset

In this section, we integrate tabular and image classifiers to perform multimodal classification. We investigate whether or not multimodal methods outperform image-only classifiers.

First, we employ the ResNet50 architecture in conjunction with tabular classifiers mentioned in 3.2 to create a multimodal classifier. Table 3 shows how each combination performs.

Next, using the best-performing tabular model, we evaluated the performance of three versions of ResNet: ResNet34, ResNet50, and ResNet101. In addition, we tested the inclusion of Squeeze-and-excitation blocks [30] in ResNet50. For multi-modal classification, we used the same networks as the feature extractor block in Fig. 1. In each case, the best model is chosen from a 5-fold cross-validation training, and the evaluation is done on a fixed test set containing 669 samples. Table 4 compares the performance of these networks.

Overall, ResNet50 performs best for both imaging and multimodal data. Across all network architectures, we perceive improvements by using multimodal data. Fig. 3 illustrates the confusion matrices of DR classification using imaging and multimodal data. ResNet50 backbone was used for both cases. As we can see, there are clear improvements for all classes. In addition, with multimodal data, falsely classified samples (non-diagonal values) are closer to the right class (closer to the diameter). As discussed before, the classification results are poor for the PDR class due to the shortage of PDR training samples.

## 3.4. Statistical comparison

Tables 3 and 4 indicate that a combination of an MLP tabular classifier and a ResNet50 image classifier (MLPResNet50) produce the best overall result for DR grading in our test set. In this section, we present a statistical comparison between MLPResNet50 and a few other combinations presented in Section 3.3. We investigate whether or not the reported performance difference is statistically significant.

To that end, a Bhapkar test [31] is performed between each pair of methods to examine the proportion of samples labeled the same way by both methods for each class. Bhapkar's test is a marginal
<table><caption>**Table 2**

Classification metrics for multiple algorithms performing DR grading using tabular data.</caption><thead><tr><th>Algorithm</th><th>Accuracy</th><th>Kappa</th><th>Weighted F1</th><th>ROC AUC one-vs-rest</th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>MLP</td><td>0.432</td><td>0.530</td><td>0.425</td><td>0.721</td><td>0.602</td><td>0.641</td><td>0.674</td><td>0.638</td><td>0.921</td></tr><tr><td>Residual</td><td>0.459</td><td>0.569</td><td>0.470</td><td>0.719</td><td>0.657</td><td>0.656</td><td>0.627</td><td>0.635</td><td>0.918</td></tr><tr><td>FTTransformer</td><td>0.450</td><td>0.553</td><td>0.473</td><td>0.719</td><td>0.645</td><td>0.677</td><td>0.680</td><td>0.570</td><td>0.939</td></tr><tr><td>GANDALF</td><td>0.414</td><td>0.581</td><td>0.433</td><td>0.739</td><td>0.650</td><td>0.668</td><td>0.651</td><td>0.628</td><td>0.938</td></tr><tr><td>TabPFN</td><td>0.414</td><td>0.540</td><td>0.443</td><td>0.736</td><td>0.665</td><td>0.655</td><td>0.652</td><td>0.573</td><td>0.915</td></tr><tr><td>Random Forest</td><td>0.454</td><td>0.625</td><td>0.463</td><td>0.739</td><td>0.739</td><td>0.665</td><td>0.613</td><td>0.619</td><td>0.938</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** The confusion matrix for the output of (a) TabPFN and (b) Random Forest, when compared to the ground truth.</figcaption></figure>

<table><caption>**Table 3**

Classification metrics for multiple algorithms performing DR grading using both tabular and imaging data.</caption><thead><tr><th>Algorithm</th><th>Image</th><th>Accuracy</th><th>Kappa</th><th>Weighted F1</th><th>ROC AUC one-vs-rest</th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>MLP</td><td>ResNet50</td><td>0.634</td><td>0.807</td><td>0.626</td><td>0.892</td><td>0.786</td><td>0.803</td><td>0.869</td><td>0.578</td><td>0.986</td></tr><tr><td>Residual</td><td>ResNet50</td><td>0.599</td><td>0.799</td><td>0.605</td><td>0.877</td><td>0.732</td><td>0.766</td><td>0.830</td><td>0.693</td><td>0.957</td></tr><tr><td>FTTransformer</td><td>ResNet50</td><td>0.629</td><td>0.766</td><td>0.613</td><td>0.881</td><td>0.793</td><td>0.796</td><td>0.853</td><td>0.538</td><td>0.979</td></tr><tr><td>GANDALF</td><td>ResNet50</td><td>0.574</td><td>0.784</td><td>0.587</td><td>0.869</td><td>0.755</td><td>0.789</td><td>0.876</td><td>0.524</td><td>0.968</td></tr><tr><td>TabPFN</td><td>ResNet50</td><td>0.593</td><td>0.788</td><td>0.592</td><td>0.897</td><td>0.755</td><td>0.759</td><td>0.845</td><td>0.614</td><td>0.970</td></tr><tr><td>Random Forest</td><td>ResNet50</td><td>0.589</td><td>0.774</td><td>0.586</td><td>0.850</td><td>0.730</td><td>0.740</td><td>0.761</td><td>0.520</td><td>0.974</td></tr></tbody></table>

<table><caption>**Table 4**

Classification metrics for multiple network backbone structures performing DR grading using imaging and tabular data. When used in a multimodal setting, each network is paired with an MLP tabular block.</caption><thead><tr><th rowspan="2">Backbone</th><th colspan="3">Image</th><th colspan="3">Multimodal</th></tr><tr><th>Accuracy</th><th>Kappa</th><th>Weighted F1</th><th>Accuracy</th><th>Kappa</th><th>Weighted F1</th></tr></thead><tbody><tr><td>ResNet34</td><td>0.629</td><td>0.748</td><td>0.591</td><td>0.616</td><td>0.772</td><td>0.612</td></tr><tr><td>ResNet50</td><td>0.613</td><td>0.785</td><td>0.603</td><td>0.634</td><td>0.807</td><td>0.626</td></tr><tr><td>ResNet50 with SE</td><td>0.570</td><td>0.761</td><td>0.575</td><td>0.637</td><td>0.781</td><td>0.622</td></tr><tr><td>ResNet101</td><td>0.522</td><td>0.737</td><td>0.532</td><td>0.625</td><td>0.775</td><td>0.608</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** The confusion matrix for the output of (a) ResNet50 using imaging data and (b) MLPResNet50 using multimodal data, when compared to the ground truth.</figcaption></figure>
Table 5

Statistical comparison between MLPResNet50 fusion model and other combinations of classifiers. Reported numbers are p-values of a Bhapkar test done on a test set of 669 samples for each pair. Values smaller than 0.05 indicate a statistically significant difference between the two methods.

<table><thead><tr><th rowspan="2">Algorithm</th><th rowspan="2">Image</th><th colspan="2">Bhapkar test</th></tr><tr><th>Statistic</th><th>p-value</th></tr></thead><tbody><tr><td>Tabular</td><td>Image</td><td></td><td></td></tr><tr><td>-</td><td>ResNet50</td><td>44.09</td><td>2.22E-08</td></tr><tr><td>MLP</td><td>ResNet34</td><td>29.34</td><td>2.00E-05</td></tr><tr><td>MLP</td><td>ResNet50/SE</td><td>8.19</td><td>1.46E-01</td></tr><tr><td>MLP</td><td>ResNet101</td><td>33.60</td><td>3.00E-06</td></tr><tr><td>Residual</td><td>ResNet50</td><td>22.40</td><td>4.39E-04</td></tr><tr><td>FTTransformer</td><td>ResNet50</td><td>34.93</td><td>2.00E-06</td></tr><tr><td>GANDALF</td><td>ResNet50</td><td>81.25</td><td>4.60E-16</td></tr><tr><td>TabPFN</td><td>ResNet50</td><td>80.51</td><td>6.58E-16</td></tr></tbody></table>

homogeneity test that measures the agreement between two raters in categorical data. The associated *p*-value shows whether the two raters differ significantly from each other or whether the difference is the result of random chance.

Table 5 reports the p-values associated with the Bhapkar test for MLPResNet50 and a few selected algorithms. Values lower than 0.05 indicate a statistically significant difference between the classification results of the two methods.

As we can see, when compared to an image classifier, we get a *p*-value of $2.22 \times 10^{-8}$, which shows a statistically significant change for the fusion model. In addition, when compared to MLPResNet50/SE, we get a *p*-value of 0.146, which is greater than 0.05, indicating that the addition of SE blocks does not significantly change the performance of ResNet50 in the fusion model.

### 3.5. Late stage vs. intermediate fusion

As reported in Section 3.3, the fusion model MLPResNet50 performs the best for DR classification. This model performs intermediate fusion by concatenating the tabular and image feature vectors. Normally, we set the intermediate feature length $M = 64$. In this section, we report our experiments to find the best value for $M$, as well as a comparison with late stage fusion. For late fusion, we get an output vector of 6 values for our 6 classes from each modality and combine them for the final prediction.

Table 6 reports the classification metrics calculated on the test set for an MLPResNet50 fusion model in different configurations. Late fusion shows a much lower accuracy and F1-score value but similar kappa hinting at many wrongly classified samples but with small differences to the ground truth. Overall, considering multiple metrics, intermediate fusion convincingly outperforms late fusion. Table 6 also shows that $M = 64$ performs the best across multiple metrics. This may be due to the number of tabular features in our dataset 76. It is natural to assume that the optimal value for $M$ may depend both on the length of the image feature vector and the number of tabular features. As a result of this experiment, we have used $M = 64$ as the default value for the MLPResNet50 fusion model.

### 3.6. Inference times

In the field of medical imaging, it is always critical to develop methods with low inference times and to ensure that they run on basic hardware. In this section, we report the inference times for the MLPResNet50 fusion model and compare it with a unimodal image classifier. The times are calculated both on GPU (Nvidia Tesla V100) and CPU (AMD EPYC 7402). Table 7 reports these values. As we can see, the addition of tabular data adds only a small overhead in inference time. In addition, even on CPU the inference times remain practical for any real-world scenarios.

### 3.7. Evaluation on external datasets

In previous sections, we demonstrated that tabular clinical data can enhance the performance of DR classification when used in conjunction with imaging data. To further evaluate this claim, we train and test our model using two external datasets BRset and mBRset [26,27]. As described in Section 3.1.2, these two datasets contain multiple tabular features that can be used in conjunction with their imaging data to improve DR classification. Since BRset and mBRset contain a different number of tabular features than our dataset, we need to retrain our models on both datasets separately. It should be noted that the BRset contains many samples from non-diabetic patients. We have used only the portion of the dataset with diabetic patients.

Similar to the previous experiments, we trained an image classifier (ResNet50) and a fusion model (MLPResNet50 as described in Section 2.3) on each dataset separately. Tables 8 and 9 show the evaluation results for both datasets. For BRset, we can see an improvement of 0.065 for kappa and 0.017 for accuracy when the fusion model is used. Similarly for mBRset, kappa increases by 0.048 and accuracy by 0.012. We also statistically compared the two models using Bhapkar's test [31]: for BRset we get a p-value of $7.22 \times 10^{-3}$ and for mBRset a p-value of $2.81 \times 10^{-3}$. In both cases, the p-value is lower than 0.05 signaling a statistically significant difference between the image classifier and the fusion model.

In contrast to the EVIRED dataset, both BRset and mBRset have considerable number of PDR cases. As we can see in Tables 8 and 9 our fusion model shows very strong performance for the PDR class with 0.913 AUC for BRset and 0.958 for mBRset. This suggests that the low and variable performance of our model for the PDR class on the EVIRED dataset can be improved in the future when more PDR samples are added to the dataset.

To summarize, we were able to reproduce our results using two external datasets containing fundus RGB images and tabular clinical data. Once again the integration of tabular and imaging features outperformed an image-only approach.

### 3.8. Model interpretability and feature importance

In this section, we use multiple interpretability methods in conjunction with our model to investigate how the classification process works. Depending on the ground truth class, the model attributes different importance to different parts of the input data (image and tabular features). For imaging data, we can visualize this as a positive attribution heat-map, and for the tabular data, we show which features contribute the most to the predicted class.

The following methods have been used to explain the outputs of our model:

* **GradCAM:** Gradient-weighted Class Activation Mapping (Grad-CAM) [32] is an attribution technique for CNNs allowing for visualization of internal convolution layers. Usually, it is applied to the last layer of the network and provides high-level visual explanations. It uses the gradient of a target class flowing into the last layer to highlight important areas of the image for predicting that class [32].

* **Integrated Gradients:** Integrated Gradients (IG) [33] is an attribution method that aims to explain a model's output based on its input features. The algorithm computes the integral of gradients of the model function along a path from a baseline input to the actual input (in our case, from a black image to an RGB retina image). This integral is approximated by using a sum in *m* steps. In this work, we set *m* = 100.

* **SHAP:** SHapley values [34] is a popular concept from cooperative game theory that has been used recently to explain ML models [35]. We use the SHAP (SHapley Additive exPlanations) python package [36] to measure and visualize the importance of each tabular feature for the overall classification output.
<table><caption>**Table 6**
Comparison between late stage and intermediate fusion for the MLPResNet50 model. For intermediate fusion, results are reported for several different feature vector lengths *M*.</caption><thead><tr><th>Algorithm</th><th></th><th>Accuracy</th><th>Kappa</th><th>Weighted F1</th><th>ROC AUC one-vs-rest</th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td rowspan="2">Fusion</td><td>M</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>16</td><td>0.593</td><td>0.778</td><td>0.586</td><td>0.867</td><td>0.748</td><td>0.749</td><td>0.861</td><td>0.664</td><td>0.971</td></tr><tr><td rowspan="4">Intermediate</td><td>32</td><td>0.623</td><td>0.778</td><td>0.595</td><td>0.868</td><td>0.769</td><td>0.758</td><td>0.822</td><td>0.590</td><td>0.965</td></tr><tr><td>64</td><td>0.634</td><td>0.807</td><td>0.626</td><td>0.892</td><td>0.786</td><td>0.803</td><td>0.869</td><td>0.578</td><td>0.986</td></tr><tr><td>128</td><td>0.592</td><td>0.777</td><td>0.572</td><td>0.854</td><td>0.754</td><td>0.746</td><td>0.819</td><td>0.687</td><td>0.970</td></tr><tr><td>256</td><td>0.638</td><td>0.796</td><td>0.613</td><td>0.889</td><td>0.776</td><td>0.774</td><td>0.851</td><td>0.463</td><td>0.964</td></tr><tr><td colspan="2">Late</td><td>0.405</td><td>0.780</td><td>0.341</td><td>0.841</td><td>0.703</td><td>0.737</td><td>0.867</td><td>0.818</td><td>0.968</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Impact (negative or positive) on model output based on feature values for the top 10 features. Each dot represents a sample. (a) class 0: No DR; (b) class 2: Moderate DR; (c) class 5: PRP.</figcaption></figure>

<table><caption>**Table 7**
Average inference times per image for ResNet50 image classifier and MLPRes-Net50 fusion model on two devices.</caption><thead><tr><th>Model</th><th>GPU</th><th>CPU</th></tr></thead><tbody><tr><td>Image classifier</td><td>0.116 s</td><td>0.316 s</td></tr><tr><td>Image+Tabular fusion</td><td>0.118 s</td><td>0.316 s</td></tr></tbody></table>

**Fig. 6** shows the outputs of GradCAM and IG for a few different input images based on a ResNet50 classifier using the Captum Python library [37]. These generated interpretations have been reviewed by clinical experts in our team (AC, SB). They have found that the highlighted areas have a significant overlap with the regions of the image they consider when grading DR. The outputs of IG mostly highlight microaneurysms, exudates (in case of moderate or severe DR), and laser burns (for class PRP) present in the image. Similarly, GradCAM shows which high-level parts of the image are positively attributed to the classification result. As we can see, both methods generally agree in their assessments.

In case of incorrect classifications, the model may mistakenly focus on artifacts or normal features of the retina (optic disk, fovea or healthy vessels). An example of this is given in the last row of Fig. 6. In this example, the ground truth is *No DR* but the model classified it as a *moderate* case of DR. As we can see, the attributions mostly point to the optic disk and healthy vessels.

To analyze the contribution of each tabular feature, we have used Shapley values [34,38] in conjunction with the Random Forest tabular classifier.

The Shapley value $S_{ijk}$ is generated for the combination of sample $x_i$, feature $f_j$ ($j \in \{0, ..., 76\}$) and class $c_k$ ($k \in \{0, ..., 5\}$). Positive (Negative) values signify that the feature $f_j$ contributes positively (negatively) towards class $c_k$ being predicted for sample $x_i$. Higher $|S_{ijk}|$ indicate higher significance of feature $f_j$.

**Fig. 4** shows the contribution of the top ten most significant features to the model output for a specific class. In each chart, samples are shown by a dot: the blue dots have low feature values, while the red ones have high feature values. In the case of binary features, red means 1, and blue means 0. Dots on the right of the center show a positive impact on the output, and dots on the left show a negative impact. Overall, we can see the relationship between the value of a feature and the output for that class. For example, for the class *No DR* in Fig. 4(a), we can see that high visual acuity has a strong positive impact on the model output in favor of *No DR*. On the other hand, high duration of diabetes has a negative impact for class *No DR* while having a positive impact for class *Moderate* (Fig. 4(b)) and *PRP* (Fig. 4(c)).

Overall, across all six classes, we can rank the features from the most important to the least. Table 10 mentions the top ten most significant features. The results reported in Fig. 4, and Table 10 are in line with the previously reported risk factors of DR. According to Yau et al. [39] and more recently Bek [40], diabetes duration, $HbA1c$ and blood pressure are the three major risk factors of DR. Moreover, they report that DR is more prevalent among patients with diabetes type 1 than type 2 which is in agreement with our findings. Nonetheless, it is important to note that each feature's contribution can result from a real-world causal relationship with DR or simply due to a sampling bias present in the dataset. For example, Age in Fig. 4(a) appears to have an inverse relationship to the *No DR* class. We know that this is the result of not having many young patients with no DR in the dataset. This is why careful study of model explanations can help correct the participation requirements of future studies.

Similarly, we can examine a specific sample and its corresponding Shapley values to figure out which features contributed the most to the decision of the model. As an example, Fig. 5 illustrates how the top three most important features (for this sample) contribute to the model output (*No DR*). The blue (red) color indicates a positive (negative) attribution, and the thickness of each line represents how big the
<table><caption>**Table 8**
Classification metrics for DR grading on the BRset [26] dataset. A comparison between an image classifier and a fusion model using tabular and imaging data.</caption><thead><tr><th>Algorithm</th><th></th><th>Accuracy</th><th>Kappa</th><th>Weighted F1</th><th>ROC AUC one-vs-rest</th></tr></thead><tbody><tr><td>Tabular</td><td>Image</td><td></td><td></td><td></td><td></td></tr><tr><td>-</td><td>ResNet50</td><td>0.642</td><td>0.699</td><td>0.622</td><td>0.936</td><td>0.833</td><td>0.782</td><td>0.741</td><td>0.901</td></tr><tr><td>MLP</td><td>ResNet50</td><td>0.659</td><td>0.764</td><td>0.659</td><td>0.912</td><td>0.836</td><td>0.777</td><td>0.910</td><td>0.913</td></tr></tbody></table>

<table><caption>**Table 9**
Classification metrics for DR grading on the mBRset [27] dataset. A comparison between an image classifier and a fusion model using tabular and imaging data.</caption><thead><tr><th>Algorithm</th><th></th><th>Accuracy</th><th>Kappa</th><th>Weighted F1</th><th>ROC AUC one-vs-rest</th></tr></thead><tbody><tr><td>Tabular</td><td>Image</td><td></td><td></td><td></td><td></td></tr><tr><td>-</td><td>ResNet50</td><td>0.620</td><td>0.723</td><td>0.590</td><td>0.848</td><td>0.646</td><td>0.866</td><td>0.944</td><td>0.962</td></tr><tr><td>MLP</td><td>ResNet50</td><td>0.632</td><td>0.771</td><td>0.601</td><td>0.850</td><td>0.683</td><td>0.843</td><td>0.969</td><td>0.958</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Impact (negative or positive) on model output based on top three feature values for one specific sample. The Blue (red) color indicates a positive (negative) impact.</figcaption></figure>

<table><caption>**Table 10**
Overall top ten most important tabular features according to the Random Forest classifier and Shapley values.</caption><thead><tr><th>Rank</th><th>Feature name</th><th>Stats in the dataset</th></tr></thead><tbody><tr><td>1</td><td>Panretinal photocoagulation</td><td>14.2% yes; 81.5% no</td></tr><tr><td>2</td><td>Duration of diabetes</td><td>Min: 0; Max: 64; Median: 20</td></tr><tr><td>3</td><td>Age</td><td>Min: 21; Max: 97; Median: 60</td></tr><tr><td>4</td><td>Diabetes type</td><td>25.3% type 1; 67.2% type 2</td></tr><tr><td>5</td><td>Hemoglobin A<sub>1C</sub></td><td>Min: 4.6%; Max: 18.4%; Median: 7.5%</td></tr><tr><td>6</td><td>Visual acuity</td><td>Min: 0.01; Max: 2.0; Median: 1.0</td></tr><tr><td>7</td><td>Triglycerides</td><td>Min: 0.11; Max: 9.0; Median: 0.86</td></tr><tr><td>8</td><td>Total cholesterol</td><td>Min: 0.38; Max: 5.38; Median: 1.25</td></tr><tr><td>9</td><td>Serum creatinine</td><td>Min: 10.0; Max: 1096.0; Median: 64.0</td></tr><tr><td>10</td><td>Systolic arterial pressure</td><td>Min: 79.0; Max: 230.0; Median: 125.0</td></tr></tbody></table>

contribution is. For this sample, the classification is made correctly with the help of mainly two features: normal visual acuity and low duration of diabetes. In general, this chart can help explain not only correct classifications but also why the model fails on certain samples.

### 3.9. Tabular feature selection

In Section 3.8, we reported the most important tabular features for DR classification. In this section, we investigate the variation in performance of the fusion model if we limit the number of features to only the top-f most important ones. The feature selection is performed using the SHAP values as described in Section 3.8.

We train the MLPResNet50 fusion model using top-5, top-10 and top-20 features and compare them to the performance using the full set. Table 11 reports the classification metrics for all four cases. As we can see, classification using limited sets of features does not match the performance of the full set. However, classification with the top-10 features come close with similar accuracy but lower kappa, F1 score, and AUC values. This shows that while achieving good performance is possible with a limited set of clinical features, it is still vital to assemble an extensive set of features to reach optimum performance.

### 4. Conclusions and future works

In this work, we presented an extensive overview of different methods to integrate imaging and tabular data. Specifically, we showed that the use of clinical data in conjunction with UWF-CFP images significantly improves the accuracy of DR severity assessment. This was done by performing an ablation study on DR classification using each modality alone compared to their fusion. We examined various tabular classification models, and image feature extraction backbones to find an optimal combination for multimodal classification. In addition, several data pre-processing steps were introduced to optimize the tabular data for deep learning models.

Furthermore, various interpretation methods were used to showcase the inner workings of the network on UWF-CFP images. This can help build trust with the practitioners by illustrating the reasons behind the network's predictions. With regard to the tabular data, this study highlights which features are useful to help assess DR severity and
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** (a),(d),(g),(j) input image; (b),(e),(h),(k) GradCAM output for the last convolutional layer blended with the input image; (c),(f),(i),(l) Integrated Gradients output blended with the input image. All cases work with a ResNet50 classifier. The first three rows show cases where the model predicted the correct class while the last row shows a failed case.</figcaption></figure>

<caption>**Table 11**</caption>

Classification metrics performing DR grading for the MLPResNet50 fusion model using the top-$f$ most important tabular features. $f = 76$ represents the full set of features.

<table><thead><tr><th>$f$</th><th>Accuracy</th><th>Kappa</th><th>Weighted F1</th><th colspan="6">ROC AUC one-vs-rest</th></tr></thead><tbody><tr><td>5</td><td>0.614</td><td>0.785</td><td>0.593</td><td>0.838</td><td>0.712</td><td>0.771</td><td>0.796</td><td>0.632</td><td>0.976</td></tr><tr><td>10</td><td>0.637</td><td>0.785</td><td>0.604</td><td>0.841</td><td>0.702</td><td>0.774</td><td>0.774</td><td>0.568</td><td>0.980</td></tr><tr><td>20</td><td>0.599</td><td>0.788</td><td>0.588</td><td>0.877</td><td>0.685</td><td>0.730</td><td>0.837</td><td>0.648</td><td>0.962</td></tr><tr><td>76</td><td>0.634</td><td>0.807</td><td>0.626</td><td>0.892</td><td>0.786</td><td>0.803</td><td>0.869</td><td>0.578</td><td>0.986</td></tr></tbody></table>
<table><caption>**Table A.12**
List of 76 tabular features gathered and measured for each patient.</caption><thead><tr><th>Type</th><th>No.</th><th>Features</th></tr></thead><tbody><tr><td>Demographic data</td><td>3</td><td>Age; Sex; Ethnicity.</td></tr><tr><td>Medical history</td><td>30</td><td>Diabetes type, duration and treatment; High blood pressure; Lipid-lowering therapy; Smoking; Cardiovascular events; Plantar foot ulcer/amputation; Nephropathy; Neuropathy; Sleep Apnea; Pregnancy; Bariatric surgery; Historic HbA1c measurements; Dialysis or kidney transplant; Pancreatic transplant; Hypertension treatment.</td></tr><tr><td>Ophthalmological history</td><td>13</td><td>Cataract; IOP lowering drops; Steroid intravitreal treatment; Dexamethasone/Fluocinolone implant; Corticoid injections; Anti-VEGF intravitreal treatment; Vitrectomy; Complication of PDR/Macular surgery; Focal laser treatment; Non-central DME threatening the fovea; Presence of macular microaneurysms.</td></tr><tr><td>General medical exam</td><td>12</td><td>Weight; Height; BMI; No. of severe hypoglycemias; Continuous glucose monitoring; Time spent at the glycemic target; SAP; DAP; HbA1c; Serum creatinine; Total cholesterol; Triglycerides.</td></tr><tr><td>Ophthalmological exam</td><td>12</td><td>Visual acuity; Refraction sphere, Cylinder and axis; Intraocular pressure; Iris rubeosis; Pseudophakia; Significative lens opacity; Capsulotomy; Macular laser burns; Panretinal photocoagulation; Intravitreal haemorrhage.</td></tr><tr><td>Health questionnaire</td><td>6</td><td>Ability to perform tasks: mobility, self-care, usual activities; Amount of pain, anxiety; General health status.</td></tr></tbody></table>

which are not. This knowledge could potentially help streamline DR screening programs, without compromising patient's safe diagnosis.

Currently, there are no similar public datasets that include both UWF-CFP images and rich tabular data. Nevertheless, we have trained and evaluated our model using two similar external datasets that contain regular CFP images (as opposed to ultra-widefield images in our dataset) and a limited set of clinical tabular data. We established that in both datasets our fusion model outperforms an equivalent image classifier.

The use of clinical tabular data alongside image deep learning models in ophthalmology is still a new concept. Therefore, it is imperative that future studies further explore different fusion frameworks similar to other medical domains [41]. Furthermore, this study focused only on one imaging device (CLARUS 500), so compatibility with other devices (such as Optos, Eidon, etc.) can be a useful topic to explore [42]. In the wider context of multimodal DR analysis, other modalities such as hyperspectral imaging have also proved a valuable source of information [43]. Exploring the use of these modalities in conjunction with patient tabular data can be another future avenue of study.

As mentioned in Section 3.8, some relationships found between features and the output of the model are solely the result of population biases present in the current clinical study. Therefore, future studies should focus on population independence to generalize this work [44].

The EVIRED project is an ongoing project which monitors the progression of DR by performing followups on each patient every year. For that reason, the dataset still does not contain many cases of Proliferate DR (class 4 PDR in our classifications). This has led to low and varying classification performance for the PDR class across all of our models. We addressed this limitation by training and evaluating our model on two external datasets: the BRset and mBRset datasets, which have a considerable number of PDR samples. Our fusion model performs strongly in this class for both datasets showing the possibility of higher performance for the PDR class on the EVIRED dataset if more PDR cases become available in the future.

Lastly, this study can be extended by taking into account previous medical examinations instead of just the current one. In this way, the classification problem proposed here turns into a progression prediction problem based on several previous data points [45,46].

### CRediT authorship contribution statement

**Alireza Rezaei:** Writing – original draft, Visualization, Methodology, Investigation, Data curation. **Sarah Matta:** Writing – review & editing, Data curation. **Rachid Zeghlache:** Writing – review & editing. **Pierre-Henri Conze:** Writing – review & editing, Supervision.

**Capucine Lepicard:** Project administration. **Pierre Deman:** Funding acquisition, Conceptualization. **Laurent Borderie:** Funding acquisition, Conceptualization. **Deborah Cosette:** Funding acquisition, Conceptualization. **Sophie Bonnin:** Resources, Formal analysis. **Aude Couturier:** Resources, Project administration, Funding acquisition, Formal analysis, Conceptualization. **Béatrice Cochener:** Resources, Formal analysis. **Mathieu Lamard:** Writing – review & editing, Supervision. **Mostafa El Habib Daho:** Writing – review & editing, Supervision, Methodology. **Gwenolé Quellec:** Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Aude Couturier reports financial support was provided by Agence Nationale de la Recherche. Pierre Deman and Laurent Borderie report a relationship with Evolucare Technologies that includes: employment. Deborah Cosette reports a relationship with Carl Zeiss Meditec Inc that includes: employment. Beatrice Cochener reports a relationship with Carl Zeiss Meditec Inc that includes: consulting or advisory. Gwenole Quellec reports a relationship with Evolucare Technologies that includes: consulting or advisory. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The work takes place in the framework of EVIRED, an ANR RHU project. This work benefits from State aid managed by the French National Research Agency under "Investissement d'Avenir" program bearing the reference ANR-18-RHUS-0008.

### Appendix. Clinical contextual data

See Table A.12.

### Data availability

The data that has been used is confidential.
References

[1] Z.L. Teo, Y.-C. Tham, M. Yu, M.L. Chee, T.H. Rim, N. Cheung, M.M. Bikbov, Y.X. Wang, Y. Tang, Y. Lu, et al., Global prevalence of diabetic retinopathy and projection of burden through 2045: systematic review and meta-analysis, Ophthalmol 128 (11) (2021) 1580-1591.

[2] K.-K. Maninis, J. Pont-Tuset, P. Arbeláez, L. Van Gool, Deep retinal image understanding, in: Medical Image Computing and Computer-Assisted Intervention-MICCAI 2016: 19th International Conference, Athens, Greece, October 17-21, 2016, Proceedings, Part II 19, Springer, 2016, pp. 140-148.

[3] N.A. El-Hag, A. Sedik, W. El-Shafai, H.M. El-Hoseny, A.A. Khalaf, A.S. El-Fishawy, W. Al-Nuaimy, F.E. Abd El-Samie, G.M. El-Banby, Classification of retinal images based on convolutional neural network, Microsc. Res. Tech. 84 (3) (2021) 394-414.

[4] D.S.W. Ting, C.Y.-L. Cheung, G. Lim, G.S.W. Tan, N.D. Quang, A. Gan, H. Hamzah, R. Garcia-Franco, I.Y. San Yeo, S.Y. Lee, et al., 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) (2017) 2211-2223.

[5] E. Ho, E. Wang, S. Youn, A. Sivajohan, K. Lane, J. Chun, C.M. Hutnik, Deep ensemble learning for retinal image classification, Transl. Vis. Sci. Technol. 11 (10) (2022) 39.

[6] P.S. Silva, A.J. Dela Cruz, M.G. Ledesma, J. van Hemert, A. Radwan, J.D. Cavalleroano, L.M. Aiello, J.K. Sun, L.P. Aiello, Diabetic retinopathy severity and peripheral lesions are associated with nonperfusion on ultrawide field angiography, Ophthalmol 122 (12) (2015) 2465-2472, http://dx.doi.org/10.1016/j.ophtha.2015.07.034.

[7] S. Cai, T.A. Liu, The role of ultra-widefield fundus imaging and fluorescein angiography in diagnosis and treatment of diabetic retinopathy, Curr. Diabetes Rep. 21 (9) (2021) 30.

[8] M. El Habib Daho, Y. Li, R. Zeghlache, Y.C. Atse, H.L. Boité, S. Bonnin, D. Cosette, P. Deman, L. Borderie, C. Lepicard, R. Tadayoni, B. Cochener, P. Conze, M. Lamard, G. Quellec, Improved automatic diabetic retinopathy severity classification using deep multimodal fusion of UWF-CFP and OCTA images, in: OMIA Workshop, in: LNCS, vol. 14096, Springer, 2023, pp. 11-20, http://dx.doi.org/10.1007/978-3-031-44013-7_2.

[9] J.H. Friedman, Greedy function approximation: a gradient boosting machine, Ann. Stat. (2001) 1189-1232.

[10] Y. Gorishniy, I. Rubachev, V. Khrulkov, A. Babenko, Revisiting deep learning models for tabular data, Adv. Neural Inf. Process. Syst. 34 (2021) 18932-18943.

[11] M. Joseph, H. Raj, GANDALF: gated adaptive network for deep automated learning of features, 2022, arXiv preprint arXiv:2207.08548.

[12] K. Cho, B. Van Merriënboer, D. Bahdanau, Y. Bengio, On the properties of neural machine translation: Encoder-decoder approaches, 2014, arXiv preprint arXiv:1409.1259.

[13] N. Hollmann, S. Müller, K. Eggensperger, F. Hutter, TabPFN: A transformer that solves small tabular classification problems in a second, 2022, arXiv preprint arXiv:2207.01848.

[14] N. Hollmann, S. Müller, L. Purucker, A. Krishnakumar, M. Körfer, S.B. Hoo, R.T. Schirmeister, F. Hutter, Accurate predictions on small data with a tabular foundation model, Nature 637 (8045) (2025) 319-326.

[15] Y. Li, M. El Habib Daho, P.-H. Conze, R. Zeghlache, H. Le Boité, R. Tadayoni, B. Cochener, M. Lamard, G. Quellec, A review of deep learning-based information fusion techniques for multimodal medical image classification, Comput. Biol. Med. 177 (2024) 108635, http://dx.doi.org/10.1016/j.combiomed.2024.108635, URL https://www.sciencedirect.com/science/article/pii/S0010482524007200.

[16] M. Artzi, E. Redmard, O. Tzemach, J. Zeltser, O. Gropper, J. Roth, B. Shofty, D.A. Kozyrev, S. Constantini, L. Ben-Sira, Classification of pediatric posterior fossa tumors using convolutional neural network and tabular data, IEEE Access 9 (2021) 91966-91973.

[17] T.N. Wolf, S. Pölsterl, C. Wachinger, A.D.N. Initiative, et al., Daft: A universal module to interweave tabular data and 3D images in CNNs, NeuroImage 260 (2022) 119505.

[18] M.-Y. Hsu, J.-Y. Chiou, J.-T. Liu, C.-M. Lee, Y.-W. Lee, C.-C. Chou, S.-C. Lo, E. Kornelius, Y.-S. Yang, S.-Y. Chang, et al., Deep learning for automated diabetic retinopathy screening fused with heterogeneous data from EHRs can lead to earlier referral decisions, Transl. Vis. Sci. Technol. 10 (9) (2021) 18.

[19] L. Breiman, Random forests, Mach Learn. 45 (2001) 5-32.

[20] N.V. Chawla, K.W. Bowyer, L.O. Hall, W.P. Kegelmeyer, SMOTE: synthetic minority over-sampling technique, J. Artificial Intelligence Res. 16 (2002) 321-357.

[21] T. Li, W. Bo, C. Hu, H. Kang, H. Liu, K. Wang, H. Fu, Applications of deep learning in fundus images: A review, 2021, arXiv:2101.09864. URL https://arxiv.org/abs/2101.09864.

[22] G. Quellec, H. Al Hajj, M. Lamard, P.-H. Conze, P. Massin, B. Cochener, ExplAln: Explanatory artificial intelligence for diabetic retinopathy diagnosis, Med. Image Anal. 72 (2021) 102118, http://dx.doi.org/10.1016/j.media.2021.102118.

[23] F. Shamshad, S. Khan, S.W. Zamir, M.H. Khan, M. Hayat, F.S. auz Khan, H. Fu, Transformers in medical imaging: A survey, Med. Image Anal. 88 (2023) 102802, http://dx.doi.org/10.1016/j.media.2023.102802.

[24] M. El Habib Daho, Y. Li, R. Zeghlache, H.L. Boité, P. Deman, L. Borderie, H. Ren, N. Mannivanan, C. Lepicard, B. Cochener, A. Couturier, R. Tadayoni, P.-H. Conze, M. Lamard, G. Quellec, DISCOVER: 2-D multiview summarization of optical coherence tomography angiography for automatic diabetic retinopathy diagnosis, Artif. Intell. Med. 149 (2024) 102803, http://dx.doi.org/10.1016/j.artmed.2024.102803.

[25] L.F. Nakayama, M. Goncalves, L. Zago Ribeiro, H. Santos, D. Ferraz, F. Malerbi, L.A. Celi, C. Regatieri, A Brazilian multilabel ophthalmological dataset (BRSET), PhysioNet (2023) 2, http://dx.doi.org/10.13026/xcxw-8198.

[26] L.F. Nakayama, D. Restrepo, J. Matos, L.Z. Ribeiro, F.K. Malerbi, L.A. Celi, C.S. Regatieri, BRSET: a Brazilian multilabel ophthalmological dataset of retina fundus photos, PLOS Digit. Hea! 3 (7) (2024) e0000454.

[27] L.F. Nakayama, L. Zago Ribeiro, D. Restrepo, et al., mBRSET, a mobile Brazilian retinal dataset, 2024.

[28] C. Wu, D. Restrepo, L.F. Nakayama, L. Zago Ribeiro, Z. Shuai, N.S. Barboza, M.L.V. Sousa, R.D. Fitterman, A.D.A. Pereira, C.V.S. Regatieri, et al., A portable retina fundus photos dataset for clinical, demographic, and diabetic retinopathy prediction, Sci. Data 12 (1) (2025) 323.

[29] I. Loshchilov, F. Hutter, Decoupled weight decay regularization, 2017, arXiv preprint arXiv:1711.05101.

[30] J. Hu, L. Shen, G. Sun, Squeeze-and-excitation networks, in: IEEE/CVF Conf Comput Vis Pattern Recog, 2018, pp. 7132-7141, http://dx.doi.org/10.1109/CVPR.2018.00745.

[31] V.P. Bhapkar, A note on the equivalence of two test criteria for hypotheses in categorical data, J. Amer. Statist. Assoc. 61 (313) (1966) 228-235.

[32] R.R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, D. Batra, Gradcam: Visual explanations from deep networks via gradient-based localization, in: Proceedings of the IEEE International Conference on Computer Vision, 2017, pp. 618-626.

[33] M. Sundararajan, A. Taly, Q. Yan, Axiomatic attribution for deep networks, in: International Conference on Machine Learning, PMLR, 2017, pp. 3319-3328.

[34] L.S. Shapley, A value for n-person games, Contrib. To Theory Games 2 (1953).

[35] M. Sundararajan, A. Najmi, The many Shapley values for model explanation, in: International Conference on Machine Learning, PMLR, 2020, pp. 9269-9278.

[36] S.M. Lundberg, G. Erion, H. Chen, A. DeGrave, J.M. Prutkin, B. Nair, R. Katz, J. Himmelfarb, N. Bansal, S.-I. Lee, From local explanations to global understanding with explainable AI for trees, Nat. Mach. Intell. 2 (1) (2020) 2522-25839.

[37] N. Kokhlikyan, V. Miglani, M. Martin, E. Wang, B. Alsallakh, J. Reynolds, A. Melnikov, N. Kliushkina, C. Araya, S. Yan, et al., Captum: A unified and generic model interpretability library for pytorch, 2020, arXiv preprint arXiv: 2009.07896.

[38] L. Merrick, A. Taly, The explanation game: Explaining machine learning models using shapley values, in: Machine Learning and Knowledge Extraction: International Cross-Domain Conference, CD-MAKE 2020, Springer, 2020, pp. 17-38.

[39] J.W. Yau, S.L. Rogers, R. Kawasaki, E.L. Lamoureux, J.W. Kowalski, T. Bek, S.-J. Chen, J.M. Dekker, A. Fletcher, J. Grauslund, et al., Global prevalence and major risk factors of diabetic retinopathy, Diabetes Care 35 (3) (2012) 556-564.

[40] T. Bek, Systemic risk factors contribute differently to the development of proliferative diabetic retinopathy and clinically significant macular oedema, Diabetologia 63 (2020) 2462-2470.

[41] D. Duenias, B. Nichyporuk, T. Arbel, T.R. Raviv, HyperFusion: A hypernetwork approach to multimodal integration of tabular and medical imaging data for predictive modeling, 2024, arXiv preprint arXiv:2403.13319.

[42] M. El Habib Daho, Y. Li, R. Zeghlache, A. Rezaei, H. Le Boité, A. Couturier, S. Bonnin, S. Magazzeni, A. Le Guilcher, F. Potevin, et al., Cross-device AI fusion: Enhancing diabetic retinopathy diagnosis with combined clarus and optos images, Investig. Ophthalmol. Vis. Sci. 65 (7) (2024) 5630.

[43] C.-Y. Wang, A. Mukundan, Y.-S. Liu, Y.-M. Tsao, F.-C. Lin, W.-S. Fan, H.-C. Wang, Optical identification of diabetic retinopathy using hyperspectral imaging, J. Pers. Med. 13 (6) (2023) 939.

[44] S. Matta, M. Lamard, P. Zhang, A.L. Guilcher, L. Borderie, B. Cochener, G. Quellec, A systematic review of generalization research in medical image classification, 2024, arXiv preprint arXiv:2403.12167.

[45] L. Dai, B. Sheng, T. Chen, Q. Wu, R. Liu, C. Cai, L. Wu, D. Yang, H. Hamzah, Y. Liu, et al., A deep learning system for predicting time to progression of diabetic retinopathy, Nature Med. 30 (2) (2024) 584-594.

[46] R. Zeghlache, P.-H. Conze, M. El Habib Daho, Y. Li, A. Rezaei, H.L. Boité, R. Tadayoni, P. Massin, B. Cochener, I. Brahim, et al., L-MAE: Longitudinal masked auto-encoder with time and severity-aware encoding for diabetic retinopathy progression prediction, Comput. Biol. Med. (2025) http://dx.doi.org/10.1016/j.combiomed.2024.109508.