Interpretable surrogate models to approximate the predictions of convolutional neural networks in glaucoma diagnosis

Two different datasets have been used in this work. Some experiments were carried out with our own dataset composed of 606 eye fundus images with which the Deep Learning models were trained and tested. From these images, 399 correspond to glaucoma cases and 207 to healthy subjects. The images were captured in infrarred (815 nm) in the Hospital Universitario de Canarias with a Spectralis SD-OCT. The fact of using infrared images is only a matter of convenience to facilitate their acquisition due to the logistics of the hospital. Images in the visible could also have been used for the experiments we have performed. This dataset will be called 'Dataset IR' throughout the rest of this paper.

Patients with a diagnosis of primary or secondary open-angle glaucoma with untreated intraocular pressure greater than 21 mmHg were included in the study. The diagnosis of glaucoma was not only based on the observation of the eye fundus images but on the presence of reproducible defects in the white-white perimetry and/or morphological criteria based on a spectral domain optical coherence tomograph SD-OCT Spectralis with the Glaucoma Premium and Posterior Pole module. Both eyes were included if they met the inclusion criteria. Subjects with concomitant ocular pathology other than glaucoma, lower visual acuity of 20/40, refractive error greater than five diopters of spherical equivalent or three diopters of astigmatism, level of false positives, negatives and fixation errors equal to or greater than 25% in the visual field were excluded from the study. Patients with hypoplastic or oblique optic nerves were also excluded.

As mentioned in section 1, a simplified version of the eye fundus images showing only the contours of the disc and the cup on a gray background was also considered, as shown in an example in figure 2. More specifically, the pixels in the background of the image were assigned the value 128 and the pixels in the edges the value 255. These contours have been traced by an ophthalmologist specialized in the diagnosis of glaucoma, with more than 15 years of experience. In what follows, this variant of the first dataset will be called 'Dataset C'.

Zoom In Zoom Out Reset image size

The other dataset used in this work for training and testing the neural networks is an extended version of the previous one that includes, apart from our samples, other samples from very popular and publicly available image sets such as RIM–ONE DL, REFUGE and DRISHTI–GS. Thus, with this extension we have 1271 images corresponding to healthy subjects and 721 images corresponding to subjects with glaucoma. This set of images has only been used in its simplified version with the only information of the disc and cup contours, and will be called 'Dataset XC' in the following sections.

RIM–ONE DL [39] results from combining the three previous versions of the popular RIM-ONE dataset [40], [41]. This new version consists of 313 retinographies from normal subjects and 172 retinographies from patients with glaucoma. All of these images have been assessed by two experts and include a manual segmentation of the disc and cup. The REFUGE challenge database [42] consists of 1200 retinal images. About 10% of the dataset (120 samples) corresponds to glaucomatous subjects, including Primary Open Angle Glaucoma and Normal Tension Glaucoma. Ground truth segmentation of disc and cup is also included. This set of images is provided in three parts: training, validation and test, each with 400 images. In this work only the training and test partitions have been used, the validation one was ignored. DRISHTI-GS consists of 50 training and 51 testing images of which 70 correspond to subjects with glaucoma and 31 to healthy subjects [43]. For each image, manual segmentations were collected for both OD and cup region from four different human experts with varying clinical experience.

Each of the datasets used in this study underwent a partitioning process into training and testing sets. In the case of our dataset, across all its variations, we performed a random split, allocating 80% of the data for training and 20% for testing. For the other datasets, we adhered to their respective predefined splits, as published in the literature, to establish training and testing subsets.

Table 1 provides a concise overview of the datasets composition, detailing the number of samples per dataset and class. Additionally, it specifies the count of samples allocated for training and testing within each dataset.

Glaucoma diagnosis is a difficult problem due to the lack of reliable biomarkers so the diagnosis is usually based on a combination of information including medical history and both structural and functional longitudinal tests. In spite of this, in the case of diagnosis based on fundus imaging, some patterns have been identified that aid in the discrimination between healthy and glaucomatous eyes. As shown in figure 2, the optic disc usually appears in the image as a vertically oval shaped area. The cup is another important region inside the disc which usually appears as a brighter intensity area, as a result of the loss of nerve fibers caused by the disease. The part that remains between the cup and the disc is known as the neuroretinal rim.

