Package torch

Currently, only models created by torch::nn_sequential are accepted. However, the most popular standard layers and activation functions are available:

• Linear layers: nn_linear()

• Convolutional layers: nn_conv1d(), nn_conv2d() (but only with padding_mode = "zeros" and numerical padding)

• Max-pooling layers: nn_max_pool1d(), nn_max_pool2d() (both only with default arguments for padding = 0, dilation = 1, return_indices = FALSE and ceil_mode = FALSE)

• Average-pooling layers: nn_avg_pool1d(), nn_avg_pool2d() (both only with default arguments for padding = 0, return_indices = FALSE and ceil_mode = FALSE)

• Batch-normalization layers: nn_batch_norm1d(), nn_batch_norm2d()

• Flatten layer: nn_flatten()

• Skipped layers: nn_dropout()

• Activation functions: nn_relu, nn_leaky_relu, nn_softplus, nn_sigmoid, nn_softmax, nn_tanh (open an issue if you need any more)

📝 Notes

In a torch model, the shape of the inputs is not stored; hence it must be specified with the argument input_dim within the initialization of the Converter object.

library(torch)

torch_model <- nn_sequential(
  nn_conv2d(3, 5, c(2, 2), stride = 2, padding = 3),
  nn_relu(),
  nn_avg_pool2d(c(2, 2)),
  nn_flatten(),
  nn_linear(80, 32),
  nn_relu(),
  nn_dropout(),
  nn_linear(32, 2),
  nn_softmax(dim = 2)
)

# For torch models the optional argument `input_dim` becomes a necessary one
converter <- convert(torch_model, input_dim = c(3, 10, 10))

Package keras

Keras models created by keras_model_sequential or keras_model are accepted. Within these functions, the following layers are allowed to be used:

• Input layer: layer_input()

• Linear layers: layer_dense()

• Convolutional layers: layer_conv_1d(), layer_conv_2d()

• Pooling layers: layer_max_pooling_1d(), layer_max_pooling_2d(), layer_average_pooling_1d(), layer_average_pooling_2d(), layer_global_average_pooling_1d(), layer_global_average_pooling_2d(), layer_global_max_pooling_1d(), layer_global_max_pooling_2d()

• Batch-normalization layer: layer_batch_normalization()

• Flatten layer: layer_flatten()

• Merging layers: layer_add(), layer_concatenate() (but it is assumed that the concatenation axis points to the channel axis)

• Padding layers: layer_zero_padding_1d(), layer_zero_padding_2d()

• Skipped layers: layer_dropout()

• Activation functions: The following activation functions are allowed as character argument (activation) in a linear and convolutional layer: "relu", "softplus", "sigmoid", "softmax", "tanh", "linear". But you can also specify the activation function as a standalone layer: layer_activation_relu(), layer_activation_softmax(). But keep in mind that an activation layer may only follow a dense, convolutional or pooling layer. If you miss an activation function, feel free to open an issue on GitHub.

library(keras)

# Create model
keras_model_seq <- keras_model_sequential()
keras_model_seq <- keras_model_seq %>%
  layer_dense(10, input_shape = c(5), activation = "softplus") %>%
  layer_dense(8, use_bias = FALSE, activation = "tanh") %>%
  layer_dropout(0.2) %>%
  layer_dense(4, activation = "softmax")

converter <- convert(keras_model_seq)

library(keras)

input_image <- layer_input(shape = c(10, 10, 3))
input_tab <- layer_input(shape = c(20))

conv_part <- input_image %>%
  layer_conv_2d(5, c(2, 2), activation = "relu", padding = "same") %>%
  layer_average_pooling_2d() %>%
  layer_conv_2d(4, c(2, 2)) %>%
  layer_activation(activation = "softplus") %>%
  layer_flatten()

output <- layer_concatenate(list(conv_part, input_tab)) %>%
  layer_dense(50, activation = "relu") %>%
  layer_dropout(0.3) %>%
  layer_dense(3, activation = "softmax")

keras_model_concat <- keras_model(inputs = list(input_image, input_tab), outputs = output)

converter <- convert(keras_model_concat)

Package neuralnet

Using nets from the package neuralnet is very simple and straightforward, because the package offers much fewer options than torch or keras. The only thing to note is that no custom activation function can be used. However, the package saves the names of the inputs and outputs, which can, of course, be overwritten with the arguments input_names and output_names when creating the converter object.

library(neuralnet)
data(iris)
set.seed(42)

# Create model
neuralnet_model <- neuralnet(Species ~ Petal.Length + Petal.Width, iris,
  linear.output = FALSE
)

# Convert model
converter <- convert(neuralnet_model)
# Show input names
converter$input_names
# Show output names
converter$output_names

Model as named list

Besides models from the packages keras, torch and neuralnet it is also possible to pass a self-defined model in the form of a named list to the Converter class/convert() function. This enables the interpretation of networks from other libraries with all available methods provided by the innsight package.

If you want to create a custom model, your list (e.g., model <- list()) needs at least the keys model$input_dim and model$layers. However, other optional keys that can be used to name the input features and output nodes or to test the model for correctness. In summary:

• input_dim
  The model input dimension excluding the batch dimension in the format “channels first”. If there is only one input layer, it can be specified as a vector. Otherwise, use a list of the shapes of the individual input layers in the correct order.

  Examples

  1. one dense input layer with five features:

  model$input_dim <- c(5)

  2. two input layers consisting of one dense layer of shape 10 and one convolutional layer of shape 10 × 10 × 3:

  input_dim <- list(c(10), c(3, 10, 10)) # channels have to be first!

  
  
• input_nodes
  An integer vector specifying the indices of the input layers from the layers list. If there are multiple input layers, the indices of the corresponding layers must be in the same order as in the input_dim and input_names arguments. If this argument is not set, a warning is printed and it is assumed that the first layer in layers is the only input layer.
  
  

• input_names (optional)
  The names for each input dimension excluding the batch axis in the format “channels first”. If there is only one input layer, it can be specified as a list of character vectors or factors for each input axis. Otherwise, use a list of those for each input layer in the correct order.

  Examples

  1. one dense layer with two features:

  model$input_names <- c("Feature_1", "Feature_2")

  2. one convolutional layer with shape 3 × 3 × 2:

  model$input_names <- list(c("C1", "C2"), c("H1", "H2"), c("W1", "W2"))

  3. two input layers consisting of one dense layer of shape 2 and one convolutional layer of shape 3 × 3 × 2:

  model$input_names <- list(
    list(c("Feature_1", "Feature_2")),
    list(c("C1", "C2"), c("H1", "H2"), c("W1", "W2"))
  )

• output_dim (optional)
  An integer vector or a list of vectors with the model output dimension without the batch dimension analogous to input_dim. This value does not need to be specified and will be calculated otherwise. However, if it is set, the calculated value will be compared with it to avoid errors while converting the model.

• output_nodes
  An integer vector specifying the indices of the output layers from the layers list. If there are multiple output layers, the indices of the corresponding layers must be in the same order as in the output_dim and output_names arguments. If this argument is not set, a warning is printed and it is assumed that the last layer in layers is the only output layer.
  
  

• output_names (optional)
  A list or a list of lists with the names for each output dimension for each output layer analogous to input_names. By default (NULL), the names are generated.

• layers
  (see next subsection)

