Interpreting Deep Neural Networks with the Package innsight
The R package innsight offers a general toolbox for revealing variable-wise
interpretations of deep neural networks' predictions with so-called feature
attribution methods. Aside from the unified and user-friendly framework, the
package stands out in three ways: It is generally the first R package
implementing feature attribution methods for neural networks. Secondly, it
operates independently of the deep learning library allowing the interpretation
of models from any R package, including keras, torch, neuralnet, and even
custom models. Despite its flexibility, innsight benefits internally from the
torch package's fast and efficient array calculations, which builds on LibTorch
$-$ PyTorch's C++ backend $-$ without a Python dependency. Finally, it offers a
variety of visualization tools for tabular, signal, image data or a combination
of these. Additionally, the plots can be rendered interactively using the
plotly package.
Niklas Koenen, Marvin N. Wright
# Interpreting Deep Neural Networks with the Package innsight

###### Abstract

The \(\mathsf{R}\) package **innsight** offers a general toolbox for revealing variable-wise interpretations of deep neural networks' predictions with so-called feature attribution methods. Aside from the unified and user-friendly framework, the package stands out in three ways: It is generally the first \(\mathsf{R}\) package implementing feature attribution methods for neural networks. Secondly, it operates independently of the deep learning library allowing the interpretation of models from any \(\mathsf{R}\) package, including **keras**, **torch**, **neuralnet**, and even custom models. Despite its flexibility, **innsight** benefits internally from the **torch** package's fast and efficient array calculations, which builds on **LibTorch** - **PyTorch**'s \(\mathsf{C++}\) backend - without a Python dependency. Finally, it offers a variety of visualization tools for tabular, signal, image data or a combination of these. Additionally, the plots can be rendered interactively using the **plotly** package.

_Keywords_: neural networks, feature attribution, interpretable machine learning, explainable, artificial intelligence, XAI, IML, **torch**, **keras**, \(\mathsf{R}\).

## 1 Introduction