As noted in [44, 45] glaucoma is characterized by progressive narrowing of the neuroretinal rim, taking the form of diffuse narrowing, localized notching, or a combination of both. Figure 3 shows the progression of glaucomatous damage at the optic disc. In this respect, there is a certain consensus when it comes to identifying certain areas of the rim as those that are usually the main affected in the case of suffering from glaucoma, depending on the state of the disease. Therefore, the geometry of the disc, cup and rim can provide very valuable information for the diagnosis of glaucoma.

Zoom In Zoom Out Reset image size

Taking into account the above and taking as a reference what is exposed in [45], 14 features related to the shape and dimensions of the optic disc and cup have been calculated for this work, which are detailed in table 2. It should be noted that both the disc and the cup have been manually delimited by an expert with more than 15 years of experience in the field. Features have been extracted from these contours using standard image processing techniques.

Table 2. List of the 14 features used in this work related to the shape and dimensions of the optic disc, cup and neuroretinal rim. Their number, mathematical symbol and abbreviation are also indicated, as used throughout the rest of the article.

No.	Symbol (Abbr.)	Feature
1	$a_\mathrm{OD}$ (OD size)	Optic disc size
2	$s_\mathrm{OD}$ (OD shape)	Optic disc shape
3	$a_\mathrm{NR}$ (NR size)	Neuroretinal rim size
4	I-ISNT	Neuroretinal rim shape (I—ISNT rule)
5	S-ISNT	Neuroretinal rim shape (S—ISNT rule)
6	N-ISNT	Neuroretinal rim shape (N—ISNT rule)
7	T-ISNT	Neuroretinal rim shape (T—ISNT rule)
8	$w^\mathrm{min.rim}_\mathrm{NR}$ (NR width)	Neuroretinal rim shape (narrowest rim width)
9	$w_\mathrm{RDR}$ (RD ratio)	Rim-to-disc ratio
10	$a_\mathrm{OC}$ (OC size)	Optic cup size
11	$s_\mathrm{OC}$ (OC shape)	Optic cup shape
12	$a_\mathrm{CDR}$ (CDR size)	Optic cup size in relation to the optic disc size
13	$v_\mathrm{CDR}$ (V CDR)	Vertical cup-to-disc ratio
14	$h_\mathrm{CDR}$ (H CDR)	Horizontal cup-to-disc ratio

The size of the disc ($a_\mathrm{OD}$), cup ($a_\mathrm{OC}$) and neuroretinal rim ($a_\mathrm{NR}$) was calculated as the area in number of pixels of the corresponding region in each image. The optic cup size in relation to the optic disc size was also calculated:

$$\begin{align} a_\mathrm{CDR}: = \frac{a_\mathrm{OC}}{a_\mathrm{OD}}. \end{align} \tag{ 1 }$$

Quantifying the shape of an optic disc requires the delimitation of a bounding box (figure 4(a)) where the horizontal and vertical diameters are defined (h_d , v_d ). The corresponding diameters for the optic cup are similarly defined (h_c , v_c ). Then a shape parameter for the optic disc and the optic cup ($s_\mathrm{OD}$, $s_\mathrm{OC}$) can be calculated as:

Zoom In Zoom Out Reset image size

$$\begin{align} s_\mathrm{OD} & : = \frac{h_d}{v_d} \nonumber s_\mathrm{OC} & : = \frac{h_c}{v_c}. \end{align} \tag{ 2 }$$

To capture the shape of the neuroretinal rim, different features were calculated. On the one hand, the thickness of the rim was determined at different points according to the ISNT rule [46] where I represents the Inferior sector of the optic disc, S the Superior sector, N the Nasal sector, and T the Temporal sector (figure 4(b)).

On the other hand, if $A^{\text{min.rim}}_\mathrm{OD}$ represents the axis of the optic disc with the narrowest rim, its length ($L^{\text{min.rim}}_\mathrm{OD}$) and the width of the neuroretinal rim in the same axis ($w^{\text{min.rim}}_\mathrm{NR}$) (figure 4(b)) are obtained, and used to calculate the rim-to-disc ratio ($w_\mathrm{RDR}$) [47]:

$$\begin{align} w_\mathrm{RDR}: = \frac{w^{\text{min.rim}}_\mathrm{NR}}{L^{\text{min.rim}}_\mathrm{OD}}. \end{align} \tag{ 3 }$$

Finally, two other measures were calculated, the ratio of both vertical and horizontal cup and disc diameters ($v_\mathrm{CDR}$, $h_\mathrm{CDR}$) [48] defined as:

$$\begin{align} \begin{split} v_\mathrm{CDR} & : = \frac{v_c}{v_d} \\ h_\mathrm{CDR} & : = \frac{h_c}{h_d} . \end{split} \end{align} \tag{ 4 }$$

The vCDR is, without doubt, the most popular and studied geometrical feature for diagnosing glaucoma on fundus images. Orlando et al [42] show that using only this measure it is possible to train a classifier that achieves quite competitive results on the REFUGE dataset. On the other hand, the European Glaucoma Society [44] discourages its use as the sole measure for diagnosis.

The datasets described in section 3.1 were used to train CNN models with different architectures. In this section, we present the methodology followed to train all these models. The results achieved are contained in section 4.1.

Four well–known CNN architectures were considered, VGG19 [49], ResNet50 [50], InceptionV3 [51], and Xception [52], available in the Keras module of the Tensorflow v2 package. The main motivation for using four different CNN architectures was to evaluate the influence of the specific network architecture on the results obtained. To adapt the models to the problem, we replaced the top layer of each network with a GlobalAveragePooling2D layer, followed by a DropOut layer and a Dense layer with 2 outputs using the SoftMax activation function. The input layer was set to 224 × 224 × 3 for ResNet50 and VGG19 and 299 × 299 × 3 for InceptionV3 and Xception. The DropOut value was set to 0.5 for VGG19 and 0.2 for InceptionV3, ResNet50, and Xception. Additionally, we kept the BatchNormalization layers in inference mode for InceptionV3, ResNet50, and Xception, to avoid updating the non-trainable weights of these layers during the training phase.

In all cases, to further ensure the robustness of our models, we employed a 5-fold cross-validation technique on the training set. We trained all the models using the same strategy. First, we froze the pre-trained base model and trained the new top layer for 200 epochs using an RMSprop optimizer with a learning rate of $1\times 10^{-6}$ and categorical cross-entropy loss. Second, we unfroze the base model and trained the entire model end-to-end for 250 epochs using the same optimizer, learning rate, and loss function. We used a batch size of 32 for all models, except for Xception, where we used a batch size of 16 due to memory constraints.

We applied the pre-processing function included in Keras for each of the models. For ResNet50 and VGG19, we transformed each image from RGB to BGR and centered each channel to zero by subtracting the mean values of the ImageNet dataset. For InceptionV3 and Xception, we scaled each image between −1 and 1.

To address class imbalance, we set class weights when training the models. The weight for each class was computed per fold using the formula in equation (5), where W_c is the weight for class c, N is the total number of data samples, and N_c is the number of samples of class c.

$$\begin{align} W_c = \frac{N}{2 * N_c} \end{align} \tag{ 5 }$$

To prevent overfitting, we applied data augmentation to the input samples, consisting of random rotations $(-15^{\circ}, 15^{\circ})$, horizontal flip and random brightness adjustments in the range [0.5, 1.5].

We selected the best-performing model for each fold per network architecture based on the epoch that maximized the validation accuracy averaged across the 5 folds. In total, we obtained 20 models per dataset, 5 different models per architecture.

The methodology used to evaluate how well the black box model f, in our case a neural network model, is approximated by the interpretable surrogate model g is the one proposed in [16]. Accordingly, the following steps have to be accomplished to obtain the surrogate model (see algorithm 1):

• 1.  
  Choose a dataset X, which is a matrix where each column represents one of the geometrical features outlined in section 3.2 and each row an image. These features are extracted from the set of images that was used for training the neural network models.
• 2.  
  Obtain predictions from the trained neural network models for the selected set of images. The ouput probability of being a glaucomatous eye was taken as the target for this problem.