📝 Notes
The arguments for the input and output names are optional. By default (NULL), they are generated, i.e.,

• the output names are c("Y1", "Y2", "Y3", ... ) for each output layer.

• the input names are

  • c("X1", "X2", "X3", ...) for tabular input layers,
    
  • list(c("C1", "C2", ...), c("L1", "L2", ...)) for 1D input layers and
    
  • list(c("C1", "C2", ...), c("H1", "H2", ...), c("W1", "W2", ...)) for 2D input layers.

# Define dense layer
dense_layer <- list(
  type = "Dense",
  input_layers = 0, # '0' means input layer
  output_layers = 2,
  weight = matrix(rnorm(5 * 2), 2, 5),
  bias = rnorm(2),
  activation_name = "tanh",
  dim_in = 5, # optional
  dim_out = 2 # optional
)

# 1D convolutional layer
conv_1D <- list(
  type = "Conv1D",
  input_layers = 1,
  output_layers = 3,
  weight = array(rnorm(8 * 3 * 2), dim = c(8, 3, 2)),
  bias = rnorm(8),
  padding = c(2, 1),
  activation_name = "tanh",
  dim_in = c(3, 10), # optional
  dim_out = c(8, 9) # optional
)

# 2D convolutional layer
conv_2D <- list(
  type = "Conv2D",
  input_layes = 3,
  output_layers = 5,
  weight = array(rnorm(8 * 3 * 2 * 4), dim = c(8, 3, 2, 4)),
  bias = rnorm(8),
  padding = c(1, 1, 0, 0),
  dilation = c(1, 2),
  activation_name = "relu",
  dim_in = c(3, 10, 10) # optional
)

# 2D average pooling layer
avg_pool2D <- list(
  type = "AveragePooling2D",
  input_layers = 1,
  output_layers = 3,
  kernel_size = c(2, 2)
)

# batch normalization layer
batchnorm <- list(
  type = "BatchNorm",
  input_layers = 1,
  output_layers = 3,
  num_features = 3,
  eps = 1e-4,
  gamma = c(1.1, 0.0, -0.3),
  beta = c(1, -3, -1.4),
  run_mean = c(-1.9, 3, 2.3),
  run_var = c(1, 2.1, 3.5)
)

# flatten layer
flatten <- list(
  type = "Flatten",
  input_layers = 1,
  output_layers = 3,
  start_dim = 2, # optional
  end_dim = 4, # optional
  dim_in = c(3, 10, 10), # optional
  out_dim = c(300) # optional
)

# global MaxPooling layer
global_max_pool2D <- list(
  type = "GlobalPooling",
  input_layers = 1,
  output_layers = 3,
  method = "max",
  dim_in = c(3, 10, 10), # optional
  out_dim = c(3) # optional
)

# global AvgPooling layer
global_avg_pool1D <- list(
  type = "GlobalPooling",
  input_layers = 1,
  output_layers = 3,
  method = "average",
  dim_in = c(3, 10), # optional
  out_dim = c(3) # optional
)

# padding layer
padding <- list(
  type = "Padding",
  input_layers = 1,
  output_layers = 3,
  padding = c(2, 4),
  mode = "constant",
  value = 1,
  dim_in = c(3, 10), # optional
  out_dim = c(3, 16) # optional
)

# concatenation layer
concat <- list(
  type = "Concatenation",
  input_layers = c(1, 3),
  output_layers = 5,
  dim = 2,
  dim_in = list(c(5), c(3)), # optional
  out_dim = c(8) # optional
)

# adding layer
add <- list(
  type = "Add",
  input_layers = c(1, 3),
  output_layers = 5,
  dim_in = list(c(3, 10, 10), c(3, 10, 10)), # optional
  out_dim = c(3, 10, 10) # optional
)

Argument dtype

This argument defines the numerical floating-point number’s precision with which all internal calculations are performed. Accepted are currently 32-bit floating point ("float" the default value) and 64-bit floating point numbers ("double"). All weights, constants and inputs are then converted accordingly into the data format torch_float() or torch_double().

📝 Note
At this point, this decision is especially crucial for exact comparisons, and if the precision is too inaccurate, errors could occur. See the following example:

torch_manual_seed(123)
A <- torch_randn(10, 10)
B <- torch_randn(10, 10)

# result of first row and first column after matrix multiplication
res1 <- torch_mm(A, B)[1, 1]
# calculation by hand
res2 <- sum(A[1, ] * B[, 1])

# difference:
res1 - res2

torch_manual_seed(123)
A <- torch_randn(10, 10, dtype = torch_double())
B <- torch_randn(10, 10, dtype = torch_double())

# result of first row and first column after matrix multiplication
res1 <- torch_mm(A, B)[1, 1]
# calculation by hand
res2 <- sum(A[1, ] * B[, 1])

# difference:
res1 - res2

Argument save_model_as_list

As already described in the introduction vignette, a given model is first converted to a list and then the torch model is created from it. By default, however, this list is not stored in the Converter object, since this requires a lot of memory for large models and is otherwise not used further. With the logical argument save_model_as_list, this list can be stored in the field model_as_list for further investigations. For example, this list can again be used as a model for a new Converter instance.

Argument converter

The Converter object from the first step is one of the crucial elements for the application of a selected method because it converted the original model into a torch structure necessary for innsight in which the methods are pre-implemented in each layer.

Argument data

In addition to the converter object, the input data is also essential as it will be analyzed and explained using the methods provided at the end. Accepted are data as:

• Base R data types like matrix, array, data.frame or other array-like formats of size (batch_size, input_dim). These formats can be used mainly when the model has only one input layer. Internally, the data is converted to an array using the as.array function and stored as a torch_tensor in the given dtype afterward.

• torch_tensor: The converting process described in the last point can also be skipped by directly passing the data as a torch_tensor of size (batch_size, input_dim).

• list: You can also pass a list with the corresponding input data according to the upper points for each input layer.

📝 Note
The argument data is a necessary argument only for the local interpretation methods. Otherwise, it is unnecessary, e.g., the global variant of the Connection Weights method can be used without data.

Argument channels_first

This argument tells the package where the channel axis for images and signals is located in the input data. Internally, all calculations are performed with the channels in the second position after the batch dimension (“channels first”), e.g., c(10,3,32,32) for a batch of ten images with three channels and a height and width of 32 pixels. Thus input data in the format “channels last” (i.e., c(10,32,32,3) for the previous example) must be transformed accordingly. If the given data has no channel axis, use the default value TRUE.

Argument output_idx

These indices specify the model’s output nodes for which the method is to be applied. For the sake of models with multiple output layers, the method object gives the following possibilities to select the indices of the output nodes in the individual output layers:

• A vector of indices: If the model has only one output layer, the values correspond to the indices of the output nodes, e.g., c(1,3,4) for the first, third and fourth output node. If there are multiple output layers, the indices of the output nodes from the first output layer are considered.

• A list of index vectors: If the method is to be applied to output nodes from different layers, a list can be passed that specifies the desired indices of the output nodes for each output layer. Unwanted output layers have the entry NULL instead of a vector of indices, e.g., list(NULL, c(1,3)) for the first and third output node in the second output layer.

• NULL (default): The method is applied to all output nodes in the first output layer but is limited to the first ten as the calculations become more computationally expensive for more output nodes.

Argument output_label