Throughout the past decade, neural networks have unleashed a tremendous surge of attention and infiltrated almost all conceivable domains of science, industry, and public life. Mainly, their increasing popularity is due to their natural ability to extract patterns and knowledge from vast amounts of structured raw data thanks to modern computing capacities and deliver outstanding performance (Krizhevsky _et al._, 2017; LeCun _et al._, 2015; Silver _et al._, 2016; Bengio _et al._, 2021). However, the intelligently learned decision-making process of the neural network remains inscrutable and hidden from the user due to its enormous complexity. Interpretations cannot be inferred as straightforward from this so-called _black box_ as, for example, the coefficients of linear models. As a consequence, the gain in predictive accuracy and model flexibility generally comes at the price of an increasingly opaque and intricate machine learning model, as was already noted by Gunning and Aha (2019). Nevertheless, it is precisely this question of interpretability - or, informally speaking: Why did a network make a certain prediction? - that is becoming more and more relevant for applications with high-stake decisions and possibly becoming a legal requirement, e.g., in autonomous systems (O'Sullivan _et al._, 2019), healthcare (Schneeberger _et al._, 2020) or data processing in general (European Union, 2016;
###### Abstract

We propose a novel approach to the classification of the joint and temporal
is produced in the standard forward pass and the class to be explained is selected. In the subsequent method-specific backward pass, each input variable is assigned a relevance score to the chosen output class and can be visualized in a heat map or bar chart depending on the input type. An overview of the most popular feature attribution methods can be found in Section 2.

To make feature attribution methods widely accessible to users and to provide them with a unified and easy-to-use interface, several software packages have been developed in the last few years. The most popular packages are **investigate**(Alber _et al._, 2019) for the deep learning library **Keras**(Chollet _et al._, 2015), and **captum**(Kokhlikyan _et al._, 2020) and **zennit**(Anders _et al._, 2021) for **PyTorch**(Paszke _et al._, 2019). However, all these packages are Python-exclusive and only support networks of specific deep learning libraries. In order to make feature attribution methods easily available to the R community, we contribute the software package **innsight** pursuing the following goals:

* _First feature attribution R package_: **innsight** is the first R package that implements the most popular feature attribution methods for neural networks unified in a single user-friendly package.
* _Computationally efficient_: The powerful **torch**(Falbel and Luraschi, 2023) package is utilized internally for all calculations, which builds on **LibTorch**, the C++ variant of **PyTorch**(Paszke _et al._, 2019), and doesn't rely on a Python dependency.
* _Deep-learning-library-agnostic_: The passed trained models are not limited to a specific deep learning library. The package supports models from the R packages **keras**(Kalinowski _et al._, 2023), **torch** and **neuralnet**(Gunther and Fritsch, 2010), but under some constraints, an arbitrary model can be passed as a list to be fully flexible.
* _Visualization tools_: **innsight** offers several visualization methods for individual or summarized results regardless of whether it is tabular, signal, image data or a combination of these. Additionally, interactive plots can be created based on the **plotly** package (Sievert _et al._, 2022).

The **innsight** package is available from the Comprehensive R Archive Network (CRAN) at [https://CRAN.R-project.org/package=innsight](https://CRAN.R-project.org/package=innsight) or from our GitHub repository at [https://github.com/bips-hb/innsight/](https://github.com/bips-hb/innsight/).

Figure 1: General procedure of feature attribution methods: First, an input instance \(\mathbf{x}\) flows through the model \(f\) to obtain a prediction \(\mathbf{\hat{y}}\). Then, the desired output node or class \(\hat{y}_{c}\) to be explained is selected. Finally, the relevance \(R^{c}_{i}\) of the individual input variables \(i\) at the selected output \(c\) is calculated in a backward pass.
The rest of the paper is structured as follows: First, we overview the most popular feature attribution methods for neural networks in Section 2. Then, in Section 3, we elaborate on the package's design, functionality, and capabilities. Next, the package is applied to a basic example on a penguin dataset and an advanced example for melanoma detection based on image and tabular data as input types. In the concluding Section 5, the obtained package's outputs are compared and validated with the already mentioned Python equivalents.

## 2 Methodology of feature attribution

Feature attribution methods for neural networks are a group of local interpretation methods that assign to each input variable the contribution or impact to a chosen model output. For example, suppose an input instance \(\mathbf{x}\in\mathbb{R}^{p}\) with \(p\in\mathbb{N}\) variables is fed forward through a neural network \(f:\mathbb{R}^{p}\to\mathbb{R}^{C}\) resulting in an output \(f(\mathbf{x})=\mathbf{\hat{y}}\in\mathbb{R}^{C}\) with \(C\in\mathbb{N}\) classes or regression outputs. In this case, a feature attribution method assigns relevance scores \(R^{c}_{1},\ldots,R^{c}_{p}\) to each of the input features \(x_{1},\ldots,x_{p}\) of \(\mathbf{x}\) on a chosen output class or node \(\hat{y}_{c}\) of the prediction \(\mathbf{\hat{y}}\) to be explained as already described in Figure 1.

### Gradient-based methods

Gradient-based methods are the fastest and most straightforward interpretation methods because they operate on the default techniques of the high-level deep learning libraries for computing gradients during the training loop. However, these methods - in a sense - calculate the derivatives of chosen output to the input variables instead of the derivatives of the loss value to the model parameters during gradient descent. Although the terminology of gradient-based methods can often be interpreted more broadly, we only consider techniques that use the default gradient methods. For example, Ancona _et al._ (2018) showed that variants of the _layerwise relevance propagation (LRP)_ and _deep learning important features (DeepLift)_, discussed later in Section 2.2 and 2.3, can approximately be considered gradient-based. Regardless, they are not mentioned in this section, since, on the one hand, the fundamental purpose of these variants is different, and, on the other hand, they require an overwriting of the standard gradient methods.

#### Vanilla Gradient

One of the first and most intuitive methods for interpreting neural networks is the _Gradient_ method introduced by Simonyan _et al._ (2014), also known as _Vanilla Gradients_ or _Saliency maps_. This method computes the gradients of the selected output with respect to the input variables. Therefore, the resulting relevance values indicate those variables that can be locally perturbed the least to change the outcome the most. Mathematically, this method can be described by the following formula for the input variable \(x_{i}\) with \(\mathbf{x}\in\mathbb{R}^{p}\), the model \(f:\mathbb{R}^{p}\to\mathbb{R}^{C}\) and the output \(\hat{y}_{c}=f(\mathbf{x})_{c}\):

\[\text{Gradient}(\mathbf{x})_{i}^{c}=\frac{\partial\,f(\mathbf{x})_{c}}{\partial\,x_{ i}}=\frac{\partial\,\hat{y}_{c}}{\partial\,x_{i}}.\]

Assuming that the function behaves linearly overall, increasing \(x_{i}\) by one raises the output by \(\text{Gradient}(\mathbf{x})_{i}^{c}\). In general, however, neural networks are highly nonlinear, so this interpretation is only valid for very small changes of \(x_{i}\).
### Smoothed gradients

The _smoothed gradients (SmoothGrad)_ method, introduced by Smilkov _et al._ (2017), addresses a significant problem of the basic Gradient method. As described in the previous subsection, gradients locally assume a linear behavior, but this is generally no longer the case for deep neural networks with many layers. These have large fluctuations and abruptly change their gradients, making the interpretations of the Gradient method noisier and potentially misleading. Smilkov _et al._ (2017) proposed to compute the gradients of randomly perturbed copies of \(x_{i}\) and determine the average gradient from that instead of calculating only the gradient in \(x_{i}\). For obtaining relevance values with the SmoothGrad method for the individual components \(x_{i}\in\mathbb{R}\) of an instance \(\mathbf{x}\in\mathbb{R}^{p}\), first \(K\in\mathbb{N}\) realizations of a \(p\)-dimensional multivariate Gaussian distribution \(q=\mathcal{N}(0,\sigma^{2}\mathrm{I}_{p})\) are generated describing the random perturbations, i.e., \(\mathbf{\varepsilon}^{1},\ldots,\mathbf{\varepsilon}^{K}\sim q\). Then the empirical mean of the gradients for variable \(x_{i}\) and output index \(c\) are calculated as follows:

\[\mathrm{SmoothGrad}(\mathbf{x})_{i}^{c}=\frac{1}{K}\sum_{j=1}^{K}\frac{\partial\, f(\mathbf{x}+\mathbf{\varepsilon}^{j})_{c}}{\partial\,x_{i}+\varepsilon_{i}^{j}}\approx \mathbb{E}_{\mathbf{\varepsilon}\sim q}\left[\frac{\partial\,f(\mathbf{x}+\mathbf{ \varepsilon})_{c}}{\partial\,x_{i}+\varepsilon_{i}}\right].\]

The number of perturbations \(K\) and the variance \(\sigma^{2}\) are hyperparameters. With the value of \(K\), the estimation accuracy for the mean gradient can be increased, but this goes hand in hand with a higher computational effort. The second parameter \(\sigma^{2}\) is mostly specified indirectly via a noise level \(\lambda\geq 0\) determining the percentage of the total range of the input domain that is covered by the standard deviation \(\sigma\), i.e., \(\lambda=\frac{\sigma}{x_{\text{max}}-x_{\text{min}}}\). Especially for images, this argument can be used to control the visual smoothness of the explanation.

### Gradient\(\times\)Input and SmoothGrad\(\times\)Input

A simple modification can change both previously discussed methods to the methods _Gradient\(\times\)Input_ and _SmoothGrad\(\times\)Input_. The gradients are calculated as in the respective sections and then multiplied by the corresponding input instance. The _Gradient\(\times\)Input_ method was introduced by Shrikumar _et al._ (2017b) and relies on a well-grounded mathematical background despite its simple idea: The basic concept is decomposing the output prediction \(\hat{y}_{c}\) according to its relevance to each input variable \(x_{i}\), i.e., into variable-wise additive effects

\[\hat{y}_{c}=f(\mathbf{x})_{c}=\sum_{i=1}^{p}R_{i}. \tag{1}\]

Mathematically, this method is based on the first-order Taylor decomposition. Assuming that a function \(g:\mathbb{R}^{p}\rightarrow\mathbb{R}\) is continuously differentiable in \(\mathbf{x}\in\mathbb{R}^{p}\), a remainder term \(\varepsilon(g,\mathbf{z},\mathbf{x}):\mathbb{R}^{p}\rightarrow\mathbb{R}\) with \(\lim_{\mathbf{z}\rightarrow\mathbf{x}}\varepsilon(g,\mathbf{z},\mathbf{x})=0\) exists such that

\[g(\mathbf{z}) =g(\mathbf{x})+\nabla_{\mathbf{x}}g(\mathbf{x})\cdot(\mathbf{z}-\mathbf{x})^{\top}+ \varepsilon(g,\mathbf{z},\mathbf{x})\] \[=g(\mathbf{x})+\sum_{i=1}^{p}\frac{\partial\,g(\mathbf{x})}{\partial\,x _{i}}(z_{i}-x_{i})+\varepsilon(g,\mathbf{z},\mathbf{x}),\quad\mathbf{z}\in\mathbb{R}^{p}.\]

The first-order Taylor formula thus describes a linear approximation of the function \(g\) at the point \(\mathbf{x}\) since only the first derivatives are considered. Consequently, a highly nonlinear and
continuous function \(g\) is well approximated only in a small neighborhood around \(\mathbf{x}\). For larger distances from \(\mathbf{x}\), sufficient small values of the residual term are not guaranteed anymore. The Gradient\(\times\)Input method considers the data point \(\mathbf{x}\) and sets \(\mathbf{z}=\mathbf{0}\). In addition, the residual term \(\varepsilon(f_{c},\mathbf{x},\mathbf{0})\) and the summand \(f(\mathbf{0})_{c}\) are ignored, which then results in the following approximation of the prediction \(f(\mathbf{x})_{c}\) in variable-wise relevances:

\[f(\mathbf{x})_{c}\approx\sum_{i=1}^{p}\frac{\partial\,f(\mathbf{x})_{c}}{ \partial\,x_{i}}\cdot x_{i},\quad\text{hence}\] \[\text{Gradient}\times\text{Input}(\mathbf{x})_{i}^{c}=\frac{ \partial\,f(\mathbf{x})_{c}}{\partial\,x_{i}}\cdot x_{i}.\]

Analogously, this multiplication can be applied to all gradients in the summation of the SmoothGrad method in order to compensate for local fluctuations:

\[\text{SmoothGrad}\times\text{Input}(\mathbf{x})_{i}^{c}=\frac{1}{K}\sum_{j=1}^{K} \frac{\partial\,f(\mathbf{x}+\mathbf{\varepsilon}^{j})_{c}}{\partial\,x_{i}+\varepsilon _{i}^{j}}\cdot(x_{i}+\varepsilon_{i}^{j}),\quad\mathbf{\varepsilon}^{1},\ldots, \mathbf{\varepsilon}^{K}\sim\mathcal{N}(0,\sigma^{2}\text{I}_{p}).\]

### Layer-wise relevance propagation (LRP)

The _layer-wise relevance propagation (LRP)_ method was introduced by Bach _et al._ (2015) and has a similar goal as the Gradient\(\times\)Input approach explained in the previous section: decomposing the output into variable-wise relevances conforming to Equation 1. The distinguishing aspect is that the prediction \(\hat{g}_{c}\) is redistributed layer by layer from the output node back to the inputs according to the layer's weights and intermediate values. The entire procedure is accomplished by rule-based relevance messages defining how to redistribute the upper-layer relevance to the lower layer. Consequently, the LRP method is applied to a model with \(L\) layers as follows, irrespective of the selected rule: As a first step, the model output score \(\hat{g}_{c}\) is used as the upper-layer relevance \(R_{1}^{L}\) to calculate the relevances for the layer of index \(L-1\) before the output layer via the selected relevance message. After that, the obtained values \(R_{i}^{L-1}\) are used as the upper-layer relevances for the layer preceding that layer of index \(L-1\), and so on. Repeating this procedure layer by layer until the input layer is reached, yields the relevance values \(R_{i}^{1}\) for each input variable \(i\) of the model input instance and, thus, the results of the LRP method. A visual overview of these steps is given in Figure 2.

In the following, we briefly overview the most popular variations of relevance messages flowing from a node of index \(j\) in layer \(l+1\) to node \(i\) in the preceding layer:

* **The simple rule:** The fundamental rule on which all other variations of relevance messages are more or less based is the _simple rule_ (also known as _LRP-0_). The relevances are redistributed to the lower layers according to the ratio between local and

Figure 2: Backward pass of the LRP method.
global pre-activation. Let \(\mathbf{x}\) the inputs of the preceding layer, \(\mathbf{w}\) the weight matrix and \(\mathbf{b}\) the bias vector of layer \(l\), and \(R_{j}^{l+1}\) the upper-layer relevance; then \(x_{i}\,w_{i,j}\) is the local and \(z_{i}=b_{j}+\sum_{k}x_{k}\,w_{k,j}\) the global pre-activation defining the simple rule as \[r_{i\gets j}^{(l,\,l+1)}=\frac{x_{i}\,w_{i,j}}{z_{j}}\,R_{j}^{l+1}.\]
* **The \(\varepsilon\)-rule:** One issue with the simple rule is that it is numerically unstable when the global pre-activation \(z_{j}\) vanishes and causes a division by zero. The \(\varepsilon\)_-rule_ (also known as _LRP-\(\varepsilon\)_) tackles those situations by adding a stabilizer \(\varepsilon>0\) that moves the denominator away from zero, i.e., \[r_{i\gets j}^{(l,\,l+1)}=\frac{x_{i}\,w_{i,j}}{z_{j}+\text{sign}(z_{j})\, \varepsilon}\,R_{j}^{l+1}.\] This inserted value \(\varepsilon\) absorbs some of the relevance and can, therefore, be utilized to achieve sparser and less noisy results for the explanation. As \(\varepsilon\) increases, a greater portion of the relevance is intercepted, sustaining only the most salient relevances for this relevance message.

Both variants have in common that they distribute the upper-layer relevance proportionally downward regarding the local and global pre-activations, i.e., \(x_{i}w_{i,j}\) and \(z_{j}\). Even though the \(\varepsilon\)-rule avoids division by zero, numerical inconsistencies can occur in both variants for very deep models. Since the pre-activations are not necessarily guaranteed to be positive, the local pre-activations may take on substantial positive or negative values that cancel out in the global pre-activation leading to magnified values in the preceding layer. As a result, larger relevances in the lower layers potentially accumulate in deep models and increasingly reach the limits of computational representation of floating point numbers. To prevent this blow-up of relevances, the authors introduced the \(\alpha\)_-\(\beta\)-rule_, which treats positive and negative pre-activations separately:

* **The \(\alpha\)-\(\beta\)-rule:** The \(\alpha\)-\(\beta\)_-rule_ was introduced to avoid numerical instabilities and enable a weighting between positive and negative relevances depending on the user's focus. This relevance message applies the simple rule to the positive and negative parts of the pre-activations, respectively, and takes the weighted sum of both. The weighting can be regulated by the hyperparameters \(\alpha,\beta\in\mathbb{R}\) satisfying \(\alpha+\beta=1\). Mathematically formulated, the rule is defined as follows: \[r_{i\gets j}^{(l,\,l+1)}=\left(\alpha\frac{(x_{i}\,w_{i,j})^{+}}{z_{j}^{+ }}+\beta\,\frac{(x_{i}\,w_{i,j})^{-}}{z_{j}^{-}}\right)\,R_{j}^{l+1}\] with \[z_{j}^{\pm}=(b_{j})^{\pm}+\sum_{k}(x_{k}\,w_{k,j})^{\pm},\quad(\cdot)^{+}= \max(\cdot,0),\quad(\cdot)^{-}=\min(\cdot,0).\]

For any of the rules described above, the relevance of the lower-layer node \(R_{i}^{l}\) is determined by summing up all incoming relevance messages \(r_{i\gets j}^{(l,\,l+1)}\) into the respective node of index \(i\), i.e.,

\[R_{i}^{l}=\sum_{j}r_{i\gets j}^{(l,\,l+1)}. \tag{2}\]
Since the bias vector is included in the computation of the global pre-activations in all presented variants, this term absorbs a certain amount of the upper-layer relevance. Consequently, the LRP methods approximate the output prediction rather than providing an accurate representation of the targeted decomposition in Equation 1.

There are even more variants of relevance messages discussed in the literature suitable for various situations or layer types: For example, the _deep Taylor decomposition_ (also called \(z^{+}\)-_rule_) in ReLU models - also achieved with the \(\alpha\)-\(\beta\)-rule with \(\alpha=1\) - allows filtering out only positive relevances (Montavon _et al._, 2017), or the _\(\gamma\)-rule_ favoring positive over negative relevances (Montavon _et al._, 2019). Moreover, some rules are specifically designed for the input layer (Montavon _et al._, 2017). Due to the rule independence of how the lower-layer relevances are computed from the relevance messages in Equation 2, the rules can also be set individually for each layer, called _composite-rule_(Montavon _et al._, 2019; Kohlbrenner _et al._, 2020).

### Deep learning important features (DeepLift)

One method that, to some extent, echoes the idea of LRP is the so-called _deep learning important features (DeepLift)_ method introduced by Shrikumar _et al._ (2017). It behaves similarly to LRP in a layer-by-layer backpropagation fashion from a selected output node back to the input variables considering the simple rule. However, it incorporates a reference value \(\mathbf{\tilde{x}}\) to compare the relevances with each other. Hence, the relevances of DeepLift represent the relative effect of the outputs of the instance to be explained \(f(\mathbf{x})_{c}\) and the output of the reference value \(f(\mathbf{\tilde{x}})_{c}\). By taking the difference, the bias term is eliminated in the relevance messages, preventing the relevance absorption and leading to an exact variable-wise decomposition of the difference-from-reference output \(\Delta\hat{g}_{c}=f(\mathbf{x})_{c}-f(\mathbf{\tilde{x}})_{x}\), i.e.,

\[\Delta\hat{g}_{c}=f(\mathbf{x})_{c}-f(\mathbf{\tilde{x}})_{x}=\sum_{i=1}^{d}R_{i}.\]

In contrast to the LRP method, DeepLift defines a multiplier layer by layer, starting from the output layer and propagating to the input layer instead of directly determining the relevances in each intermediate stage. Based on these multipliers, the contribution of an arbitrary variable to the difference-from-reference output can be obtained by multiplying it by the corresponding difference-from-reference input. For an arbitrary layer with the layer's input \(\mathbf{x}\), reference input \(\mathbf{\tilde{x}}\) and multiplier \(m_{\Delta\mathbf{x}\Delta\hat{g}_{c}}\), this means:

\[\sum_{i}m_{\Delta x_{i}\Delta\hat{g}_{c}}\left(x_{i}-\tilde{x}_{i}\right)=m_{ \Delta\mathbf{x}\Delta\hat{g}_{c}}\cdot\left(\Delta\mathbf{x}\right)^{\top}=\Delta \hat{g}_{c}. \tag{3}\]

The multipliers fulfill a chain rule allowing the computation of the multiplier for the preceding layer given the already calculated one \(m_{\Delta\mathbf{t},\Delta\hat{g}_{c}}\), i.e.,

\[m_{\Delta x_{i}\Delta\hat{g}_{c}}=\sum_{j}m_{\Delta x_{i}\Delta t_{j}}\,m_{ \Delta t_{j}\Delta\hat{g}_{c}}. \tag{4}\]

In other words, the chain rule justifies defining the multipliers for each layer or part of a layer separately before combining them with the upper-layer multipliers. The authors distinguish between the linear and nonlinear components of a layer and provide definitions of the multipliers for each of them:
* **Linear rule:** For the linear components of a layer, such as matrix multiplication in dense or convolution layers, the weights of the corresponding layer are used as the multipliers, i.e., \(m_{\Delta x_{i}\Delta z_{j}}=w_{i,j}\).
* **Rescale rule:** This rule can be used for all nonlinear parts of a layer that can be reduced to a one-dimensional function \(\sigma\), e.g., all activations such as ReLU, tanh, or sigmoid. In this case, the ratio between the difference-from-reference activation \(\Delta\sigma(z)_{j}=\sigma(z_{j})-\sigma(\tilde{z}_{j})\) and the pre-activation \(\Delta z_{j}=z_{j}-\tilde{z}_{j}\) gives the multiplier, i.e., \(m_{\Delta z_{j}\Delta\sigma(z)_{j}}=\frac{\Delta\sigma(z)_{j}}{\Delta z_{j}}\). To avoid numerical instability caused by a vanishing denominator, the gradient of \(\sigma\) at \(z_{j}\) is used instead of the multiplier when \(z_{j}\) is close to its reference value \(\tilde{z}_{j}\).
* **RevealCancel rule:** This rule is designed for non-linearities \(\sigma\) to propagate meaningful relevances for saturated activations and discontinuous gradients through the layers' activation part, even when activations like ReLU eliminate the values. Similar to the normal pre-activations in the \(\alpha\)-\(\beta\)-rule for LRP, the positive \(\Delta z_{j}^{+}\) and negative \(\Delta z_{j}^{-}\) difference-from-reference pre-activations are considered separately, ensuring the propagation of expressive contribution scores. Descritively, the _RevealCancel_ rule can be explained in a way that the multiplier for the positive part \(m_{\Delta z_{j}^{+}\Delta y_{j}^{+}}\) is the ratio between the average effect of \(\Delta z_{j}^{+}\) after activating, before and after the negative part \(\Delta z_{j}^{-}\) has been added, and the positive difference-from-reference pre-activation \(\Delta z_{j}^{+}\). In the same way, the negative multiplier \(m_{\Delta z_{j}^{-}\Delta y_{j}^{-}}\) is given by the ratio of the average impact of \(\Delta z_{j}^{-}\) after activating, before and after the positive part \(\Delta z_{j}^{+}\) has been added, to \(\Delta z_{j}^{-}\). Mathematically, the rule is defined as \[m_{\Delta z_{j}^{\pm}\Delta y_{j}^{\pm}}=\frac{\frac{1}{2}\left(\sigma(\tilde{ z}_{j}+\Delta z_{j}^{\pm})-\sigma(\tilde{z}_{j})+\sigma(\tilde{z}_{j}+\Delta z_{j}^{ \pm}+\Delta z_{j}^{\mp})-\sigma(\tilde{z}_{j}+\Delta z_{j}^{\pm})\right)}{ \Delta z_{j}^{\pm}}.\]

These rules, along with the chain rule (Equations 3-4), enable the successive computation of the input variables' contributions \(R_{i}\) to the difference-from-reference output \(\Delta\hat{g}_{c}\) in a single backward pass.

The reference value is the only crucial hyperparameter for the DeepLift method, apart from the rule for non-linearities. This choice depends significantly on the application and usually requires proficient domain-specific knowledge. Nevertheless, the authors suggest asking oneself the question what one wants to measure an effect against. For example, taking the background color or blurred versions of the original picture as the reference values for images are reasonable choices. In many cases, zeros as a baseline are also used. Ancona _et al._ (2018) showed that using the Rescale rule with activations crossing the origin (i.e., \(\sigma(0)=0\)) and a zero baseline as reference value \(\tilde{x}\) coincides with the Gradient\(\times\)Input method discussed in Section 2.1.3.

### Connection weights

One of the earliest methods specifically designed for neural networks is the _connection weights (CW)_ method invented by Olden _et al._ (2004), resulting in a global relevance score for each input variable. The basic idea of this approach is to multiply all path weights for each possible connection between an input variable \(x_{i}\) and the output node or class \(\hat{y}_{c}\) and then calculate
the sum of all of them. However, this method ignores all bias vectors and all activation functions during calculation. Analogously to the previous methods, CW can also be defined layer by layer, deriving the relevance for layer \(l\) from the upper layer as follows:

\[R_{i}^{l}=\sum_{j}w_{i,j}R_{j}^{l+1}.\]

Since only the model weights are used, this method is independent of input data and, thus, a global interpretation method. Inspired by the method Gradient\(\times\)Input (see Sec. 2.1.3), it can also be extended into a local method by taking the point-wise product of the global CW method and the input data.

## 3 Functionality and usage

The R package **innsight** combines all the methods discussed in the previous section in a user-friendly structure and a unified step-based workflow from the trained model to the visualization of the relevances of a feature attribution method. For efficient high-dimensional array calculations, the package uses the R package **torch**[18], which builds on **LibTorch** (the C++ variant of **PyTorch**[17]), and consequently runs without a Python dependency (see Fig. 3). The following three steps yield the requested results regardless of the class of the passed model or the chosen feature attribution method:

* Step 1: Convert the model
* Step 2: Apply selected method
* Step 3: Get or visualize results.

Apart from the utilized packages for the internal workflows, calculations, and visualizations discussed in the following sections, the packages **checkmate**[10] and **cli**[11] are generally used for all argument verifications, internal checks, and terminal outputs of messages, warnings, and errors.

### Step 1 - Convert the model

The key step that turns the **innsight** package into a deep-learning-library-agnostic approach and unlocks the provided **torch** toolbox to all methods is this first step, which essentially analyzes a passed model and creates a **torch**-based replication. For the user, however, the internal processes remain hidden, and the entire conversion step is accomplished by creating a new instance of the class Converter:

Converter$new(model,  input_dim = NULL, input_names = NULL, output_names = NULL,  dtype = "float", save_model_as_list = FALSE)

Figure 3: **innsight** utilizes the package **torch**, which builds directly on the C++ library **LibTorch** without a Python dependency.
This object is implemented using the object-oriented R6 class imported from the equally named **R6** package [12]. The only necessary argument is the passed model, which can be either an nn_sequential object from **torch**, a keras_model object from **keras**[13], a neuralnet object from **neuralnet**[13], or a named list in a specific style. The other arguments input_dim, input_names and output_names are optional - except input_dim in combination with **torch** models - and are used for internal validation of the copied model or to assign labels to the input and output nodes used for the visualizations in Step 3 explained in Section 3.3. In addition, the arguments dtype and save_model_as_list specify the calculations' numerical precision and save the entire model as a named list in the instance's field model_as_list, which is created as an intermediate step during the conversion process and are explained in more detail in the next paragraph and Figure 4.

To be as flexible as possible and to interpret almost arbitrary models from any R package, a conversion method is implemented for each of the model classes mentioned above of the three packages **torch**, **keras** and **neuralnet**, summarizing all decisive components and layers of the passed model in an ordered and unified way into a list. Then, a **torch**-based model ConvertedModel (i.e., a subclass of nn_module) is created internally from this list. In addition, the interpretation methods described in Section 2 are pre-implemented for each valid layer type which can be called layer by layer in the following step. Since the creation of the converted model is consequently independent of the class of the given model, the conversion call can be bypassed by directly passing the desired model as a list. Hence, custom wrappers for other packages' models can be written, allowing an interpretation of models not being created by the packages **torch**, **keras** or **neuralnet**. An overview of the individual steps that are performed internally when initializing a new instance of the Converter class is summarized in Figure 4. In addition to the fields shown in Figure 4, there are also fields containing the labels ($input_names, $output_names) and shapes ($input_dim, $output_dim) of the input and output layers in a unified list structure. What kind of list structure is required for a model passed as a list, which layers are generally accepted and even more is explained in detail in the vignette "In-depth explanation" (see vignette("detailed_overview", package = "innsight") or the online documentation available at [https://bips-hb.github.io/innsight/articles/detailed_overview.html](https://bips-hb.github.io/innsight/articles/detailed_overview.html)) and is only referred to at this point.

### Step 2 - Apply selected method

As previously mentioned, the **innsight** package provides the most popular feature attribution techniques in a unified framework. Besides the individual method-specific variations, the

Figure 4: Internal conversion process executed during initialization of an Converter object.
overall structure of each method is the same. Internally, this unification is achieved by the R6 super class InterpretingMethod, from which all methods intended for users inherit and only add method-specific arguments to those of the super class. The rudimentary call of initializing a new method object looks like this:

InterpretingMethod$new(converter, data,  channels_first = TRUE, output_idx = NULL, ignore_last_act = TRUE,  verbose = interactive(), dtype = "float")

The key arguments for every method are the converter object from the first step (see Sec. 3.1), containing the torch-converted model, and the data to be interpreted. The data can be passed in any format as long as the R base method as.array() can convert it into an array and it matches the expected input dimension of the model. In addition, it is common for image or signal data to place either the channel axis directly after the batch axis or at the last position. However, this placement can generally not be extracted unambiguously from the data, which is why the channels_first argument specifies where the channel axis is located, allowing the use of both formats in **innsight**. The remaining arguments output_idx, ignore_last_act, verbose and dtype set which output nodes or classes are to be explained, whether the last activation function is ignored, whether a progress bar is displayed, or change the numerical precision for the calculations.

The feature attribution techniques designed for the package user's regular use cases and applications are inheritors of the super class InterpretingMethod and extend it by method-specific arguments. How the methods from Section 2 are realized, is summarized below:

* The methods _Gradient_ and _Gradient\(\times\)Input_ are implemented as the R6 class Gradient, which has times_input as the only additional argument apart from the inherited ones. This argument switches between the usual gradients (times_input=FALSE) and the gradients multiplied by the corresponding inputs (times_input=TRUE). With **innsight*
* they are applied with the following R code:
* Gradient$new(converter, data, times_input = FALSE,...)
* Gradient$new(converter, data, times_input = TRUE,...)
* Similarly, the methods _SmoothGrad_ and _SmoothGrad\(\times\)Input_ are realized in the R6 class SmoothGrad containing the arguments n for the number of perturbations and noise_level for the noise scale in addition to the times_input argument. The call in **innsight*
* is as follows:
* SmoothGrad$new(converter, data, times_input = FALSE,  n = 50, noise_level = 0.1,...)
* SmoothGrad$new(converter, data, times_input = TRUE,  n = 50, noise_level = 0.1,...)
* The _LRP_ method, including the simple rule ("simple"), \(\varepsilon\)-rule ("epsilon"), \(\alpha\)-\(\beta\)-rule ("alpha_beta"), and a composition of these rules, is implemented in the R6 class LRP. The rule and its corresponding parameter (if available) are set with the arguments rule_name and rule_param. For both arguments, named lists can also be passed to assign a rule or parameter to each layer type separately. Since many zeros are produced in a maximum pooling layer during the backward pass due to the selection of the maximum value in the pooling kernel, the argument winner_takes_all can be used to treat a maximum as an average pooling layer in the backward pass instead. The overall call in **innsight** is the following: # LRP (with defaults) LRP$new(converter, data, rule_name = "simple", rule_param = NULL, winner_takes_all = TRUE,...)

 # LRP with average pooling in the backward pass and # rule "alpha_beta" with alpha = 1 LRP$new(converter, data, rule_name = "alpha_beta", rule_param = 1, winner_takes_all = FALSE,...)

* Analogously, the method _DeepLift_ is realized in the R6 class DeepLift including the argument rule_name for selecting the _Rescale_ ("rescale") or _RevealCancel_ ("reveal_cancel") rule for non-linearities. The reference value is set with x_ref de-faulting to a baseline of zeros. DeepLift can also run into problems in maximum pooling layers since the maximum values in the pooling kernel from the normal and reference input generally do not coincide. Hence, with the winner_takes_all argument, this layer type can be treated as an average pooling layer in a backward pass. With **innsight** DeepLift is applied with the following code: # DeepLift (with defaults) DeepLift$new(converter, data, rule_name = "rescale", x_ref = NULL, winner_takes_all = TRUE,...)

 # DeepLift with average pooling in the backward pass and # rule "reveal_cancel" DeepLift$new(converter, data, rule_name = "reveal_cancel", winner_takes_all = FALSE,...)

* The last method provided by **innsight** is the _connection weights (CW)_ method realized in the R6 class ConnectionWeights. The argument times_input specifies whether the global result of the CW method is calculated or whether it is additionally multiplied by the inputs to obtain local instance-wise explanations. A notable aspect, in this case, is that the data argument is not needed for the global variant but it is required for the local one. The call in **innsight** is as follows: # Connection weights (global) ConnectionWeights$new(converter, times_input = FALSE,...)

 # Connection weights (local) ConnectionWeights$new(converter, data, times_input = TRUE,...)
### Step 3 - Get and visualize the results

After creating an object of a selected method, the third step is to extract the results and, if required, present them in a descriptive and visual way. For this purpose, the **innsight** package provides three generic methods get_result(), plot() and boxplot() that either return the results as an R object (such as an array, torch_tensor or data.frame) or create visualizations for individual instances or aggregated results over the whole passed dataset. All three generic functions call the respective class methods in the InterpretingMethod super class, which are inherited by all the interpreting methods from the second step by design, e.g., for a method's result method, plot(method,...) is the same as method$plot(...).

#### Generic function get_result()

The function get_result() can be used to obtain the results in various forms, whatever is favored according to the user's subsequent workflow or application. This method has only the argument type (besides the method object) which determines the representation of the returned results. By default (type="array"), the result is returned as an R base array, including the input and output names in the corresponding dimensions specified in the first step in the converter object (see Sec. 3.1). The shape of the array is composed of the input shape including the batch size and the number of computed output nodes, i.e., for a tabular input with ten instances and four input variables, the shape is 10\(\times\)4\(\times\)3 if the method was applied to three output nodes in Step 2. In the same way, type="torch_tensor" returns a torch_tensor object having the same shape as the array, but without dimension labels. However, both variants can also return a list or list of lists with the related results as an array or torch_tensor for models with multiple input or output layers. The third and last format of the results is an R base data.frame obtained with type="data.frame".

Included are columns

Figure 5: Overview of the generic get_result() function yielding the results of a method object as either a named array (left), a torch_tensor (middle), or a data.frame (right) depending on the argument type.
for the input instance ("data"), the input and output layer of the model ("model_input" and "model_output"), the input variable ("feature") - possibly also a second one for images ("feature_2") and the channel for signal and image data ("channel") - the output node or class ("output_node"), and the relevance ("value") for the corresponding values. A visual overview of this generic function is given in Figure 5.

#### Generic function plot()

The generic function plot() visualizes individual instances of the result of the method applied before based on the graphic package ggplot2(Wickham, 2016) or the package **plotly**(Sievert et al., 2022) for interactive graphics if the corresponding argument as_plotly is set. The call is executed as follows:

ggplot2-based plot plot(method, data_idx = 1, output_idx = NULL, aggr_channels = "sum", as_plotly = FALSE, same_scale = FALSE)

interactive plotly-based plot plot(method, data_idx = 1, output_idx = NULL, aggr_channels = "sum", as_plotly = TRUE, same_scale = FALSE)

The key arguments for this function are data_idx and output_idx, which specify the indices for the dataset instances and for the desired output nodes or classes whose result is to be visualized. By default, the first data instance and the first computed output node are used. In the argument output_idx, no arbitrary indices can be passed, but only those for which the results were calculated previously in the second step. The further argument aggr_channels can be used to define how the channels are aggregated for image and signal data. Since visualization depends on the data type, it is internally distinguished between tabular/signal data and image data; accordingly, a bar chart or a raster chart is created, as shown in Figure 6 on the left. The relevances in the bars or the pixels are also scaled by color, facilitating a visual comparison; red means positive, blue negative, and white the absence of relevance. Since, in general, the scales vary significantly for the selected output class or data instance, the plots are scaled separately for each value in output_idx and data_idx. When several input layers are to be visualized, the remaining argument same_scale can be used to select whether the individual input layers are also scaled separately in terms of color. This decision depends on the use case illustrated in the melanoma example in Section 4.2. Furthermore, instead of returning objects from the ggplot2 or **plotly** packages, instances of the S4 classes innsight_ggplot2 and innsight_plotly are produced, which are examined in the Section 3.3.4 for advanced visualizations.

#### Generic function boxplot()

Global behavioral patterns and insights into the model's decision-making process can be derived from the results of multiple instances by appropriately summarizing and aggregating them. The generic function boxplot() visualizes these global interpretations over the whole or parts of the given dataset based on the graphics package ggplot2(Wickham, 2016) or the package **plotly**(Sievert et al., 2022) for interactive charts analogous to the previously discussed function plot(). The call for this function is the following:
boxplot(method, output_idx = NULL, data_idx = "all", ref_data_idx = NULL, aggr_channels = "sum", preprocess_FUN = abs, as_plotly = FALSE, same_scale = FALSE,...)

In addition to the identical arguments output_idx, aggr_channels, as_plotly, and same_scale for the plot() function, options for selecting the data points to be aggregated (data_idx), for drawing a reference data point (ref_data_idx) and a pre-process function of the results (preprocess_FUN) are added. Again, the visualization style depends on the type of input data, as in the plot() function; tabular and signal data are displayed as box plots, whereas only a raster plot with the pixel-wise median is rendered for image data due to the high dimensionality. However, if the chart is **plotly**-based, there is a slider to select which quantile to display. Basic examples and an overview of the boxplot() function is given in Figure 6 on the right. Despite the creation of ggplot2 or **plotly** graphs, instances of the S4 class innsight_ggplot2 or innsight_plotly are returned, which are explained in the following section.

### Advanced visualization

The previous two sections have already explained the basic plot() and boxplot() functions. As mentioned, these functions create either an object of the S4 class innsight_ggplot2 (if as_plotly = FALSE) or one of the S4 class innsight_plotly (if as_plotly = TRUE). These functions are intended to generalize the usual ggplot2 or **plotly** objects since, with these packages, the limits of clear visualization possibilities for models with multiple input layers are quickly reached. For example, two charts with different scales for each output node or class need to be generated in a model with images and tabular data as inputs. In this case, a ggplot2-based or **plotly**-based plot is generated for each single input instance and output node and then combined into one large visualization using arrangeGrob() from **gridExtra**(Auguie, 2017) or subplot() from **plotly**, respectively. In contrast, the S4 class innsight_ggplot2 behaves as a wrapper for the ggplot2 object for ordinary models with only one input or

Figure 6: Overview of the visualization tools plot() and boxplot() provided by the **innsight** package depending on the type of input and the argument as_plotly.
output layer. Nevertheless, instances of the innsight_ggplot2 class can be treated and modified as regular ggplot2 objects providing a ggplot2-typical usage by adding, for example, themes, scales or geometric objects; hence the intermediate step via this class is generally not noticeable to the user. For example, the following code is valid:

plot(method) +  ggplot2::theme_bw() +  ggplot2::xlab("My new x label") +  ggplot2::scale_y_continuous(trans = "pseudo_log") +  ggplot2::geom_text(ggplot2::aes(label = signif(value))

Conveniently, all ggplot2 objects are based on the same data.frame, which is also obtained via the get_result() method (see Sec. 3.3.1), i.e., the corresponding column names can be used as variables in the ggplot2 objects, as can be seen in the last line of the code chunk above. For objects of the innsight_plotly class, the entire plot is always created using the plotly::subplot() function. However, this has the consequence that individually assigned modifications are partially overwritten by the grouping, which is why the usual **plotly**-typical adaptations can only be performed after the innsight_plotly object has been printed and returned by the generic print() function for this class, i.e.,

print(plot(method, as_plotly = TRUE)) %>%  plotly::hide_color() %>%  plotly::layout(axis = list(title = "My new x label"))

In addition, generic functions for both S4 classes are implemented, which provide a deeper and more detailed examination of an already created plot through indexing or indexed modification. Section 4.2 demonstrates the application and illustration of some of these generic methods using visualized explanations of a model that takes tabular and image data as inputs. However, for a more detailed description and usage of these classes, please refer to the vignette "In-depth explanation" (see vignette("detailed_overview", package = "innsight") or the online documentation at [https://bips-hb.github.io/innsight/articles/detailed_overview.html](https://bips-hb.github.io/innsight/articles/detailed_overview.html)).

## 4 Illustrations

To exemplify the methods and step-by-step execution of the **innsight** package, a standard dataset with only numerical tabular inputs on a simple model and a more complex dataset with image and tabular data on an extensive non-sequential network are analyzed in the following. The penguin dataset from the **palmerpenguins** package (Horst _et al._, 2020) is used as the simple dataset taking only the numerical variables of bill length and depth, flipper length, and body weight as inputs. The melanoma dataset (Rotemberg _et al._, 2020) of the Kaggle competition1 is taken as the second dataset, which classifies the malignancy or benignity of the skin cell based on images of skin lesions and moles, and patient-level contextual information.
### Example 1: Penguin dataset

In the first example, the penguin dataset provided by the **palmerpenguins** package (Horst _et al._, 2020) is used and a neural network consisting of a dense layer is trained using the **neuralnet** package (Gunther and Fritsch, 2010). Before the **innsight** package can be used, the dataset must be processed and the neural network must be trained on the modified dataset. As a first pre-processing step, only the variables with the species, bill length and depth, flipper length, and body weight are selected, cleaned of missing values, and numerical variables are normalized:

R> library("palmerpenguins") R> data <- na.omit(penguins[, c(1, 3, 4, 5, 6)]) R> data[, 2:5] <- scale(data[, 2:5]) Next, the dataset is divided into training data and test data at a ratio of 75% to 25%:

R> train_idx <- sample.int(nrow(data), as.integer(nrow(data) * 0.75)) R> train_data <- data[train_idx, ] R> test_data <- data[-train_idx, -1] In the second pre-processing step, a network with 128 units in a single hidden layer and the logistic function as activation is fitted on the training data train_data:

R> library("neuralnet") R> model <- neuralnet(species -., + data = train_data, hidden = 128, act.fct = "logistic", + err.fct = "ce", linear.output = FALSE) Now, we follow the three steps that provide and visualize an explanation of the model model on the test data test_data, which were described in detail in Section 3. As a reminder, the first step uses the R6 class Converter to convert the given model to a **torch**-based model with the pre-implemented methods in each layer:

R> library("innsight") R> conv <- Converter$new(model) Then, in the second step, the desired method is selected and applied to the test data test_data via the corresponding R6 class. In this example, the LRP method is used with the \(\alpha\)-\(\beta\)-rule with \(\alpha=2\):

R> lrp <- LRP$new(conv, test_data, rule_name = "alpha_beta", rule_param = 2) In the last step, the results are visualized in two ways: Using the plot() function, the relevances of one instance of the species Adelie (data index 1) and one of the species Gentoo (data index 76) are displayed for both corresponding classes (output nodes 1 and 3). Secondly, the results for the two classes, Adelie and Gentoo, are aggregated over the entire test data and box plots are generated using the boxplot() function without pre-processing and including the first data point as a reference. As mentioned in Section 3, these two variants can be treated and modified like ordinary ggplot2 objects, e.g., adding themes or rotating the x-axis labels. Both visualizations are executed by the following code and can be viewed in Figure 7:
R> library("ggplot2") R> plot(lrp, data_idx = c(1, 76), output_idx = c(1, 3)) + + theme_bw() + + theme(axis.text.x = element_text(angle = 45, vjust = 0.6)) R> boxplot(lrp, output_idx = c(1, 3), preprocess_FUN = identity, + ref_data_idx = 1) + + theme_bw() + + theme(axis.text.x = element_text(angle = 45, vjust = 0.6)) In Figure 6(a) it can be seen that the bill length for the chosen penguin of the Adelie class (index 1 in the dataset test_data) is highly relevant - based on the trained model - for this particular class. However, at the same time, this feature also argues against the Gentoo class due to its strong negative relevance. For the Gentoo penguin, the bottom row in Figure 6(a) reveals that the bill depth is decisively in favor of the Gentoo class and concurrently against the Adelie species. Besides these instance-wise explanations, the boxplot() function provides aggregate insights across the entire test data test_data, summarized in Figure 6(b). The box plots show that the bill length has high positive relevance for the Adelie class and consequently strongly influences it. Simultaneously, however, it also negatively affects the Gentoo class in general. It further emerges that the bill depth and flipper length are crucial features for the Gentoo class.

Figure 7: Generated visualizations of LRP results with the \(\alpha\)-\(\beta\)-rule (\(\alpha=2\)) on the penguin dataset. Sub-figure (a) shows the individual results from data points 1 and 76 from the test data test_data for the Adelie and Gentoo classes. In contrast, the summarized results as box plots across the whole test data for the same two classes can be found in (b), including the individual result of the first data point with the red line.
### Example 2: Melanoma dataset

The second example examines the melanoma dataset (Rotemberg _et al._, 2020) from the Kaggle challenge2 in 2020, issued by the society of imaging informatics in medicine (SIIM) and based on the international skin imaging collaboration (ISIC) archive, the most extensive publicly available collection of quality-controlled dermoscopic images of skin lesions. This dataset consists of 33 126 labeled images with associated patient-level contextual information, such as the age, gender, and image location of the skin lesion or mole.

Footnote 2: See the following link for the official dataset description [https://www.kaggle.com/competitions/sitim-isic-melanoma-classification/overview/description](https://www.kaggle.com/competitions/sitim-isic-melanoma-classification/overview/description)

Due to the complexity and high dimensionality of the data, training a neural network is not straightforward and overall not the main focus of this paper, thus reference is made to the GitHub repository for reproduction ([https://github.com/bips-hb/JSS_innsight/](https://github.com/bips-hb/JSS_innsight/)), and only the most notable points are summarized in the following: The tabular input part's numerical and one-hot encoded categorical variables are fed into a sequential model of dense layers. On the other hand, we use an architecture based on the established residual layers (He _et al._, 2016) considering skip connections between convolutional layers for the image data. Afterward, the two outputs of the respective input parts are merged by concatenation and finally flow in a sequential model with only dense layers to obtain a prediction probability for the skin lesion status. The coarse structure is summarized in Figure 8, where additional dropout layers are used between dense layers. Furthermore, the numerical variable age and the one-hot encoded variables gender and location yield ten features as inputs for the tabular model, and the images are resized to 224\(\times\)224\(\times\)3 for the image model. This model architecture is trained on the melanoma dataset with a validation split of 20% and a batch size of 256 instances using the **Keras** library (Chollet _et al._, 2015) with stochastic gradient descent (SGD) as the optimizer and class-weighted binary cross-entropy as the loss function. However, the best model is selected based on the highest value of the area under the ROC curve (AUC) on the validation data. This metric is chosen because the dataset is highly imbalanced with only 584 of the 33 126 images containing a malignant skin lesion. Since the model is trained from scratch and the image model has significantly more parameters than the tabular one, the training starts with 300 warm-up epochs on the image model using the image data only. Then, the image model is joined with the tabular and the dense output model. Afterward, training continues on the image and tabular data, saving the model with the highest value of the AUC metric on the validation data. In addition, the initial learning rate of 0.01 is reduced by a factor of 0.1 after 20 epochs without a validation AUC improvement, and training is terminated after 40 unimproved epochs. With this approach, a AUC value of 87.71% and an accuracy of 84.19% on the validation data are achieved and the model to be interpreted is selected.

Figure 8: Model architecture for the melanoma dataset.
Based on this model, the obtained predictions can now be explained using the 3-step approach of **innsight**: In the first step, the trained model is loaded and converted to a **torch**-based model with the pre-implemented feature attribution methods using the Converter class. However, since **keras** models do not include names of the input variables and output nodes, these can be passed along when initializing the converter to preserve meaningful labels of the input and output variables in the visualizations. Thus, the first step is executed by the following R code:

R> library("keras") R> library("innsight") R> model <- load_model_tf("path/to/model") R> input_names <- list( + list(paste0("C", 1:3), paste0("H", 1:224), paste0("W", 1:224)), + list(c("Sex: Male", "Sex: Female", "Age", + "Loc: Head/neck", "Loc: Torso", "Loc: Upper extrm.", + "Loc: Lower extrem.", "Loc: Palms/soles", "Loc: Oral/genital", + "Loc: Missing"))) R> output_name <- c("Probability of malignant lesion") R> converter <- Converter$new(model, input_names = input_names, + output_names = output_name)

Next, the LRP method with composite rules is applied, which selects the propagation rule depending on the layer type. For convolutional layers, the \(\alpha\)-\(\beta\)-rule with \(\alpha=1.5\) is used to favor the positive over the negative relevances. In addition, the \(\varepsilon\)-rule with \(\varepsilon=0.01\) is performed on all dense layers and the simple rule - the rule used by default - on average pooling layers. This second step is performed with **innsight** as follows:

R> rule_name <- list(Conv2D_Layer = "alpha_beta", Dense_Layer = "epsilon") R> rule_param <- list(Conv2D_Layer = 1.5, Dense_Layer = 0.01) R> res <- LRP$new(converter, input, channels_first = FALSE, + rule_name = rule_name, rule_param = rule_param)

For the sake of simplicity, the loading of the input data input is omitted in the above code snippet and can be found in the reproduction material together with the whole example. In addition, the channels axis of the images is located at the last position, which is why the argument channels_first must be set to FALSE. The results can be visualized using the implemented plot() function for all interpretability methods. By default, the results are scaled using colors (red for positive and blue for negative relevances) for each instance, each considered output node and each input layer individually. This behavior is especially appropriate for models with multiple input layers consisting of images mixed with tabular data. Because even if the relevances are the same at the end of the tabular and image model before merging, they are further propagated to only ten input variables for the tabular and 224\(\times\)224\(\times\)3 variables for the image model leading to potentially different relevance scales. The following code produces the plot object based on the S4 class innsight_ggplot2 for the first three dataset instances, which can be treated as a ggplot2 object:

R> library("ggplot2") R> p <- plot(res, data_idx = c(2, 3, 1)) + theme_bw()
The order of the indices of the data instances also specifies the ordering of the appearances in the plot, i.e., in the above code, the second instance is visualized first, then the third, and finally, the first instance. Since this model has no standard architecture and the visualization is more extensive, the suggested packages **gridExtra**(**Aguuie 2017**) and **gtable**(**Wickham and Pedersen 2023**) are required. However, each individual plot in the object p can now be modified individually based on the ggplot2 syntax. The indexing works as the objects are plotted in a matrix-wise fashion provided by the facet rows and columns. It is pointed out that each plot object is based on the same dataset, which is also created by the method get_result(type = "data.frame"), i.e., the same column names can be used within the ggplot2 syntax. In the following code snippet, the facet and the tabular x-axis labels are changed manually and the plot is visualized, which can be found in Figure 9:

R> p[1, 1] <- p[1, 1, restyle = FALSE] + + facet_grid(cols = vars(model_input), + labelller = as_labelller(c(Input_1 = "Image input"))) R> p[1, 2] <- p[1, 2, restyle = FALSE] + + facet_grid(rows = vars(data), cols = vars(model_input), + labelller = as_labelller(c(data_2 = "ISIC_6535558 (87.82%)", + Input_2 = "Tabular input"))) R> p[2:3, 2] <- p[2:3, 2, restyle = FALSE] + + facet_grid(rows = vars(data), + labelller = as_labelller(c(data_3 = "ISIC_7291021 (0.05%)", + data_1 = "ISIC_0946787 (47.47%)"))) + + theme(axis.text.x = element_text(angle = 45, vjust = 0.6)) R> plot(p, heights = c(0.31, 0.31, 0.38)) The argument restyle is set when indexing the innsight_ggplot2 object, ensuring that the subplots are extracted in the same way as they are displayed in the whole plot. Otherwise, the entire plot's corresponding facet stripes and axis labels are transferred to the selection. In addition, the arguments in the generic function plot() for innsight_ggplot2 objects are forwarded to the function arrangeGrob() when the plot is finally rendered. This feature allows adjusting the relative heights and widths, demonstrated in the last line of code, to slightly compensate for the increased vertical space of the rotated axis labels.

The three instances in Figure 9 describe different explanatory approaches to the trained model's predictions: The top image in Figure 9a of a malignant lesion was recorded on the torso of a 65-year-old female patient. In the associated interpretation generated by **innsight** (top row in Fig. 9b), it can be observed that, on the one hand, the model identifies the lesioned skin area. On the other hand, the darker and patchy pigmentation and the ragged borders positively influenced the prediction of \(87,82\%\) for melanoma. This observation is also consistent with the official ABCD checklist for melanoma (**Friedman et al.** 1985**), which states that asymmetry, irregular borders, varying color, and large diameters are indicative of a malignant skin lesion. However, the patient's age also positively affected the prediction, as evident from the tabular patient-level information explanation in Figure 9b. A complimentary picture results from the middle image in Figure 9a, showing a benign mole located on the lower extremities of a 40-year-old man. The model saw melanoma for only 0.05% and explains its decision with the symmetrical shape, uniform color pigmentation, and lack of notched borders. In addition, the age of 40 also has a slightly negative influence on the
prediction (middle row in Fig. (b)b). The last instance exemplifies a situation where the model is uncertain whether it is a malignant or benign skin lesion. The truly malignant skin area originates from the lower extremities of a 90-year-old woman (bottom image in Fig. (a)a). Especially the image input explanation in the last row in Figure (b)b shows the model's uncertainty because the mole's upper part looks very regular, arguing for a healthy lesion and consistently highlighted with negative relevance by the model's explanation. In contrast, the lower part contains some notches potentially favoring melanoma, which the model also correctly identified. Furthermore, the high age of the 90-year-old patient has a solid positive relevance to the model's prediction, demonstrating the strong effect of the feature age.

## 5 Validation and runtime

To evaluate the validity and computational performance of **innsight**, the results of the presented feature attribution methods on simulated models and data are compared with the results of the Python implementations **zennit**(Anders _et al._, 2021), **investigate**(Alber _et al._, 2019), **captum**(Kokhlikyan _et al._, 2020) and **deeplift**(Shrikumar _et al._, 2017). The packages **deeplift** and **investigate** are based on the high-level machine learning library **Keras**(Chollet _et al._, 2015) and utilize **TensorFlow**(Abadi _et al._, 2015) as the backend for all calculations. In addition, both packages initially create a replication of the passed model with the interpretations methods pre-implemented in the individual layers similar to **innsight**. In contrast, the packages **zennit** and **captum** use **PyTorch**(Paszke _et al._, 2019) and run without a conversion

Figure 9: The image part of the instance of the melanoma dataset to be explained and the associated visualization generated by **innsight**. Figure (a) shows a (top) malignant lesion image of a 65-year-old female, (middle) benign lesion of a 40-year-old male and (bottom) malignant lesion of a 90-year-old female patient. Figure (b) displays the LRP explanation of the patients from (a) created with the plot() function and subsequent minor modifications such as facet and x-axis labels.
step since hooks are used to modify the automated backward pass according to the applied method on the fly. However, this only enables the application of methods that can be considered independent of the preceding and following layers, which complicates, for example, an implementation of DeepLift with the RevealCancel rule. Furthermore, not every package supports all methods. The gradient-based methods Gradient and Gradient\(\times\)Input are provided by all packages. In contrast, DeepLift with Rescale rule is only implemented in **deeplift** and **captum**, and the RevealCancel rule only in **deeplift**. The LRP methods are available in all packages except for **deeplift**. However, **captum** does not accept the \(\alpha\)-\(\beta\)-rule.

### Validity comparison

For the validation, shallow untrained dense and convolutional models with the most commonly used layer types - such as 2D convolution, 2D maximum/average pooling and dense layers - and normally distributed input data are generated. More specifically, 32 different architectures are considered, using ReLU and hyperbolic tangent to include both constrained and unconstrained activation functions, with and without bias vectors, with varying pooling layers and a different number of output nodes. From each of these architectures, 50 randomly initiated models are created, resulting in 1600 distinct models, and evaluated on normally distributed datasets with 32 input instances each. The experimental details can be found in the Appendix A.1. Moreover, all figures and results are reproducible using the code in the reproduction material or on GitHub ([https://github.com/bips-hb/JSS_innsight](https://github.com/bips-hb/JSS_innsight)).

As a measure of quality, the mean absolute error (MAE) between the result of **innsight** and the corresponding reference implementation over all input variables and output nodes is considered. Consequently, for each combination of method, model, input instance, and output node, a value for this quality measure is derived, leading to box plots for visualizing the differences. In addition to the box plots, the acceptable error range of up to \(10^{-6}\) is highlighted in light gray to distinguish numerical tolerated differences caused by calculations of single-precision floating point numbers according to the IEEE 754 standard (IEEE 2019)

Figure 10: Comparison of feature attribution methods’ results of **innsight** and the reference implementations **captum**, **zennit**, **investigate**, and **deeplift** regarding the mean absolute error as box plots over different model architectures and repetitions. It shows the results separated into (a) gradient-based methods, (b) DeepLift and (c) LRP. The shaded gray area indicates the error tolerance of \(10^{-6}\).
from abnormal discrepancies. The results are summarized in Figure 10. For the gradient-based methods Gradient and Gradient\(\times\)Input, the method's results coincide precisely for **PyTorch**-based packages, and they differ at most by \(10^{-7}\) for the **Keras**-based packages (see Fig. 9(a)). The main reason for this discrepancy is that **innsight**, **zennit**, and **captum** compute gradients via **PyTorch** or **PyTorch**'s C++ library **LibTorch**, and **investigate** and **deeplift** use **TensorFlow** for the calculations. A similar picture results for the DeepLift method with the Rescale and RevealCancel rules but with few outliers with a maximum error of up to \(10^{-3}\) for the Rescale rule (see Fig. 9(b)). However, all outliers with an error exceeding \(10^{-6}\) originate from models with the hyperbolic tangent as activation and can thus be explained by numerical inaccuracies due to the saturated activation. In addition, minor discrepancies are probably caused by different treatments of vanishing denominators in the multipliers or numerical uncertainties between the backends **PyTorch**/**LibTorch** and **TensorFlow** in general. For the LRP methods, a few adjustments are needed for the **investigate** and **captum** packages since they use only the simple rule for average pooling layers, which is modified in **innsight** using composite rules. Apart from that, the results from **innsight** compared to **captum** or **zennit** for the simple, \(\varepsilon\)-rule and \(\alpha\)-\(\beta\)-rule differ negligibly and are far below the maximally tolerated error of \(10^{-6}\) (see Fig. 9(c)). For the simple and \(\varepsilon\)-rule, **investigate** is consistent with **innsight** except for a few deviations. Again, 96.8% of the cases with errors exceeding \(10^{-6}\) are caused by a saturated hyperbolic tangent activation, lower errors on different stabilizers for the denominators in the relevance messages, and general numerical inaccuracies between their backends. However, significant discrepancies can be observed using the \(\alpha\)-\(\beta\)-rule, which only occur in models with a bias vector (see Fig. 9(c) bottom). The reason for this is a different interpretation of the positive or negative part of the bias vector, which is discussed in more detail in the Appendix B.

### Runtime comparison

In addition to comparing whether **innsight**'s results are consistent with the reference implementations, a runtime comparison is also conducted concerning the number of output nodes, hidden units or filters, hidden layers, batch size, and for images, the size of the input images. It must be noted again that the packages based on **Keras** and **innsight** first convert the passed model, and the **PyTorch**-based packages use hooks to overwrite the automated backward pass while executing, making them considerably faster. Therefore, in the results, only the execution time excluding the conversion step - as far as possible - is presented and not the total time. For comparisons of the total time needed to calculate an explanation or conversion time only, see Appendix A.3. However, the **investigate** package has a special characteristic in this regard since the entire conversion

Figure 11: Model architectures for the run-time comparison. The hyperparameters for the number of outputs (C), number of hidden units or filter size (U), number of hidden layers (L), batch size (B) and the height/width of the input images (W) were varied in each case.
process and the construction of the underlying graph only happens during the analysis of the first batch of input data3. For this reason, conversion times are almost hardly present in the results. Since this simulation assumes that an interpretation method is being applied for the first time to a model and only to a single batch of input instances, the results of **investigate** are slightly biased and would be notably quicker if the same model is employed with more input batches. Analogously to the comparison from Section 5.1, untrained dense and convolutional neural networks with the architectures shown in Figure 11, and normally distributed input data are used for the time comparisons. Depending on the type of time comparison, the hyperparameters for the number of output nodes (C), number of hidden units or filter size (U), number of hidden layers (L), batch size (B), and the size of the input images (W) are varied. The hyperparameters not considered in the respective comparison remain unchanged and take default values, i.e., one output node, 10 hidden units for the tabular model and 5 filters for the image model, two layers in total, a batch size of 16 and an input image size of 10. In addition, 20 replicates of each architecture are created to compensate for potential numerical fluctuations. For a more detailed simulation description or analysis of the results including the conversion time, please refer to Appendix A.2 and A.3, and for a reproduction of the results, see the reproduction material or the code in the GitHub repository at [https://github.com/bips-hb/JSS_innsight/](https://github.com/bips-hb/JSS_innsight/).

Footnote 3: See the GitHub issue #50 and #129 for the **investigate** package

In general, comparing the runtimes of the different packages reveals that **innsight** is faster than **ininvestigate** and **deeplift** (which are based on **Keras**), but slower than **captum** and **zennit** (which are based on **PyTorch**). This overall trend is particularly evident when adding more layers to the models since **innsight** is 10-15 times slower than **captum** and **zennit**, but still for the same order of magnitude faster than **investigate** and **deeplift**. One notable exception is that for deeper models using LRP with \(\alpha\)-\(\beta\)-rule, **innsight** is 20 times faster than **ininvestigate**, and only slightly faster than **deeplift** using the DeepLift method (see Fig. 12a). However, there are both positive and negative deviations from this trend, especially when varying the number of output nodes. In these cases, **innsight** stands out with the LRP and DeepLift methods, and performs similarly or even faster than the **PyTorch**-based packages (see Fig.12b). The primary reason for this, however, is that the results for several output nodes can be calculated at once in **innsight**. In contrast, all other implementations only allow the calculation for single nodes and thus the method's results are computed by iterative execution. Although **innsight** and **investigate** become more comparable to the **PyTorch** packages with larger inputs and more filters and hidden units, a weakness of **innsight** can be stated: Even if the results on tabular data almost correspond to those of the **PyTorch**-based packages, **innsight** is slower than all reference implementations for larger image sizes in the DeepLift method and more filters in LRP and DeepLift (see Fig. 16 and 18 in the appendix). Even though the packages in this section are compared regarding runtime and some packages' weaknesses are revealed, all considered implementations provide an explanation within a reasonable runtime of a few seconds, even for deep neural networks with several large images as inputs.

## 6 Summary and discussion

In summary, we have presented **innsight**, an R package that provides the most well-known feature attribution methods for interpreting neural network predictions. After a detailed
introduction of the implemented feature attribution methods, the internal structure utilizing **torch**'s fast array calculations and the R6 class Converter demonstrated how the deep-learning-model-agnostic approach was implemented to enable the analysis of models from any R package in an efficient way. This flexibility is complemented by a unified 3-step approach from model to plotted results, including multiple visualization tools based on ggplot2 or **plotly** for interactive plots. The step-wise procedure was illustrated using a model on the tabular penguin dataset and a deep neural network on the melanoma dataset consisting of structured patient-level information and images. Furthermore, the results of the simulation study show that **innsight** returns identical feature-wise explanations to the reference implementations **captum**, **zennit**, **investigate**, and **deeplift** in Python, except for negligible numerical inaccuracies. In terms of runtime, the package shows that it is generally faster than the **Keras**-based packages but slower than **captum** and **zennit**, and it only suffers during convolution if the image height/width or the number of filters is large. Apart from that, the package also has some limitations that could be improved in the future. For example, only converting sequential models (i.e., nn_sequential) from the **torch** package is possible because no structured network graph can be extracted from an arbitrary nn_module. Nevertheless, passing a model

Figure 12: Package’s average evaluation time in seconds over 20 repetitions for applying different feature attribution methods on models with (a) a varying number of hidden layers and (b) a varying number of output nodes (only image data).
as a list allows the user to do the conversion step on their own in such cases. Furthermore, an activation function is assigned to a linear or convolutional layer only if it is defined in the layer itself or immediately after the layer. This behavior is especially relevant for the RevealCancel rule in the DeepLift method, because **insight** handles separated activations with the Rescale rule, which is the case, for example, with the layer sequence of convolution, batch normalization, and activation. Moreover, even if it is possible in the **torch** package, the **insight** package currently only supports computations on CPUs and not on GPUs.

## Computational details

A 64-bit Linux platform running Ubuntu 20.04 with an AMD Ryzen Threadrippper 3960X (24 cores, 48 threads) CPU including 256 gigabyte RAM and two NVIDIA Titan RTX GPUs was used for all computations. All comparisons and calculations with the reference implementations were performed in a separate session - created by **callr**(Csardi and Chang, 2022) - using only a single CPU thread per job. An exception was the neural network training on the melanoma dataset using a single GPU, which was also the only code executed by Python and not from R 4.3.0 (R Core Team, 2023). Due to the package requirement mismatch, separate environments were created for each of the **Keras**-based and the **PyTorch**-based packages, i.e.,

* **investigate**2.0.2: Using Python 3.8.15 with **Keras**2.10.0 and **TensorFlow**2.10
* **deeplift**0.6.13: Using Python 3.6.15 with **Keras**2.2.4 and **TensorFlow**1.15
* **captum**0.6.0 and **zennit**0.5.0: Using Python 3.8.12 with **PyTorch**1.13.1 (cpu)

The corresponding environments were loaded in R and then the code was executed in Python using **reticulate**1.30(Ushey _et al._, 2023). In addition, the computer was used exclusively for the runtime measurements for the corresponding job and was not distorted by other simultaneous processes.
## Acknowledgments

This project was funded by the German Research Foundation (DFG), Emmy Noether Grant 437611051.