• 3.  
  Choose an interpretable model. In this work, we have focused on linear regression and decision tree regression models. Note the need to pose the problem as a regression since we are considering a continuous variable as the target, as mentioned in the previous step. Linear regression was chosen because of its simplicity, while decision tree regression was chosen for its ability to provide good explanations and natural visualization through its nodes and edges.
• 4.  
  Train the chosen interpretable model on the dataset X and the predictions obtained in steps 1 and 2.

Once the surrogate model g is trained, we can evaluate how well it approximates f by using the R-squared measure [16]:

$$\begin{align} R^2 = 1-\frac{\mathrm{SSE}}{\mathrm{SST}} = 1-\frac{\sum_{i = 1}^n\left(\hat{y}_*^{\left(i\right)}-\hat{y}^{\left(i\right)}\right)^2}{\sum_{i = 1}^n\left(\hat{y}^{\left(i\right)}-\bar{\hat{y}}\right)^2}. \end{align} \tag{ 6 }$$

In the equation above, $\hat{y}_*^{(i)}$ represents the prediction of the interpretable model g for the ith instance, $\hat{y}^{(i)}$ represents the prediction of f, the black box model, for the ith instance, and $\bar{\hat{y}}$ represents the mean of the predictions of f. SSE refers to the sum of squares error, while SST refers to the sum of squares total. As Molnar points out [16], the R-squared measure provides insight into how much variance can be accounted for by the surrogate model. A low SSE and an R-squared value close to 1 indicate that the surrogate model is an accurate approximation of the black box model. Conversely, a high SSE and an R-squared value close to 0 suggest that the surrogate model is unable to explain the behavior of the black box model.

Following the procedure described in section 3.3, we performed a series of experiments using the different datasets detailed in section 3.1. Our own dataset (Dataset IR), as already mentioned, comprises 606 images and was randomly partitioned into training and testing sets with an 80/20 ratio. We evaluated the performance of the trained models on this test set, consisting of 41 samples from healthy subjects and 80 samples from glaucoma subjects. The results obtained are presented in table 3, which reports the worst and best performance per network architecture in terms of balanced accuracy. Balanced accuracy, which is the arithmetic mean of sensitivity and specificity, has been shown to be a good metric for evaluating classifiers with imbalanced classes [53]. Other commonly used metrics are also shown.

The simplified version of our dataset (Dataset C) showing only the disc and cup contours was used to train another 20 models under the same conditions as before and using the same hyper-parameters except for the data augmentation technique which did not include random brightness adjustments because it makes no sense for this type of images. The results, as presented in table 4, show a slight decrease in the performance of the models with respect to those trained with the original images but still can be considered as competitive despite not having access to all the image content. This observation supports our hypothesis that the geometry of the disc and cup carries crucial diagnostic information our neural network models effectively exploit. We also compared the performance of these models with that of two expert ophthalmologists who were asked to make a diagnosis based on the simplified images alone. Interestingly, the CNN models outperformed the experts in all cases, suggesting that the neural networks do a excellent job with this limited information.

The second dataset (Dataset XC) includes only simplified images and, as discussed in section 3.1, merges our data with RIM-ONE DL, REFUGE, and DRISHTI-GS. In this case, as also described in that section, our dataset was partitioned as in the first two experiments, while the remaining datasets were split according to their respective published splits for training and testing. As a result, 763 images from healthy subjects and 511 images from subjects with glaucoma were available for training. The training of these models was also carried out under the same conditions as described for the previous experiments. The results for the test set, consisting of 508 samples from healthy subjects and 210 samples from glaucoma subjects, are displayed in table 5. As can be seen on the table, in terms of balanced accuracy, the performance achieved by the models trained on the Dataset XC is generally higher than that of the models trained on the Dataset C alone.

Overall, the trained models performed well in all three scenarios evaluated. This is crucial as a poorly performing black box model renders the interpretation of the surrogate model pointless, making the black box model itself insignificant.