These values specify the output nodes for which the method is to be applied and can be used as an alternative to the argument output_idx. Only values that were previously passed with the argument output_names in the converter can be used. In order to allow models with multiple output layers, there are the following possibilities to select the names of the output nodes in the individual output layers:

• A character vector or factor of labels: If the model has only one output layer, the values correspond to the labels of the output nodes named in the passed Converter object, e.g., c("a", "c", "d") for the first, third and fourth output node if the output names are c("a", "b", "c", "d"). If there are multiple output layers, the names of the output nodes from the first output layer are considered.

• A list of charactor/factor vectors of labels: If the method is to be applied to output nodes from different layers, a list can be passed that specifies the desired labels of the output nodes for each output layer. Unwanted output layers have the entry NULL instead of a vector of labels, e.g., list(NULL, c("a", "c")) for the first and third output node in the second output layer.

• NULL (default): The method is applied to all output nodes in the first output layer but is limited to the first ten as the calculations become more computationally expensive for more output nodes.

Argument ignore_last_act

Set this logical value to include the last activation function for each output layer, or not (default: TRUE). In practice, the last activation (especially for softmax activation) is often omitted.

Argument dtype

This argument defines the numerical precision with which all internal calculations are performed. Accepted are currently 32-bit floating point ("float" the default value) and 64-bit floating point numbers ("double"). All weights, constants and inputs are then converted accordingly into the data format torch_float() or torch_double(). See the argument dtype in the Converter object for more details.

Vanilla Gradient

One of the first and most intuitive methods for interpreting neural networks is the Gradients method introduced by Simonyan et al. (2013), 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 prediction-sensitive variables, i.e., 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 x ∈ ℝⁿ, the model f : ℝⁿ → ℝ^C and the output y_c = f(x)_c of class c: $$ \text{Gradient}(x)_i^c = \frac{\partial\ f(x)_c}{\partial\ x_i} = \frac{\partial\ y_c}{\partial\ x_i} $$

As described in the introduction of this section, the corresponding innsight-method Gradient inherits from the super class InterpretingMethod, meaning that we need to change the term Method to Gradient. Alternatively, an object of the class Gradient can also be created using the mentioned helper function run_grad(), which does not require prior knowledge of R6 objects. The only model-specific argument is times_input, which can be used to switch between the two methods Gradient (default FALSE) and Gradient×Input (TRUE). For more information on the method Gradient×Input see this subsection.

# R6 class syntax
grad <- Gradient$new(converter, data,
  times_input = FALSE,
  ... # other arguments inherited from 'InterpretingMethod'
)

# Using the helper function 
grad <- run_grad(converter, data,
  times_input = FALSE,
  ... # other arguments inherited from 'InterpretingMethod'
)

data <- matrix(c(0.45), 1, 1)

# Apply method (but don't ignore last activation)
grad <- run_grad(converter, data, ignore_last_act = FALSE)
# get result
get_result(grad)

SmoothGrad

The 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. These have large fluctuations and abruptly change their gradients, making the interpretations of the gradient worse and potentially misleading. Smilkov et al. proposed that instead of calculating only the gradient in x, compute the gradients of randomly perturbed copies of x and determine the mean gradient from that. To use the SmoothGrad method to obtain relevance values for the individual components x_i ∈ ℝ of an instance x ∈ ℝ^𝕟, we first generate K ∈ ℕ realizations of a multivariate Gaussian distribution p = 𝒩(0, σ²) describing the random perturbations, i.e., ε¹, …, ε^K ∼ p. Then the empirical mean of the gradients for variable x_i and output index c can be calculated as follows:

$$ \text{SmoothGrad}(x)_i^c = \frac{1}{K} \sum_{j = 1}^K \frac{\partial\ f(x + \varepsilon^j)_c}{\partial\ x_i + \varepsilon_i^j} \approx \mathbb{E}_{\varepsilon \sim p}\left[ \frac{\partial\ f(x + \varepsilon)_c}{\partial\ x_i + \varepsilon_i^j}\right] $$

As described in the introduction of this section, the innsight method SmoothGrad inherits from the super class InterpretingMethod, meaning that we need to change the term Method to SmoothGrad or use the helper function run_smoothgrad() for initializing an object of class SmoothGrad. In addition, there are the following three model-specific arguments:

• n (default: 50): This integer value specifies how many perturbations will be used to calculate the mean gradient, i.e., the K from the formula above. However, it must be noted that the computational effort increases by a factor of n compared to the Gradient method since the simple Gradient method is used n times instead of once. In return, the accuracy of the estimator increases with a larger n.

• noise_level (default: 0.1): With this argument, the strength of the spread of the Gaussian distribution can be given as a percentage, i.e., noise_level $=\frac{\sigma}{\max(x)-\min(x)}$.

• times_input (default: FALSE): Similar to the Gradient method, this argument can be used to switch between the two methods SmoothGrad (FALSE) and SmoothGrad×Input (TRUE). For more information on the method SmoothGrad×Input see this subsection.

# R6 class syntax
smoothgrad <- SmoothGrad$new(converter, data,
  n = 50,
  noise_level = 0.1,
  times_input = FALSE,
  ... # other arguments inherited from 'InterpretingMethod'
)

# Using the helper function 
smoothgrad <- run_smoothgrad(converter, data,
  n = 50,
  noise_level = 0.1,
  times_input = FALSE,
  ... # other arguments inherited from 'InterpretingMethod'
)

data <- matrix(c(0.6), 1, 1)

# Apply method
smoothgrad <- run_smoothgrad(converter, data,
  noise_level = 0.2,
  n = 50,
  ignore_last_act = FALSE # include the tanh activation
) 
# get result
get_result(smoothgrad)

Gradient×Input and SmoothGrad×Input

The methods Gradient×Input and SmoothGrad×Input are as simple as they sound: the gradients are calculated as in the gradient section and then multiplied by the respective input. They were introduced by Shrikumar et al. (2016) and have a well-grounded mathematical background despite their simple idea. The basic idea is to decompose the output according to its relevance to each input variable, i.e., we get variable-wise additive effects

$$ \tag{1} f(x)_c = \sum_{i = 1}^n R_i. $$

Mathematically, this method is based on the first-order Taylor decomposition. Assuming that the function f is continuously differentiable in x ∈ ℝⁿ, a remainder term ε(f, z, x) : ℝⁿ → ℝ with lim_z → xε(f, z, x) = 0 exists such that

$$ \begin{align} f(z) &= f(x) + \nabla_x\ f(x)(z-x)^T + \varepsilon(f, z, x)\\ &= f(x) + \sum_{i = 1}^n \frac{\partial\ f(x)}{\partial\ x_i} (z_i - x_i) + \varepsilon(f, z, x), \quad z\in \mathbb{R}^n. \tag{2} \end{align} $$ The first-order Taylor formula thus describes a linear approximation of the function f at the point x since only the first derivatives are considered. Consequently, a highly nonlinear function f is well approximated in a small neighborhood around x. For larger distances from x, sufficient small values of the residual term are not guaranteed anymore. The Gradient×Input method now considers the data point x and sets z = 0. In addition, the residual term and the summand f(0)_c are ignored, which then results in the following approximation of f(x)_c in variable-wise relevances

$$ f(x)_c \approx \sum_{i = 1}^n \frac{\partial\ f(x)_c}{\partial\ x_i} \cdot x_i, \quad \text{hence}\\ \text{Gradient$\times$Input}(x)_i^c = \frac{\partial\ f(x)_c}{\partial\ x_i} \cdot x_i. $$

Analogously, this multiplication is also applied to the SmoothGrad method in order to compensate for local fluctuations: $$ \text{SmoothGrad$\times$Input}(x)_i^c = \frac{1}{K} \sum_{j = 1}^K \frac{\partial\ f(x + \varepsilon^j)_c}{\partial\ x_i + \varepsilon_i^j} \cdot (x_i + \varepsilon_i^j),\quad \varepsilon^1, \ldots, \varepsilon^K \sim \mathcal{N}(0,\sigma^2). $$

Both methods are variants of the respective gradient methods Gradient and SmoothGrad and also have the corresponding model-specific arguments and helper functions for the initialization. These variants can be chosen with the argument times_input:

# the "x Input" variant of method "Gradient"
grad_x_input <- Gradient$new(converter, data,
  times_input = TRUE,
  ... # other arguments of method "Gradient"
) 

# the same using the corresponding helper function
grad_x_input <- run_grad(converter, data,
  times_input = TRUE,
  ... # other arguments of method "Gradient"
) 

# the "x Input" variant of method "SmoothGrad"
smoothgrad_x_input <- SmoothGrad$new(converter, data,
  times_input = TRUE,
  ... # other arguments of method "SmoothGrad"
)

# the same using the corresponding helper function
smoothgrad_x_input <- run_smoothgrad(converter, data,
  times_input = TRUE,
  ... # other arguments of method "SmoothGrad"
)

data <- matrix(c(0.49), 1, 1)

# Apply method
grad_x_input <- run_grad(converter, data,
  times_input = TRUE,
  ignore_last_act = FALSE # include the tanh activation
) 
# get result
get_result(grad_x_input)

data <- matrix(c(0.49), 1, 1)

# Apply method
smoothgrad_x_input <- run_smoothgrad(converter, data,
  times_input = TRUE,
  ignore_last_act = FALSE # include the tanh activation
) 
# get result
get_result(smoothgrad_x_input)

Layer-wise Relevance Propagation (LRP)

The LRP method was first introduced by Bach et al. (2015) and has a similar goal to the Gradient×Input approach explained in the last section: decompose the output into variable-wise relevances according to Eq. 1. The difference is that the prediction f(x)_c = y_c is redistributed layer by layer from the output node back to the inputs according to the weights and pre-activations. This is done by so-called relevance messages r_i ← j^{(l, l + 1)}, which can be defined by a rule on redistributing the upper-layer relevance R_j^l + 1 to the lower-layer R_i^l. In the package innsight, the following commonly used rules are defined (i is an index of a node in layer l and j an index of a node in layer l + 1):

• The simple rule (also known as LRP-0)
  This is the most basic rule on which all other rules are more or less based. The relevances are redistributed to the lower layers according to the ratio between local and global pre-activation. Let x_i the inputs, w_i, j the weights and b_j the bias vector of layer l and R_j^(l + 1) the upper-layer relevance, then the simple rule is defined as $$ r_{i \leftarrow j}^{(l, l+1)} = \frac{x_i\, w_{i,j}}{z_j} \, R_j^{l +1} \quad \text{with} \quad z_j = b_j + \sum_{k} x_k\, w_{k,j}. $$

• The ε-rule (also known as LRP-ϵ)
  One problem with the simple rule is that it is numerically unstable when the global pre-activation z_j vanishes and causes a division by zero. This problem is solved in the ε-rule by adding a stabilizer ε > 0 that moves the denominator away from zero, i.e., $$ r_{i \leftarrow j}^{(l, l+1)} = \frac{x_i\, w_{i,j}}{z_j + \text{sign}(z_j)\, \varepsilon}\, R_j^{l +1}. $$

• The α-β-rule (also known as LRP-αβ)
  Another way to avoid this numerical instability is by treating the positive and negative pre-activations separately. In this case, positive and negative values cannot cancel each other out, i.e., a vanishing denominator also results in a vanishing numerator. Moreover, this rule allows choosing a weighting for the positive and negative relevances, which is done with the parameters α, β ∈ ℝ satisfying α + β = 1. The α-β-rule is defined as $$ r_{i \leftarrow 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}\\ \text{with}\quad 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 nodes R_i^l is determined by summing up all incoming relevance messages r_i ← j^{(l, l + 1)} into the respective node of index i, i.e., R_i^l = ∑_jr_i ← j^{(l, l + 1)}.

Fig. 6: Layerwise Relevance Propagation

Analogous to the previous methods, the innsight method LRP inherits from the InterpretingMetod super class and thus all arguments. In addition, there are the following method-specific arguments for this method:

• rule_name (default: "simple"): This argument can be used to select the rule for the relevance messages. Implemented are the three rules described above, i.e., simple rule ("simple"), ε-rule ("epsilon") and α-β-rule ("alpha_beta"). However, a named list can also be passed to assign one of these three rules to each implemented layer type individually. Layers not specified in this list then use the default value "simple". For example, with list(Dense_Layer = "epsilon", Conv2D_Layer = "alpha_beta") the simple rule is used for all dense layers and the α-β-rule is applied to all 2D convolutional layers. The other layers not mentioned use the default rule. In addition, for normalization layers like 'BatchNorm_Layer', the rule "pass" is implemented as well, which ignores such layers in the backward pass. You can set the rule for the following layer types:

  • 'Dense_Layer', 'Conv1D_Layer', 'Conv2D_Layer', 'BatchNorm_Layer', 'AvgPool1D_Layer', 'AvgPool2D_Layer', 'MaxPool1D_Layer' and 'MaxPool2D_Layer'

• rule_param: The meaning of this argument depends on the selected rule. For the simple rule, for example, it has no effect. In contrast, this numeric argument sets the value of ε for the ε-rule and the value of α for the α-β-rule (remember: β = 1 − α). Passing NULL defaults to 0.01 for ε or 0.5 for α. Similar to the argument rule_name, this can also be a named list that individually assigns a rule parameter to each layer type.

• winner_takes_all: This logical argument is only relevant for models with a MaxPooling layer. Since many zeros are produced during the backward pass due to the selection of the maximum value in the pooling kernel, another variant is implemented, which treats a MaxPooling as an AveragePooling layer in the backward pass to overcome the problem of too many zero relevances. With the default value TRUE, the whole upper-layer relevance is passed to the maximum value in each pooling window. Otherwise, if FALSE, the relevance is distributed equally among all nodes in a pooling window.