Figure 5 shows the plots of the average R-squared values obtained for each of the architectures, both for our image dataset (Dataset IR) and for its simplified version (Dataset C), also with two different types of surrogate models, linear regression and decision trees. Decision tree models were trained using grid search in their parameter space, exploring all posible combinations of maximum depths of 3, 4 and 5 and minimum number of samples per leaf node of 1, 2, 3, 4, 5, 8, 10 and 15. Both linear regression and decision tree models were trained using cross validation. The averaged results of the 3 folds for each architecture are shown. On the other hand, figure 6 shows an example scatter plot to compare the predictions of one of the folds for the VGG19 architecture and the corresponding surrogate decision tree for both the original and contour images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Important differences are observed between the models trained with the original images and those trained with the equivalent contour images. While with the original images, the average R-squared values range between 0.5 and 0.6, in the case of the contour images, the values are considerably higher, even exceeding 0.9. This seems to show that the surrogate system, trained with the geometric features, finds it easier to mimic the prediction of the neural networks when the information with which they are trained is of the same nature, being only the contours of the disc and the cup. This is also reflected in the example in figure 6. On the other hand, there are slight differences between using a linear regressor or a decision tree as a surrogate system. In the case of contour images, the decision tree provides better average R-squared values for all architectures, which seems logical since it is able to generate more general nonlinear approximations. However, for the models trained with the original images, the differences on average between the two types of models are very small for all architectures, even with an advantage in some cases for linear regression.

There is no clear criteria as to when a surrogate model should be considered adequate in order to be confident that it sufficiently approximates the black box model. It is not clear what the best cut-off point for R-square is, although some publications speak of 80% of the explained variance [16]. In any case, from the experiments described above it seems clear that in the problem addressed, the most valid option is the one represented by the surrogate models based on decision trees that approximate the predictions of the CNNs trained with the disc and cup contours. For this reason, the experiments that we will describe in the remainder have been raised in this scenario. In addition, to carry out a better evaluation, we have used the neural network models trained with the extended set (Dataset XC) as explained in section 4.1. Figure 7 shows the average values of R-squared obtained for these models. It is observed how they are even higher than those shown in figure 5. Figure 8 displays a scatter plot akin to that seen in figure 6, allowing us to compare the predictions from one of the folds for the VGG19 architecture with the corresponding surrogate decision tree for Dataset XC.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As already mentioned, decision trees are often used as surrogate models because of their good properties of simplicity, clarity and intuitiveness. Feature importance in decision trees refers to the calculation of the usefulness of each feature in predicting the target variable. Decision trees make splits in the data based on the values of the input features, and the importance of each feature is determined by how much these splits improve the purity of the resulting groups of data. In our case, the importance of each feature is determined by its contribution to the reduction of the mean squared error (MSE).

When an observation from the dataset is evaluated by the splitting rule, it either satisfies the rule and goes to the left of the node, or it does not satisfy the rule and goes to the right. This process continues until a terminal node is reached, at which point the prediction for the observation is made based on the values of the target variable in that group.

Let us denote each node as $N_m, m \in {1,..,k}$ where k is the number of nodes of the tree. Every node in the decision tree, except for the leaf nodes at the final depth, has two child nodes: a left child $N_m^l$ and a right child $N_m^r$.

Figures 9–12, show decision trees fitted to the predictions made by each architecture with the extended dataset (XC). Specifically, the trees corresponding to the fold for which the highest R-squared value has been obtained for each architecture are shown. The color of each node represents the predicted value. Lower values are represented by lighter colors, while higher values are represented by darker tones of orange. In these trees, each node contains some information in the following order (from top to bottom):

• 1.  
  A splitting rule that uses the feature f and a threshold value of t. If an observation's f value is less than or equal to t, it follows the path to the left of the node, and if it is greater than t, it follows the path to the right. Leaf nodes do not have a splitting rule.
• 2.  
  The mean squared error. The statistic that is used as the splitting criteria. It is calculated with the following equations:

  $$\begin{align} \bar{y}_m & = \frac{1}{n_m}\sum_{y\in Q_m}^ny \end{align} \tag{ 7 }$$

  $$\begin{align} H\left(N_m\right) & = \frac{1}{n_m}\sum_{y\in Q_m}^n\left(y-\bar{y}_m\right)^2 \end{align} \tag{ 8 }$$

  where Q_m is the data with n_m samples at node m.
• 3.  
  The samples value represents the number of observations in the training dataset that follows the path to this node.
• 4.  
  The predicted value for the node $\bar{y}_m$ equation (7). If an observation follows the path that leads to this node, that would be the predicted value for this obervation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An essential property of a node with child nodes is its node importance equation (9). The node importance is calculated using a specific formula, where the basic idea is that if the mean squared error in the child nodes is low, then the importance of the node, and particularly the feature used in its splitting rule, is high. In other words, the lower the MSE in the children, the greater the importance of the node and the feature used to create the split.

$$\begin{align} U\left(N_m\right) = W\left(N_m\right)H\left(N_m\right)-W\left(N_m^l\right)H\left(N_m^l\right) - W\left(N_m^r\right)H\left(N_m^r\right) \end{align} \tag{ 9 }$$

where the weight $W(N_m)\in(0,1]$ is determined for every node in the tree and is found by dividing the number of observations in that node by the total number of observations in the dataset.

The importance of each feature is calculated by summing the importance of each node that uses that feature as a splitting criterion. Figure 13 shows the feature importances obtained with the trained decision trees for each architecture and fold, for the Dataset XC. Thus, what we observe is the relevance of each feature in the trained decision tree model to mimic the predictions of the corresponding CNN, i.e. the probability that the input image corresponds to an eye with glaucoma. In this regard, it is quite clear that there is one feature that stands out above the others and that is the vCDR. We have already referred to it in section 3.2, as it is a very popular measure for the diagnosis of glaucoma. The importance value assigned to it in the four network architectures considered seems to corroborate this popularity. It does not mean that the neural network is calculating the cup-to-disc ratio as such, but rather that it plays a very prominent role in the decisions made by the surrogate model, perhaps because the network learns to calculate features that may be strongly correlated with it. However, this is a somewhat controversial measure, as some authors point it out as not too reliable for the diagnosis of this pathology [44, 45, 54]. Jonas et al [45] note that the interindividual variability of cup-to-disc ratios and their dependence on optic disc size should also be taken into account. The cup-to-disc ratio is a relative measure independent of cup and disc size, and, also, of the magnification of the camera used but the optic disc size is an absolute measure. In fact, although to a much lesser extent, figure 13 shows other features that contribute to model prediction which are absolute measures such as the inferior and superior rim widths given by I-ISNT and S-ISNT, and the size of the optic cup, OC size. Some other features such as OC shape and NR width also contribute globally to the predictions but still to a lesser extent. Table 6 shows the five most important features for the surrogate decision trees per architecture, on average. Kumar et al [54] conclude in their article that the rim-to-disc ratio outperforms cup-to-disc ratio for glaucoma prescreening. However, according to the results found, in none of the 20 surrogate models considered does this characteristic seem to be of relevant importance.

Zoom In Zoom Out Reset image size

Apart from this global analysis of the importance of the features, it is also interesting to do a more local analysis based on the representation of the specific decision trees shown.

There are common patterns that are repeated in the trees of all architectures and are used in most of the training set samples. One common pattern that can be found in all trees is the use of a split at the root node which implies that the vCDR is used around a value of 0.6. This generates a left sub-tree in which the lower values are predicted and a right sub-tree in which the higher values are predicted, in general.

Within the left subtrees of all architectures, a large number of samples in the training set with low glaucoma probability values follow a pathway in which the I-ISNT, S-ISNT, and OC shape features are used. The order in which these features are used varies by architecture, but their use appears consistent across all architectures for most of the clearly negative glaucoma cases.

In the subtrees on the right, the vast majority of samples with a high probability value for glaucoma follow a trajectory marked by the use of V CDR.

There are also specific patterns with which each architecture defines paths in the tree to process a smaller set of samples, for which the probability of glaucoma has intermediate values.

The tree corresponding to the InceptionV3 architecture (figure 9) uses OC Size, OD Size and NR Width to decrease the predicted value of samples traversing the right subtree, and OC Shape and RD Ratio to slightly increase the predicted value of samples traversing its left subtree.

The tree corresponding to the VGG19 architecture (figure 10) uses the OC Shape and OC Size features in the opposite way, i.e. it uses OC Shape to decrease the predicted value in the right subtree and OC Size to increase it in the left subtree.