# R6 class syntax
lrp <- LRP$new(converter, data,
  rule_name = "simple",
  rule_param = NULL,
  winner_takes_all = TRUE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# Using the helper function for initialization
lrp <- run_lrp(converter, data,
  rule_name = "simple",
  rule_param = NULL,
  winner_takes_all = TRUE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# We can analyze multiple inputs simultaneously
data <- matrix(
  c(
    0.49, # only neuron without bias term is activated
    0.6   # neuron with bias term is activated
  ), 
  ncol = 1
)

# Apply LRP with simple rule
lrp <- run_lrp(converter, data,
  ignore_last_act = FALSE
)
get_result(lrp)

# get approximation error
matrix(lrp$get_result()) - as_array(converter$model(torch_tensor(data))[[1]])

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. in 2017. It behaves similarly to LRP in a layer-by-layer backpropagation fashion from a selected output node back to the input variables. However, it incorporates a reference value 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(x)_c and the output of the reference value f(x̃)_c, i.e., f(x)_c − f(x̃)_c. This difference eliminates the bias term in the relevance messages so that no more relevance is absorbed and we have an exact variable-wise decomposition of Δy = f(x)_c − f(x̃)_c. In addition, the authors presented two rules to propagate relevances through the activation part of the individual layers, namely Rescale and RevealCancel rule. The Rescale rule simply scales the contribution to the difference from reference output according to the value of the activation function. The RevealCancel rule considers the average impact after adding the negative or positive contribution revealing dependencies missed by other approaches.

Analogous to the previous methods, the innsight method DeepLift inherits from the InterpretingMetod super class and thus all arguments. Alternatively, an object of the class DeepLift can also be created using the helper function run_deeplift(), which does not require prior knowledge of R6 objects. In addition, there are the following method-specific arguments for this method:

• x_ref (default: NULL): This argument describes the reference input x̃ for the DeepLift method. This value must have the same format as the input data of the passed model to the converter class, i.e.,

  • an array, data.frame, torch_tensor or array-like format of size (1, input_dim) or
  • a list with the corresponding input data (according to the upper point) for each of the input layers.
  • It is also possible to use the default value NULL to take only zeros as reference input.

• rule_name (default: 'rescale'): Name of the applied rule to calculate the contributions. Use either 'rescale' or 'reveal_cancel'.

• winner_takes_all: This logical argument is only relevant for MaxPooling layers and is otherwise ignored. With this layer type, it is possible that the position of the maximum values in the pooling kernel of the normal input x and the reference input x̃ may not match, which leads to a violation of the summation-to-delta property. To overcome this problem, another variant is implemented, which treats a MaxPooling layer as an AveragePooling layer in the backward pass only, leading to a uniform distribution of the upper-layer contribution to the lower layer.

# R6 class syntax
deeplift <- DeepLift$new(converter, data,
  x_ref = NULL,
  rule_name = "rescale",
  winner_takes_all = TRUE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# Using the helper function for initialization
deeplift <- run_deeplift(converter, data,
  x_ref = NULL,
  rule_name = "rescale",
  winner_takes_all = TRUE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# Create data
x <- matrix(c(0.55))
x_ref <- matrix(c(0.1))

# Apply method DeepLift with rescale rule
deeplift <- run_deeplift(converter, x, x_ref = x_ref, ignore_last_act = FALSE)

# Get result
get_result(deeplift)

library(neuralnet)
set.seed(42)

# Crate model with package 'neuralnet'
model <- neuralnet(Species ~ ., iris, hidden = 5, linear.output = FALSE)

# Step 1: Create 'Converter'
conv <- convert(model)

# Step 2: Apply DeepLift (reveal-cancel rule)
x_ref <- matrix(colMeans(iris[, -5]), nrow = 1) # use colmeans as reference value
deeplift <- run_deeplift(conv, iris[, -5],
  x_ref = x_ref, ignore_last_act = FALSE,
  rule_name = "reveal_cancel"
)

# Verify exact decomposition
y <- predict(model, iris[, -5])
y_ref <- predict(model, x_ref[rep(1, 150), ])
delta_y <- y - y_ref
summed_decomposition <- apply(get_result(deeplift), c(1, 3), FUN = sum) # dim 2 is the input feature dim

# Show the mean squared error
mean((delta_y - summed_decomposition)^2)

Integrated Gradients

In the Integrated Gradients method introduced by Sundararajan et al. (2017), the gradients are integrated along a path from the value x to a reference value x̃. This integration results, similar to DeepLift, in a decomposition of f(x) − f(x̃). In this sense, the method uncovers the feature-wise relative effect of the input features on the difference between the prediction f(x) and the reference prediction f(x̃). This is archived through the following formula: $$ \text{IntGrad}(x)_i^c = (x - \tilde{x}) \int_{\alpha = 0}^1 \frac{\partial f(\tilde{x} + \alpha (x - \tilde{x}))}{\partial x} d\alpha $$ In simpler terms, it calculates how much each feature contributes to a model’s output by tracing a path from a baseline input x̃ to the actual input x and measuring the average gradients along that path.

Similar to the other gradient-based methods, by default the integrated gradient is multiplied by the input to get an approximate decomposition of f(x) − f(x̃). However, with the parameter times_input only the gradient describing the output sensitivity can be returned.

Analogous to the previous methods, the innsight method IntegratedGradient inherits from the InterpretingMetod super class and thus all arguments. Alternatively, an object of the class IntegratedGradient can also be created using the helper function run_intgrad(), which does not require prior knowledge of R6 objects. In addition, there are the following method-specific arguments for this method:

• x_ref (default: NULL): This argument describes the reference input x̃ for the Integrated Gradients method. This value must have the same format as the input data of the passed model to the converter class, i.e.,

  • an array, data.frame, torch_tensor or array-like format of size (1, input_dim) or
  • a list with the corresponding input data (according to the upper point) for each of the input layers.
  • It is also possible to use the default value NULL to take only zeros as reference input.

• n (default: 50): Number of steps for the approximation of the integration path along α.

• times_input (default: TRUE): Multiplies the integrated gradients with the difference of the input features and the baseline values. By default, the original definition of Integrated Gradient is applied. However, by setting times_input = FALSE only an approximation of the integral is calculated, which describes the sensitivity of the features to the output.

# R6 class syntax
intgrad <- IntegratedGradient$new(converter, data,
  x_ref = NULL,
  n = 50,
  times_input = TRUE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# Using the helper function for initialization
intgrad <- run_intgrad(converter, data,
  x_ref = NULL,
  n = 50,
  times_input = TRUE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# Create data
x <- matrix(c(0.55))
x_ref <- matrix(c(0.1))

# Apply method IntegratedGradient
intgrad <- run_intgrad(converter, x, x_ref = x_ref, ignore_last_act = FALSE)

# Get result
get_result(intgrad)

Expected Gradients

The Expected Gradients method (Erion et al., 2021), also known as GradSHAP, is a local feature attribution technique which extends the Integrated Gradient method and provides approximate Shapley values. In contrast to Integrated Gradient, it considers not only a single reference value x̃ but the whole distribution of reference values X̃ ∼ x̃ and averages the Integrated Gradient values over this distribution. Mathematically, the method can be described as follows: $$ \text{ExpGrad}(x)_i^c = \mathbb{E}_{\tilde{x}\sim \tilde{X}, \alpha \sim U(0,1)} \left[(x - \tilde{x}) \times \frac{\partial f(\tilde{x} + \alpha (x - \tilde{x}))}{\partial x} \right] $$ These feature-wise values approximate a decomposition of the prediction minus the average prediction in the reference dataset, i.e., f(x) − 𝔼_x̃[f(x̃)]. This means, it solves the issue of choosing the right reference value.

Analogous to the previous methods, the innsight method ExpectedGradient inherits from the InterpretingMetod super class and thus all arguments. Alternatively, an object of the class ExpectedGradient can also be created using the helper function run_expgrad(), which does not require prior knowledge of R6 objects. In addition, there are the following method-specific arguments for this method:

• data_ref (default: NULL): This argument describes the reference inputs x̃ for the Expected Gradients method. This value must have the same format as the input data of the passed model to the converter class, i.e.,

  • an array, data.frame, torch_tensor or array-like format of size (1, input_dim) or
  • a list with the corresponding input data (according to the upper point) for each of the input layers.
  • It is also possible to use the default value NULL to take only zeros as reference input.

• n (default: 50): Number of samples from the distribution of reference values x̃ ∼ X̃ and number of samples for the approximation of the integration path along α.

# R6 class syntax
expgrad <- ExpectedGradient$new(converter, data,
  data_ref = NULL,
  n = 50,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# Using the helper function for initialization
expgrad <- run_expgrad(converter, data,
  x_ref = NULL,
  n = 50,
  ... # other arguments inherited from 'InterpretingMethod'
) 

library(neuralnet)
set.seed(42)

# Crate model with package 'neuralnet'
model <- neuralnet(Species ~ ., iris, linear.output = FALSE)

# Step 1: Create 'Converter'
conv <- convert(model)

# Step 2: Apply Expected Gradient
expgrad <- run_expgrad(conv, iris[c(1, 60), -5],
  data_ref = iris[, -5], ignore_last_act = FALSE,
  n = 10000
)

# Verify exact decomposition
y <- predict(model, iris[, -5])
delta_y <- y[c(1, 60), ] - rbind(colMeans(y), colMeans(y))
summed_decomposition <- apply(get_result(expgrad), c(1, 3), FUN = sum) # dim 2 is the input feature dim

# Show the error between both
delta_y - summed_decomposition

DeepSHAP

The DeepSHAP method (Lundberg & Lee, 2017) extends the DeepLift technique by not only considering a single reference value but by calculating the average from several, ideally representative reference values at each layer. The obtained feature-wise results are approximate Shapley values for the chosen output, where the conditional expectation is computed using these different reference values, i.e., the DeepSHAP method decompose the difference from the prediction and the mean prediction f(x) − 𝔼_x̃[f(x̃)] in feature-wise effects. This means, the DeepSHAP method has the same underlying goal as the Expected Gradient method and, hence, also solves the issue of choosing the right reference value for the DeepLift method.

Analogous to the previous methods, the innsight method DeepSHAP inherits from the InterpretingMetod super class and thus all arguments. Alternatively, an object of the class DeepSHAP can also be created using the helper function run_deepshap()`, which does not require prior knowledge of R6 objects. In addition, there are the following method-specific arguments for this method:

• data_ref (default: NULL): The reference data which is used to estimate the conditional expectation. These must have the same format as the input data of the passed model to the converter object. This means either

  • an array, data.frame, torch_tensor or array-like format of size (1, input_dim) or
  • a list with the corresponding input data (according to the upper point) for each of the input layers.
  • It is also possible to use the default value NULL to take only zeros as reference input.

• limit_ref (default: 100): This argument limits the number of instances taken from the reference dataset data_ref so that only random limit_ref elements and not the entire dataset are used to estimate the conditional expectation. A too-large number can significantly increase the computation time.

• (other model-specific arguments already explained in the DeepLift method, e.g., rule_name or winner_takes_all).

# R6 class syntax
deepshap <- DeepSHAP$new(converter, data,
  data_ref = NULL,
  limit_ref = 100,
  ... # other arguments inherited from 'DeepLift'
) 

# Using the helper function for initialization
deepshap <- run_deepshap(converter, data,
  data_ref = NULL,
  limit_ref = 100,
  ... # other arguments inherited from 'DeepLift'
) 

library(neuralnet)
set.seed(42)

# Crate model with package 'neuralnet'
model <- neuralnet(Species ~ ., iris, linear.output = FALSE)

# Step 1: Create 'Converter'
conv <- convert(model)

# Step 2: Apply Expected Gradient
deepshap <- run_deepshap(conv, iris[c(1, 60), -5],
  data_ref = iris[, -5], ignore_last_act = FALSE,
  limit_ref = nrow(iris)
)

# Verify exact decomposition
y <- predict(model, iris[, -5])
delta_y <- y[c(1, 60), ] - rbind(colMeans(y), colMeans(y))
summed_decomposition <- apply(get_result(deepshap), c(1, 3), FUN = sum) # dim 2 is the input feature dim

# Show the error between both
delta_y - summed_decomposition

Connection Weights

One of the earliest methods specifically for neural networks was the Connection Weights method invented by Olden et al. in 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 and the output node and then calculate the sum of all of them. However, this method ignores all bias vectors and all activation functions during calculation. Since only the weights are used, this method is independent of input data and, thus, a global interpretation method. In this package, we extended this method to a local one inspired by the method Gradient×Input (see here). Hence, the local variant is simply the point-wise product of the global Connection Weights method and the input data. You can use this variant by setting the times_input argument to TRUE and providing input data.

The innsight method ConnectionWeights also inherits from the super class InterpretingMethod, meaning that you need to change the term Method to ConnectionWeights. Alternatively, an object of the class ConnectionWeights can also be created using the helper function run_cw(), which does not require prior knowledge of R6 objects. The only model-specific argument is times_input, which can be used to switch between the global (FALSE) and the local (TRUE) Connection Weights method.

# The global variant (argument 'data' is no longer required)
cw_global <- ConnectionWeights$new(converter,
  times_input = FALSE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# The local variant (argument 'data' is required)
cw_local <- ConnectionWeights$new(converter, data,
  times_input = TRUE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# Using the helper function
cw_local <- run_cw(converter, data,
  times_input = TRUE,
  ... # other arguments inherited from 'InterpretingMethod'
) 

# Apply global Connection Weights method
cw_global <- run_cw(converter, times_input = FALSE)

# Show the result
get_result(cw_global)

# Create data
data <- array(c(0.1, 0.4, 0.6), dim = c(3, 1))

# Apply local Connection Weights method
cw_local <- run_cw(converter, data, times_input = TRUE)

# Show the result
get_result(cw_local)

Array (type = 'array')

In the simplest case, when the passed model to the converter object has only one input and one output layer, an R primitive array of dimension (batch_size, input_dim, outputs) is returned, where outputs means the number of elements from the argument output_idx. In addition, the passed or generated input and output names are added to the array.

However, this method behaves differently if the passed model has multiple input and/or output layers. In these cases, a list (or a nested list) with the corresponding input and output layers with the associated results is generated as in the simple case from before:

• n input layers and one output layer:

  List of n
    $ Input_1: array [batch_size, input_1_dim, outputs]
    $ Input_2: array [batch_size, input_2_dim, outputs]
    ...
    $ Input_n: array [batch_size, input_n_dim, outputs]

• one input layer and k output layers (outputs_i are the corresponding output indices specified previously in the argument output_idx):

  List of k
    $ Output_1: array [batch_size, input_dim, outputs_1]
    $ Output_2: array [batch_size, input_dim, outputs_2]
    ...
    $ Output_k: array [batch_size, input_dim, outputs_k]

• n input and k output layers (outputs_i are the corresponding output indices specified previously in the argument output_idx):

  List of k
    $ Output_1: List of n
      $ Input_1: array [batch_size, input_1_dim, outputs_1]
      $ Input_2: array [batch_size, input_2_dim, outputs_1]
      ...
      $ Input_n: array [batch_size, input_n_dim, outputs_1]
    $ Output_2: List of n
      $ Input_1: array [batch_size, input_1_dim, outputs_2]
      $ Input_2: array [batch_size, input_2_dim, outputs_2]
      ...
      $ Input_n: array [batch_size, input_n_dim, outputs_2]
    ...
    $ Output_k: List of n
      $ Input_1: array [batch_size, input_1_dim, outputs_k]
      $ Input_2: array [batch_size, input_2_dim, outputs_k]
      ...
      $ Input_n: array [batch_size, input_n_dim, outputs_k]

# Apply method 'Gradient x Input' for classes 1 ('setosa')  and 3 ('virginica')
tab_grad <- run_grad(tab_conv, tab_data,
  output_idx = c(1, 3),
  times_input = TRUE
)
# Get result
result_array <- tab_grad$get_result()
# You can also use the S3 function 'get_result'
result_array <- get_result(tab_grad)

# Show the result for datapoint 1 and 10
result_array[c(1, 10), , ]

# Apply method 'Gradient' for outputs 1  and 2
img_grad <- run_grad(img_conv, img_data, output_idx = c(1, 2))
# Get result
result_array <- img_grad$get_result()
# You can also use the S3 function 'get_result'
result_array <- get_result(img_grad)

# Show the result
str(result_array)

library(keras)

first_input <- layer_input(shape = c(10, 10, 2))
second_input <- layer_input(shape = c(11))
tmp <- first_input %>%
  layer_conv_2d(2, c(2, 2), activation = "relu") %>%
  layer_flatten() %>%
  layer_dense(units = 11)
output <- layer_add(c(tmp, second_input)) %>%
  layer_dense(units = 5, activation = "relu") %>%
  layer_dense(units = 3, activation = "softmax")
model <- keras_model(
  inputs = c(first_input, second_input),
  outputs = output
)

conv <- convert(model)
data <- lapply(
  list(c(10, 10, 2), c(11)),
  function(x) array(rnorm(5 * prod(x)), dim = c(5, x))
)

# Apply method 'Gradient' for outputs 1  and 2
grad <- run_grad(conv, data, output_idx = c(1, 2), channels_first = FALSE)
# Get result
result_array <- grad$get_result()
# You can also use the S3 function 'get_result'
result_array <- get_result(grad)

# Show the result
str(result_array)

library(keras)

first_input <- layer_input(shape = c(10, 10, 2))
second_input <- layer_input(shape = c(11))
tmp <- first_input %>%
  layer_conv_2d(2, c(2, 2), activation = "relu") %>%
  layer_flatten() %>%
  layer_dense(units = 11)
first_output <- layer_add(c(tmp, second_input)) %>%
  layer_dense(units = 20, activation = "relu") %>%
  layer_dense(units = 3, activation = "softmax")
second_output <- layer_concatenate(c(tmp, second_input)) %>%
  layer_dense(units = 20, activation = "relu") %>%
  layer_dense(units = 3, activation = "softmax")
model <- keras_model(
  inputs = c(first_input, second_input),
  outputs = c(first_output, second_output)
)

conv <- convert(model)
data <- lapply(
  list(c(10, 10, 2), c(11)),
  function(x) array(rnorm(5 * prod(x)), dim = c(5, x))
)

# Apply method 'Gradient' for outputs 1 and 2 in the first and
# for outputs 1 and 3 in the second output layer
grad <- run_grad(conv, data,
  output_idx = list(c(1, 2), c(1, 3)),
  channels_first = FALSE
)
# Get result
result_array <- grad$get_result()
# You can also use the S3 function 'get_result'
result_array <- get_result(grad)

# Show the result
str(result_array)

Torch Tensor (type = 'torch_tensor')

Internally all calculations are performed with the package torch. Therefore the results can also be returned as torch_tensor instead of an array analogous to the case above.

Data.Frame (type = 'data.frame')

Another convenient way to output the results, for example, to visualize them directly in your own ggplot2 object, is to return them as a data.frame. The following columns are generated:

• 'data': In this column the individual data points from the argument data are distinguished. Thereby the names 'data_1', 'data_2', … are generated according to the number of data points.

• 'model_input': This column describes the different input layers of the given model and creates the labels 'Input_1', 'Input_2', … according to the number of input layers.

• 'model_output': This column describes the different output layers of the given model and creates the labels 'Output_1', 'Output_2', … according to the number of output layers.

• 'feature': This column represents the input names of the corresponding input layer and is determined by the values from input_names passed the converter object. For image data, this column corresponds to the height of the image.

• 'feature'_2: This is only generated for image inputs and describes the image width.

• 'channel': This column is only generated for signal and image inputs describing the respective input channel.

• 'output_node': The respective output node or class for the corresponding output layer specified with the argument output_idx.

• 'value': This column contains the result of the applied method for the respective data point, input layer, output layer, output node/class and input feature.

• 'pred': The prediction of the respective input instance ('data') for the output node 'output_node' in the output layer 'model_output'.

• 'decomp_sum': The sum of all relevance values ('value') for an input instance ('data') and output node 'output_node' in the output layer 'model_output'.

• 'decomp_goal': The corresponding decomposition goal of the applied feature attribution method (if available, otherwise NAs) for an input instance ('data') and output node 'output_node' in the output layer 'model_output'.

• 'input_dimension': This column contains one of the values

  • 1: the value corresponds to an tabular input.
  • 2: the value corresponds to a signal input.
  • 3: the value corresponds to an image as input.

For example, if we train a model on the Iris dataset with the four inputs ("Sepal.Length", "Sepal.Width", "Petal.Length" and "Petal.Width") and the three classes "setosa", "versicolor" and "virginica" as outputs, we get the following data.frame:

head(get_result(tab_grad, "data.frame"), 5)

Analogously, you can also output the results for the model with image data. As already mentioned, the columns for the channel ("channel") and the image width ("feature_2") are then added:

head(get_result(img_grad, "data.frame"), 5)

library(ggplot2)
library(neuralnet)

# get the result from the tabular model
df <- get_result(tab_grad, "data.frame")

# calculate mean absolute gradient
df <- aggregate(df$value,
  by = list(feature = df$feature, class = df$output_node),
  FUN = function(x) mean(abs(x))
)

ggplot(df) +
  geom_bar(aes(x = feature, y = x, fill = class),
    stat = "identity",
    position = "dodge"
  ) +
  ggtitle("Mean over absolut values of the gradients") +
  xlab("Input feature") +
  ylab("Mean(abs(gradients))") +
  theme_bw()

Plots based on ggplot2

All ggplot2-based plots generated by the innsight package are an object of the S4 class innsight_ggplto2, a simple extension of a ggplot2 object that enables more detailed analysis and comparison of the results. In addition, it provides a way to visualize the results of models with multiple input or output layers with different scales, e.g., images and tabular data. The values in the slots of this class depend very much on whether results for multiple input or output layers (multiplot = TRUE), or for only a single input and output layer (multiplot = FALSE) are visualized as described in the following. But in general, however, the slots are only explained here for internal purposes and should not be changed in standard use cases.

# Create plot for output node 1 and 2 in the first output layer and
# data points 1 and 3
p <- plot(res_one_input, output_idx = c(1, 2), data_idx = c(1, 3))

# It's not an ggplot2 object
class(p)

# The slot 'grobs' only contains a single entry
p@grobs

# It's not a multiplot
p@multiplot

# Therefore, slots 'output_strips' and 'col_dims' are empty lists
p@output_strips
p@col_dims

|      Output 1: Node 1      |      Output 1: Node 3      |
|   Input 1   |   Input 2    |   Input 1   |   Input 2    |
|---------------------------------------------------------|-------------
|             |              |             |              |
| grobs[1,1]  |  grobs[1,2]  | grobs[1,3]  | grobs[1,4]   | data point 1
|             |              |             |              |
|---------------------------------------------------------|-------------
|             |              |             |              |
| grobs[2,1]  |  grobs[2,2]  | grobs[2,3]  | grobs[2,4]   | data point 2
| 

# Create a plot for output node 1 in the first output layer and node 2 in the
# second output layer and data points 1 and 3
p <- plot(res_two_inputs, output_idx = list(1, c(1, 2)), data_idx = c(1, 3))

# It's not a ggplot2 object
class(p)

# In this case, 'grobs' is a 2x6 matrix
p@grobs

# It's a multiplot
p@multiplot

# Slot 'output_strips' is a list with the three labels for the output nodes
# and the theme for the strips
str(p@output_strips, max.level = 1)

# Slot 'col_dims' contains the number of columns for each output node
p@col_dims

Plots based on plotly

As explained in the previous sections of Step 3, for each plot method, the argument as_plotly can be used to generate an interactive plot, which allows for more detailed analysis, especially for extensive models and results. However, these are created with the suggested package plotly, which should be installed in the R environment.
Analogous to the S4 class innsight_ggplot2, the class innsight_plotly is a simple extension of the library plotly instead of ggplot2. In this case, however, no distinction is made between a plot with multiple input or output layers except from additional strips. Hence, visualization can always be created in the same way:

1. The corresponding shapes and annotations of the slots shapes and annotations are added to each plot in plots. This also adds the strips for the output node (or input layer) at the top and, if necessary, the strips on the right side for the data points.

2. Subsequently, all individual plots are combined into a single plot object with the help of the function plotly::subplot (see ?plotly::subplot for details).

3. Lastly, the global elements from the layout slot are added and if there are multiple input layers (multiplot = TRUE), another output strip is added for the columns.

An example structure of the plot with multiple input layers is shown below:

|      Output 1: Node 1      |      Output 1: Node 3      |
|   Input 1   |   Input 2    |   Input 1   |   Input 2    |
|---------------------------------------------------------|-------------
|             |              |             |              |
| plots[1,1]  |  plots[1,2]  | plots[1,3]  | plots[1,4]   | data point 1
|             |              |             |              |
|---------------------------------------------------------|-------------
|             |              |             |              |
| plots[2,1]  |  plots[2,2]  | plots[2,3]  | plots[2,4]   | data point 2
| 

📝 Notes

• The plotly package is not a required package for innsight, so when creating interactive plotly-based plots, users should have it installed in their R environment on their own.

• The slots already mentioned and explained below, are for internal use only and generally should not be touched by users applying standard use cases.

• Due to the size of interactive plots, the following plotly-based plots are only rendered in the vignette at the pkdown website and not in the base R vignette!

Slots:

• plots : The individual plotly objects arranged as a matrix of size (N_datapoints, N_outputnodes ⋅ N_inputlayers)*.

• shapes : A list of two lists with the names shapes_strips and shapes_other. The list shapes_strips contains the shapes for the strips and may not be manipulated. The other list shapes_other contains a matrix of the same size as plots and each entry contains the shapes of the corresponding plot. For more details to shapes in plotly, we refer to the plotly website.

• annotations : A list of two lists with the names annotations_strips and annotations_other. The list annotations_strips contains the annotations for the strips and may not be manipulated. The other list annotations_other contains a matrix of the same size as plots and each entry contains the annotations of the corresponding plot. For more details to annotations in plotly, see ?plotly::add_annotations or the plotly website.

• multiplot : A logical value that indicates whether there are multiple input layers and thus whether another layer of strips should be created at the top of the plot.

• layout : This list contains all global layout options, e.g., update buttons, sliders, margins etc. (see ?plotly::layout for more details).

• col_dims : A list of two lists assigning a column index (col_idx) the corresponding output strip label (col_label).

# Create a plot for output node 1 in the first layer and output node 2 in the
# second layer and data point 1 and 3
p <- plot(res_two_inputs,
  output_idx = list(1, 2), data_idx = c(1, 3),
  as_plotly = TRUE
)

# Slot 'plots' is a 2x4 matrix (2 data points, 2 output nodes and 2 input layers)
p@plots

# Slot 'shapes' contains two 2x4 matrices with the corresponding shape objects
p@shapes

# The same for the annotations
p@annotations

# The model has multiple input layers, so the slot 'multiplot' is TRUE
p@multiplot

# The overall layout is stored in the slot 'layout'
str(p@layout, max.level = 1)

# 'col_dims' assigns the label of the additional strips to the respective column
p@col_dims

As described earlier, this S4 class generates a plotly object created by plotly::subplot, which can be treated and modified according to the syntax of the plotly package. Additionally, some generic functions are implemented to visualize individual aspects of the overall plot or to examine them in more detail. All available generic functions are listed below:

• [ : With this basic generic indexing function, the plots of individual rows and columns can be accessed and is returned as an object of the class innsight_plotly.

  Example

  # Create plot
  p <- plot(res_two_inputs,
    output_idx = list(1, 2), data_idx = c(1, 3),
    as_plotly = TRUE
  )
  
  # Show the whole plot
  p

  # Now you can select specific rows and columns for in-depth investigation,
  # e.g. only the result for output "Y2"
  p_new <- p[1:2, 3:4]
  
  # It's still an object of class 'innsight_plotly'
  class(p_new)
  #> [1] "innsight_plotly"
  #> attr(,"package")
  #> [1] "innsight"
  
  # Show the subplot
  p_new

• [[ : With this basic generic indexing function, you can select a single plot by passing the respective row and column index. This plot is returned as a plotly object.

  Example

  # Create plot
  p <- plot(res_two_inputs,
    output_idx = list(1, 2), data_idx = c(1, 3),
    as_plotly = TRUE
  )
  
  # Now you can select a single plot by passing the row and column index,
  # e.g. the plot for output "Y1", data point 3 and the second input layer
  p_new <- p[[2, 2]]
  
  # It's a plotly object
  class(p_new)
  #> [1] "plotly"     "htmlwidget"
  
  # Show the plot
  p_new

• plot, print and show : These methods provide generic visualization functions, which all behave the same in this case. They create the plot of the results and return a plotly object invisibly. For plot and print, you can also pass additional arguments forwarded to plotly::subplot.
  Examples

  # Create plot
  p <- plot(res_two_inputs,
    output_idx = list(1, 2), data_idx = 1,
    as_plotly = TRUE
  )
  
  # All methods behave the same and return a plotly object
  class(print(p))
  #> [1] "plotly"     "htmlwidget"
  class(show(p))
  #> [1] "plotly"     "htmlwidget"
  class(plot(p))
  #> [1] "plotly"     "htmlwidget"
  #>
  # You can also pass additional arguments to the method 'plotly::subplot',
  # e.g. the margins
  print(p, margin = c(0.03, 0.03, 0.01, 0.01))