The tree corresponding to the ResNet50 architecture (figure 11) prefers to use NR Size to increase the predicted value of samples traversing its left subtree instead of RD Ratio, and uses H CDR instead of OC Size to decrease the predicted value of samples traversing its right subtree.

Finally, the tree corresponding to the Xception architecture (figure 12) mainly uses OC Shape in both subtrees to adjust the predicted value. In the case of the left subtree, a threshold of 0.82 is used in this feature to increase the predicted value for samples that are below it. For the right subtree, a threshold of 0.87 is used to reduce the predicted value of the samples that exceed it.

As mentioned in the introduction, despite their popularity, some authors have expressed doubts regarding the reliability and robustness of the saliency maps for interpreting CNN models in the field of medicine [11, 14]. One of the issues they have highlighted is the poor correlation between saliency maps corresponding to models with different architectures.

To go a little deeper into this question, we have performed some experiments with our own disc and cup contour images (Dataset C), and with the 20 models trained with them described in section 3. In particular, we have calculated the degree of correlation between saliency maps corresponding to the same images but obtained from CNN models with different architectures, following the same methodology as in [55]. The saliency maps used were generated with well-known gradient backpropagation, integrated gradient, occlusion and CAM-based methods [56]. Figure 14 shows the results obtained where each boxplot has been calculated from the considered saliency methods between a pair of network architectures for the same fold. It is observed that on average the correlation is quite low in all possible combinations of architectures.

Zoom In Zoom Out Reset image size

On the other hand, figure 15 shows some saliency maps computed with the occlusion method [57] corresponding to three example images for the models of different network architectures trained with the Dataset XC referred to in the previous section. Apart from the extreme difficulty of interpreting what the networks may be doing with the contour information to make the decision, it is possible to see the important differences, in some cases, between the saliency maps.

Zoom In Zoom Out Reset image size

Therefore, in view of these results, it is worth considering the usefulness of this interpretability technique in a problem such as the diagnosis of glaucoma. Unlike other medical problems in which there may be clear and well-located biomarkers of the pathology under consideration, in the case of glaucoma we find ourselves in a quite different scenario. In accordance with what was seen in the previous section, even assuming that the neural networks may be calculating features correlated with indicators commonly used by specialists such as the vCDR, it is practically impossible to deduce this fact from the simple observation of the saliency maps. As already mentioned in section 1, the usefulness of trying to explain at the pixel level what the network may be doing to generate its predictions is questionable. In this sense, the explored alternative of using a surrogate system as a means of interpretation seems to be more useful in a problem of this nature, also offering a more coherent result between the different architectures analyzed.

So far we have seen how decision trees are able to reproduce the predictions of CNNs providing a certain interpretability mechanism for their decisions. However, as Rudin points out [6], rather than producing explanations faithful to the original model, these surrogate models show trends in the way predictions are related to features. Thus, while they may provide useful and more interpretable information than techniques such as saliency methods, as discussed in the previous section, they are not an entirely satisfactory solution.

For this reason, we decided to go a step further and evaluate on the real outcome, i.e. the original class test labels, the performance of the decision tree acting as a surrogate model that has learned to make predictions in a similar way as CNN does. Note that the surrogate model itself never sees the real outcome but the predictions of the neural networks. We would therefore be talking about a form of knowledge distillation in which we transfer knowledge from the network to the decision tree. Since the surrogate models have been trained for regression to mimic the probability value of being a glaucomatous eye, to do the classification we have simply thresholded this probability value at 0.5 to determine whether or not it is glaucoma. With this in mind, table 7 shows the balanced accuracy values found, on the Dataset XC, for the surrogate decision trees corresponding to the 20 network models considered, along with the value obtained for a decision tree trained without surrogation using the real outcome. Figure 16 visualizes the decision tree for this case in which a pattern similar to those we had already discussed is observed, with the vCDR playing a major role as the discriminant feature.

Zoom In Zoom Out Reset image size

Table 7 shows that in some cases the knowledge distillation process produces decision tree models with better results. This finding is very interesting not only because of the improvement in itself, but also because we have managed to obtain intrinsically interpretable models with performance equivalent to CNN models, thus overcoming the limitations in the explanations referred to earlier in this section.