# Chain-based Discriminative Autoencoders for Speech Recognition

###### Abstract

In our previous work, we proposed a discriminative autoencoder (DcAE) for
speech recognition. DcAE combines two training schemes into one. First, since
DcAE aims to learn encoder-decoder mappings, the squared error between the
reconstructed speech and the input speech is minimized. Second, in the code
layer, frame-based phonetic embeddings are obtained by minimizing the
categorical cross-entropy between ground truth labels and predicted triphone-
state scores. DcAE is developed based on the Kaldi toolkit by treating various
TDNN models as encoders. In this paper, we further propose three new versions
of DcAE. First, a new objective function that considers both categorical
cross-entropy and mutual information between ground truth and predicted
triphone-state sequences is used. The resulting DcAE is called a chain-based
DcAE (c-DcAE). For application to robust speech recognition, we further extend
c-DcAE to hierarchical and parallel structures, resulting in hc-DcAE and pc-
DcAE. In these two models, both the error between the reconstructed noisy
speech and the input noisy speech and the error between the enhanced speech
and the reference clean speech are taken into the objective function.
Experimental results on the WSJ and Aurora-4 corpora show that our DcAE models
outperform baseline systems.

Index Terms: discriminative autoencoder, robust speech recognition, multi-
condition training

##  1 Introduction

With the advent and breakthroughs of deep learning, current automatic speech
recognition (ASR) systems perform remarkably well in relatively clean
conditions. However, due to the mismatch between training and test conditions,
noise can greatly degrade performance in practical real-world application
scenarios. Numerous techniques have been proposed to reduce the adverse
effects of noise on acoustic modeling caused by interfering background events
(such as babble or street noise), reverberation generated from non-ideal room
acoustics, or channel distortion caused by different microphone
characteristics [1]. These techniques mainly fall into two categories: multi-
condition training (MCT) and speech enhancement (SE)-based pre-processing [2,
3, 4, 5, 6].

MCT is widely adopted to achieve model robustness against additive noise,
channel mismatch [7, 8, 2], and reverberation [9, 10]. As stated in [11],
adding extra noisy speech to the training data is a form of regularization
that can provide better generalization. In addition, MCT has been shown to
provide better performance than SE-based front-ends under unseen degradation
conditions [7, 8]. However, the performance under unseen conditions still lags
behind the performance obtained in seen conditions [8, 12].

On the other hand, empirical experience also shows that SE as a front-end can
improve ASR performance. Among the achievements in integrating SE into ASR
systems [5], Narayanan and Wang inherited the SE model in [13] and firstly
proposed joint training of a time-frequency (T-F) mask-based front-end and a
deep neural network (DNN)-based acoustic model back-end [14]. The front-end
was used to estimate clean features from noisy features, and the back-end was
used to classify phoneme labels. They were connected as a combined neural
network, and the model parameters were simultaneously optimized with back-
propagation. Related extensions can be seen in [15, 16, 17]. In addition, Soni
and Panda recently proposed another version of this type of framework in [18],
which can be used without parallel clean-noisy speech training data and
without pre-trained feature enhancement modules.

The above topology of integrated SE and phoneme classification can be simply
expressed by input features / enhancement block / classification block /
output targets. A detailed description can be seen in Fig. 2 of [19], which
shows how these blocks are trained by joint training (JT) and multi-task
learning (MTK). In this study, we extend our previous work, i.e.,
discriminative autoencoder (DcAE) [20, 21], and propose new neural
architectures, in which the classification block is wrapped in the enhancement
block for robust speech recognition. A denoising autoencoder has been a
popular DNN structure that can learn the mapping between noisy and clean
speech in the frequency or time domains [22, 23, 24, 25] from single- or
multi-channel observations [26, 27].

Figure 1: Different versions of DcAE. (a) Chain-based DcAE (c-DcAE), where the
code layer is divided into P-Code (phoneme-aware code) and R-Code (residual
code). P-Code directly contributes to senone scores. (b) Parallel c-DcAE (pc-
DcAE), where the predicted clean speech presumably comes from P-Code and
S-Code (speaker-aware code) through Decoder II, and the predicted noisy speech
presumably comes from P-Code, S-Code, and R-Code through Decoder I. (c)
Hierarchical c-DcAE (hc-DcAE), where C-Code (clean speech code) and R-Code are
used to generate noisy speech by Decoder I, and C-Code is further factorized
by Encoder II into P-Code and S-Code, which are used to generate clean speech
by Decode II. (d) Two types of Output and the associated losses.

In summary, our work has at least three major contributions:

  1. 1)

Compared to frame-level cross-entropy, few implementations have introduced the
lattice-free maximum mutual information (LF-MMI) criterion into common
autoencoder-based ASR or into joint training for robust acoustic modeling.

  2. 2)

Two variants of DcAE are proposed to integrate SE into acoustic modeling. What
they have in common is that the enhancement block not only removes noise, but
also restores the original signal. Their major difference lies in the way that
the enhancement block reconstructs noisy and clean speech: in parallel or
hierarchically. The new DcAE can not only extract cleaner and more
discriminative embeddings for triphone-state (senone) classification, but also
improve generalization. The decoder layers in DcAE can be regarded as a kind
of regularizer, as described in [20].

  3. 3)

Compared with the baseline model without decoder and code layers, DcAE
achieves relative word error rate (WER) reductions of 11.09% and 10.28% in
dev93 and eval92 on the WSJ corpus. Experimental results on the Aurora-4
corpus show that DcAE achieves a relative WER reduction of 5.86% over the
baseline model trained with MCT.

##  2 Proposed Neural Structures

An autoencoder is a learning model designed to convert an input into a
reconstructed output with minimal distortion. The embedded internal
representation, i.e., the code layer, contains the most critical information
that can be used to minimize the reconstruction error. In our previous work,
in order to obtain more accurate acoustic scores for ASR, we divided the code
layer into two sub-representations, the phoneme-aware code (P-Code) and the
residual code (R-Code), to form a basic version of DcAE [20]. P-Code contains
phonetic information and is connected to the original output layer of the
acoustic model of the ASR system. R-Code carries information unrelated to
phonetic content, which may include environmental noise, speaker identity, and
any other information deemed useless for the ASR task. P-Code and R-Code are
concatenated as input to the decoder to reconstruct the original acoustic
frame. This version of DcAE was implemented with the time-delay neural network
(TDNN) [28, 29] and its combination with long short-term memory (TDNN-LSTM)
[30] for acoustic modeling. The cross-entropy (CE) was used as the ASR loss.
In this section, we introduce the chain-based DcAE and its two variants, which
can make acoustic models more robust in noisy environments when clean-noisy
utterance pairs are sufficient for model training.

###  2.1 Chain-based DcAE (c-DcAE)

There are two main requirements for a model to be called “chain-based” in the
sense of the Kaldi toolkit111https://kaldi-asr.org/doc/chain.html [31, 32].
First, due to the reduced frame rate of the neural network output, an
unconventional hidden Markov model (HMM) topology that allows traversal of the
HMM in one state is used. Second, the sequence-level objective function, i.e.,
the log-probability of the correct sequence realized by MMI, is jointly
optimized with CE for model training. As shown in Fig. 1(a), the senone scores
(the black box) are compared with the reference senone labels using CE and LF-
MMI in the objective function.

###  2.2 Parallel Chain-based DcAE (pc-DcAE)

The first version of robust DcAE is the parallel chain-based DcAE (pc-DcAE).
As shown in Fig. 1(b), the architecture of pc-DcAE can be simply described as
Input / Encoder / Code / Decoder $\times$2 / Output $\times$3. A noisy
acoustic frame and its corresponding speaker embedding (e.g., i-vector) are
fed into the encoder, which can be implemented with various neural networks.
In the code layer (Code), the nodes are divided into three parts. The phoneme-
aware code (P-Code) that is designed to contain as much phonetic information
as possible is directly connected to the senone scores. S-Code contains the
speaker information while R-Code carries other remaining information,
especially unwanted noise. Decoder II is used to reconstruct the corresponding
clean acoustic frame of the input noisy acoustic frame according to the
concatenated P-Code and S-Code. Decoder I is used to reconstruct the noisy
acoustic frame according to the concatenated P-Code, S-Code, and R-Code. Both
decoders are implemented with fully-connected networks.

Note that S-Code is not formulated by training with actual speaker labels, but
it would make more sense if associated with a speaker-related loss. S-Code
represents speaker information, since ideal clean speech can be viewed as a
composite signal present only in orthogonal phonetic and speaker dimensions.

Table 1: WERs (relative changes over Baseline) (%) on WSJ. Both Baseline and Encoder of c-DcAE use TDNN-F. Spec denotes ratio of code sizes/depth of Decoder/weight for U-Net/no. of training epochs for warm start in c-DcAE (cf. Fig. 1). N means “not applied”. Model | Flat Start | Warm Start
---|---|---
dev93 | Spec | eval92 | Spec | dev93 | Spec | eval92 | Spec
Baseline | 4.42 | - | 2.53 | - | 4.42 | - | 2.53 | -
c-DcAE | 3.96 (10.41) | 0.5/2/N/N | 2.30 (9.09) | 0.5/6/N/N | 4.02 (9.05) | 0.5/2/N/4 | 2.34 (7.51) | 0.5/2/N/4
\+ U-Net (concat) | 3.96 (10.41) | 0.5/5/N/N | 2.30 (9.09) | 0.5/5/N/N | 4.02 (9.05) | 0.5/5/N/6 | 2.41 (4.74) | 0.5/5/N/4
\+ U-Net (sum) | 4.09 (7.47) | 0.5/2/0.8/N | 2.29 (9.49) | 0/3/0.3/N | 4.04 (8.60) | 0/3/0.3/6 | 2.30 (9.09) | 0/3/0.3/4
\+ U-Net (diff+concat) | 3.93 (11.09) | 0/3/0.3/N | 2.27 (10.28) | 0.5/3/0.5/N | 4.01 (9.28) | 0.5/3/0.5/5 | 2.36 (6.72) | 0/3/0.3/4

###  2.3 Hierarchical Chain-based DcAE (hc-DcAE)

The second version of robust DcAE is the hierarchical chain-based DcAE (hc-
DcAE). As shown in Fig. 1(c), the architecture of hc-DcAE can be simply
described as Input / Encoder / Code / Encoder / Code / Decoder / Output, where
the clean acoustic frame is hierarchically distilled from the noisy acoustic
frame through two code layers. There are two autoencoders lurking in hc-DcAE.
The first is Input (Noisy Speech) / Encoder I / Code I / Decoder I / Output
(Noisy Speech). By restoring the original noisy speech, it ensures that the
code layer (Code I) contains the salient components of the speech. The second
is C-Code / Encoder II / Code II / Decoder II / Output (Clean Speech), where
C-Code, the junction of the two autoencoders, is defined as the embedding of a
clean acoustic frame. Similar to pc-DcAE, both decoders in hc-DcAE are
implemented with fully-connected networks.

The main difference between pc-DcAE and hc-DcAE is twofold. In pc-DcAE, all
codes are generated separately in one operation by the same encoder, while in
hc-DcAE, each of the two encoders has its specific function. Therefore, the
model capacity of hc-DcAE may be larger that of pc-DcAE. Second, the
processing of clean speech in hc-DcAE during model training is deeper and more
complex.

###  2.4 Objective Functions

In addition to common ASR losses including CE and LF-MMI, the objective
function used in the proposed framework also involves reconstruction and
restoration errors.

####  2.4.1 Reconstruction and restoration errors

Suppose our proposed model $\mathcal{M}$ contains deterministic mappings
$f(\cdot)$ and $g(\cdot)$, which are responsible for observation
reconstruction and restoration, respectively. Given a set of noisy data
$\mathcal{X}$ and its corresponding clean data $\mathcal{Y}$ ready to go
through the inference-generation processes
$\mathcal{X}\overset{f}{\rightarrow}\mathcal{X}$ and
$\mathcal{X}\overset{g}{\rightarrow}\mathcal{Y}$, the average reconstruction
and restoration errors denoted by $\mathcal{L}_{rc}$ and $\mathcal{L}_{rs}$,
respectively, are given by

$\mathcal{L}_{rc}(\mathcal{X})=\frac{1}{|\mathcal{X}|}\sum_{\mathbf{x}\in\mathcal{X}}\|f(\mathbf{x})-\mathbf{x}\|^{2}_{2},$
(1)
$\mathcal{L}_{rs}(\mathcal{X})=\frac{1}{|\mathcal{X}|}\sum_{\mathbf{x}\in\mathcal{X},\mathbf{y}\in\mathcal{Y}}\|g(\mathbf{x})-\mathbf{y}\|^{2}_{2},$
(2)

where $\|\cdot\|^{2}_{2}$ is the 2-norm operator, and $|\mathcal{X}|$ is the
sample size. $f(\mathbf{x})$ and $g(\mathbf{x})$ correspond to Output (Noisy
Speech) and Output (Clean Speech) shown in Fig. 1, respectively.

####  2.4.2 ASR losses

There are two ASR losses related to Output (senone scores) in Fig. 1 in our
acoustic models. One, denoted by $\mathcal{L}_{CE}$, is the CE between the
distribution represented by the reference label and the predicted
distribution. The other, denoted by $-\mathcal{F}_{LF-MMI}$, is the negative
LF-MMI criterion between the distributions of predicted and reference word
sequences. See [31, 32, 33] for details on LF-MMI training (also known as
chain modeling).

By combining the criterion based on noisy feature reconstruction, clean
feature restoration, phoneme-aware CE, and sequence-level LF-MMI, the
objective function becomes

$\mathcal{L}=-\mathcal{F}_{LF-
MMI}+5\times\mathcal{L}_{CE}+\alpha\mathcal{L}_{rc}+\beta\mathcal{L}_{rs},$
(3)

where the weight 5 of $\mathcal{F}_{CE}$ follows most recipes for chain
modeling. $\alpha$ and $\beta$ are weights to increase or decrease the
regularization strength of $\mathcal{L}_{rc}$ and $\mathcal{L}_{rs}$,
respectively. They were determined heuristically in our experiments.

##  3 Experiments and Results

###  3.1 Datasets: WSJ and Aurora-4

We evaluated c-DcAE on clean speech recognition using the Wall Street Journal
(WSJ) dataset [34], which consists of WSJ0 (LDC93S6B) and WSJ1 (LDC94S13B).
The training set (called si284 in most Kaldi recipes) consists of 81 hours
speech, and we used dev93 and eval92 as test sets.

We evaluated pc-DcAE and hc-DcAE on robust speech recognition using the
Aurora-4 dataset. Aurora-4 [35] is a medium-level vocabulary task. The
transcriptions were based on WSJ0 (LDC93S6B) [34]. The dataset contains 16 kHz
speech data with additive noise and linear convolutional channel distortion,
which were synthetically introduced into clean speech. The training set
contains 7,138 utterance pairs from 83 speakers, where each pair consists of
clean speech and speech corrupted by one of six different noises (i.e.,
street, train station, car, babble, restaurant, and airport) at 10–20 dB SNR.
The test set was generated by the same types of noise and microphones and was
grouped into four subsets: clean, noisy, clean with channel distortion, and
noisy with channel distortion, which are referred to as A, B, C, and D,
respectively. The speech data were recorded by two microphones, one of which
was a near-field microphone. The speech recorded by this microphone can be
considered clean and noise-free. In contrast, the speech recorded by the
second microphone can be considered noisy speech.

Table 2: WERs (relative changes over Baseline) (%) on Aurora-4. Spec denotes depth of Decoder I/depth of Decoder II/weight for U-Net/depth of Encoder II (cf. Fig. 1). N means “not applied”. The total depth of Encoders I and II is equal to the depth of Baseline. Encoder | Method | U-Net | A | B | C | D | Average | Spec
---|---|---|---|---|---|---|---|---
TDNN-F | Baseline | - | 1.91 | 3.69 | 3.42 | 10.31 | 6.38 | -
pc-DcAE | - | 1.70 (10.99) | 3.50 (5.15) | 3.16 (7.60) | 10.22 (0.87) | 6.23 (2.35) | 5/5/N/N
sum | 1.66 (13.09) | 3.47 (5.96) | 3.21 (6.14) | 10.29 (0.19) | 6.25 (2.04) | 3/3/1/N
concat | 1.81 (5.24) | 3.63 (1.63) | 3.31 (6.14) | 10.37 (-0.58) | 6.37 (0.16) | 3/3/N/N
diff+concat | 1.70 (10.99) | 3.53 (4.34) | 3.03 (11.40) | 10.25 (0.58) | 6.24 (2.19) | 5/5/1/N
hc-DcAE | - | 1.61 (15.71) | 3.50 (5.15) | 3.14 (8.19) | 10.19 (1.16) | 6.21 (2.66) | 3/3/N/1
sum | 1.66 (13.09) | 3.53 (4.34) | 3.29 (3.80) | 10.36 (-0.48) | 6.31 (1.10) | 4/4/1/1
concat | 1.77 (7.33) | 3.55 (3.79) | 3.42 (0.00) | 10.29 (0.19) | 6.30 (1.25) | 3/3/N/1
diff+concat | 1.74 (8.90) | 3.55 (3.79) | 3.06 (10.53) | 10.04 (2.62) | 6.17 (3.29) | 3/3/0.8/2
CNN- TDNNF | Baseline | - | 1.72 | 3.19 | 3.40 | 9.48 | 5.80 | -
pc-DcAE | - | 1.57 (8.72) | 3.19 (0.00) | 3.03 (10.88) | 9.04 (4.64) | 5.57 (3.97) | 5/5/N/N
sum | 1.63 (5.23) | 3.19 (0.00) | 3.10 (8.82) | 8.93 (5.80) | 5.53 (4.66) | 5/5/1/N
concat | 1.64 (4.65) | 3.20 (-0.31) | 2.86 (15.88) | 9.05 (4.54) | 5.57 (3.97) | 5/5/N/N
diff+concat | 1.59 (7.56) | 3.26 (-2.19) | 2.99 (12.06) | 8.83 (6.86) | 5.51 (5.00) | 4/4/0.5/N
hc-DcAE | - | 1.63 (5.23) | 3.21 (-0.63) | 2.67 (21.47) | 8.82 (6.96) | 5.46 (5.86) | 4/4/N/3
sum | 1.57 (8.72) | 3.27 (-2.51) | 2.86 (15.88) | 8.97 (5.38) | 5.56 (4.14) | 5/5/0.8/3
concat | 1.72 (0.00) | 3.19 (0.00) | 3.10 (8.82) | 8.94 (5.70) | 5.54 (4.48) | 4/4/N/3
diff+concat | 1.59 (7.56) | 3.26 (-2.19) | 2.93 (13.82) | 8.91 (6.01) | 5.54 (4.48) | 5/5/1/1

###  3.2 Model Settings

All ASR systems were developed using a Gaussian mixture model (GMM)/DNN/HMM-
based process based on a standard Kaldi recipe222https://github.com/kaldi-
asr/kaldi/blob/master/egs/{wsj/s5/run.sh,aurora4/s5/run.sh}. The standard
four-gram and trigram language models were used in the WSJ and Aurora-4 tasks,
respectively.

Our DcAE models were built on top of two acoustic models, factorized TDNN
(TDNN-F) [36] and its convolutional neural network-based variant (CNN-TDNNF).
TDNN-F consists of 13 hidden layers, each containing 1024-dimensional hidden
nodes and 128-dimensional linear bottleneck nodes. In addition, the parameters
of the final layer are also factorized using a 192-node linear bottleneck
layer. CNN-TDNNF consists of 15 layers. The first six layers are convolutional
and the remaining nine layers are TDNN-F. The time offsets and height offsets
were set to “-1, 0, 1” for all convolutional layers. The numbers of filters
for the six convolutional layers were 48, 48, 64, 64, 64, and 128,
respectively. The parameters of the TDNN-F layers were set to the same value
as TDNN-F.

Forty high-resolution MFCCs (obtained by discrete cosine transform of 40 Mel-
frequency bins and normalized by utterance-based mean subtraction) and the
100-dimensional i-vector [37] were used as frame input for all systems.

###  3.3 Results

The structure of c-DcAE provides an opportunity for model pre-training using
acoustic features only (cf. warm start in Table 1). Experimental results on
WSJ show that training from scratch using acoustic features and phoneme labels
(cf. flat start) works better. However, we found that due to some training
schemes used in Kaldi, the training error defined in Eq. 1, although
eventually converged, did not drop much. Although worse than training from
scratch, warm start still shows its potential and superiority to the baseline
by reducing the WER by 9.28% and 9.09% in dev93 and eval92, respectively. It
is worth noting that it should make more sense to pre-train on large amounts
of unlabeled data. We will investigate this further in the future. Due to
space limitations, we do not present the results of CNN-TDNNF in Table 1, but
instead compare various U-Net operations [38] applied in c-DcAE (cf. Fig. 1).
As shown in Table 1, c-DcAE+U-Net(diff+concat), where the layers in Decoder
were concatenated by the difference between the layers of Encoder and Decoder,
achieved the best performance and reduced the WER by 11.09% and 10.28% in
dev93 and eval92, respectively. Note that the weight for U-Net denotes the
scaling strength coming from Encoder (for “sum”) or the difference between
Encoder and Decoder (for “diff+concat”), and the ratio of code sizes is
defined by “$\text{the size of R-Code}/\text{the size of P-Code}$”. The size
of P-Code was fixed to 1,024.

In robust speech recognition experiments on Aurora-4, the sizes of P-Code,
S-Code and R-Code were set to 1,024 in pc-DcAE and hc-DcAE. From Table 2, we
can see that pc-DcAE (diff+concat) achieved the best performance and reduced
the WER by 11.40% in Set C (channel-distorted clean speech), when TDNN-F was
used as Encoder. The CNN-TDNNF-based hc-DcAE could even reduce the WER by
21.47% in Set C. This result shows that DcAE is more robust to channel
mismatch than Baseline. An interesting observation is that both pc-DcAE and
hc-DcAE had less improvement over Baseline under noisy conditions (Sets B and
D). CNN-TDNNF-based pc-DcAE and hc-DcAE performed even worse than Baseline in
Set B. A possible reason might be that CNN-TDNNF itself is robust enough to
additive noise, since the first few convolutional layers can filter out non-
speech sounds to some extent. This somewhat conforms with the fact that most
Transformer-based end-to-end ASR systems tend to use CNN for front-end
subsampling [39, 40]. Furthermore, the benefit that CNN-TDNNF brings to pc-
DcAE and hc-DcAE is that all modifications by U-Net seem unnecessary. As can
be seen from the last 4 rows of Table 2, naive hc-DcAE achieved the best
performance with a relative reduction of 5.86% in average WER. This result is
also better than most studies in the literature thus far, except for CNN-Raw
[41], which achieves a WER of 5.10%333In CNN-Raw, the forced alignment of
clean speech was directly applied to the corresponding noisy speech.
Considering that clean speech corresponding to noisy speech is not always
available, we performed forced alignment on clean speech and noisy speech
separately..

##  4 Conclusion

In this paper, we have proposed a neural architecture, the chain-based
discriminative autoencoder (c-DcAE), which successfully combines chain
modeling and DcAE in acoustic modeling. For robust ASR, two variants of
c-DcAE, namely pc-DcAE and hc-DcAE, are further proposed to process noisy and
clean speech in parallel and hierarchical manners, respectively. Joint
training involves ASR-related losses, especially lattice-free MMI, and feature
reconstruction/restoration errors. Several types of connections using U-Net
(like skip-connections) are introduced into the structure of DcAE. Another
highlight of our DcAE-based models is that they do not require additional
computation time than the baseline in the inference stage, as all decoders in
DcAE can be removed at inference time. Our experimental results have shown
that the proposed models outperform the baseline models on the WSJ and
Aurora-4 tasks.

## References

  * [1] V. Mitra, H. Franco, C. Bartels, J. v. Hout, M. Graciarena, and D. Vergyri, “Speech recognition in unseen and noisy channel conditions,” in _Proc. ICASSP_ , 2017.
  * [2] D. Yu, M. L. Seltzer, J. Li, J.-T. Huang, and F. Seide, “Feature learning in deep neural networks - studies on speech recognition tasks,” in _Proc. ICLR_ , 2013.
  * [3] Y. Ephraim and D. Malah, “Speech enhancement using a minimum mean-square error log-spectral amplitude estimator,” _IEEE Trans. Audio, Speech, Sinal Process._ , vol. 33, no. 2, pp. 443–445, 1985.
  * [4] T. Nakatani, T. Yoshioka, K. Kinoshita, M. Miyoshi, and B. H. Juang, “Speech dereverberation based on variance-normalized delayed linear prediction,” _IEEE Trans. Audio, Speech, Lang. Process._ , vol. 18, no. 7, pp. 1717–1731, 2010.
  * [5] Y. Xu, J. Du, L.-R. Dai, and C.-H. Lee, “An experimental study on speech enhancement based on deep neural networks,” _IEEE Signal Process. Lett._ , vol. 21, no. 1, pp. 65–68, 2014.
  * [6] A. L. Maas, Q. V. Le, T. M. O. Neil, O. Vinyals, P. Nguyen, and A. Y. Ng, “Recurrent neural networks for noise reduction in robust ASR,” in _Proc. Interspeech_ , 2012.
  * [7] M. L. Seltzer, D. Yu, and Y. Wang, “An investigation of deep neural networks for noise robust speech recognition,” in _Proc. ICASSP_ , 2013.
  * [8] J. Du, Q. Wang, T. Gao, Y. Xu, L. Dai, and C.-H. Lee, “Robust speech recognition with speech enhanced deep neural networks,” in _Proc. Interspeech_ , 2014.
  * [9] L. Couvreur, C. Couvreur, and C. Ris, “A corpus-based approach for robust ASR in reverberant environments,” in _Proc. ICSLP_ , 2000.
  * [10] R. Hsiao, J. Ma, W. Hartmann, M. Karafiat, F. Grezl, L. Burget, I. Szoke, J. Cernocky, S. Watanabe, Z. Chen, S. Mallidi, H. Hermansky, S. Tsakalidis, and R. Schwartz, “Robust speech recognition in unknown reverberant and noisy conditions,” in _Proc. IEEE ASRU_ , 2015.
  * [11] S. Yin, C. Liu, Z. Zhang, Y. Lin, D. Wang, J. Tejedor, T. Fang, and Y. Li, “Noisy training for deep neural networks in speech recognition,” _Eurasip J. Audio Speech Music Process._ , 2015.
  * [12] Y. Qian, M. Yin, Y. You, and K. Yu, “Multi-task joint-learning of deep neural networks for robust speech recognition,” in _Proc. IEEE ASRU_ , 2015.
  * [13] Y. Wang and D. Wang, “Towards scaling up classification-based speech separation,” _IEEE Trans. Audio, Speech, Lang. Process._ , vol. 21, no. 7, pp. 1381–1390, 2013.
  * [14] A. Narayanan and D. Wang, “Joint noise adaptive training for robust automatic speech recognition,” in _Proc. ICASSP_ , 2014.
  * [15] A. Narayanan, A. Misra, and K. Chin, “Large-scale, sequence-discriminative, joint adaptive training for masking-based robust ASR,” in _Proc. Interspeech_ , 2015.
  * [16] A. Narayanan and D. Wang, “Improving robustness of deep neural network acoustic models via speech separation and joint adaptive training,” _IEEE Trans. Audio, Speech, Lang. Process._ , vol. 23, no. 1, pp. 92–101, 2015.
  * [17] T. Gao, J. Du, L. R. Dai, and C.-H. Lee, “Joint training of front-end and back-end deep neural networks for robust speech recognition,” in _Proc. ICASSP_ , 2015.
  * [18] M. Soni and A. Panda, “Label driven time-frequency masking for robust continuous speech recognition,” in _Proc. Interspeech_ , 2019.
  * [19] M. Fujimoto and H. Kawai, “One-pass single-channel noisy speech recognition using a combination of noisy and enhanced features,” in _Proc. Interspeech_ , 2019.
  * [20] P.-T. Huang, H.-S. Lee, S.-S. Wang, K.-Y. Chen, Y. Tsao, and H.-M. Wang, “Exploring the encoder layers of discriminative autoencoders for LVCSR,” in _Proc. Interspeech_ , 2019.
  * [21] M.-H. Yang, H.-S. Lee, Y.-D. Lu, K.-Y. Chen, Y. Tsao, B. Chen, and H.-M. Wang, “Discriminative autoencoders for acoustic modeling,” in _Proc. Interspeech_ , 2017.
  * [22] S. Pascual, A. Bonafonte, and J. Serra, “SEGAN: speech enhancement generative adversarial network,” in _Proc. Interspeech_ , 2017.
  * [23] D. Rethage, J. Pons, and X. Serra, “A wavenet for speech denoising,” in _Proc. ICASSP_ , 2018.
  * [24] N. Tawara, T. Kobayashi, M. Fujieda, K. Katagiri, T. Yazu, and T. Ogawa, “Adversarial autoencoder for reducing nonlinear distortion,” in _Proc. APSIPA ASC_ , 2018.
  * [25] K. Kinoshita, T. Ochiai, M. Delcroix, and T. Nakatani, “Improving noise robust automatic speech recognition with single-channel time-domain enhancement network,” in _Proc. ICASSP_ , vol. 2020-May, no. 4, 2020, pp. 7009–7013.
  * [26] X. Lu, Y. Tsao, S. Matsuda, and C. Hori, “Speech enhancement based on deep denoising autoencoder,” in _Proc. Interspeech_ , 2013.
  * [27] S. Araki, T. Hayashi, M. Delcroix, M. Fujimoto, K. Takeda, and T. Nakatani, “Exploring multi-channel features for denoising-autoencoder-based speech enhancement,” in _Proc. ICASSP_ , 2015.
  * [28] A. Waibel, T. Hanazawa, G. Hinton, K. Shikano, and K. J. Lang, “Phoneme recognition using time-delay neural networks,” _IEEE Trans. Acoustics, Speech, and Signal. Process._ , vol. 37, no. 3, pp. 328–339, 1989.
  * [29] V. Peddinti, D. Povey, and S. Khudanpur, “A time delay neural network architecture for efficient modeling of long temporal contexts,” in _Proc. Interspeech_ , 2015.
  * [30] V. Peddinti, Y. Wang, D. Povey, and S. Khudanpur, “Low latency acoustic modeling using temporal convolution and LSTMs,” _IEEE Signal Process. Lett._ , vol. 25, pp. 373–377, 2018.
  * [31] K. Vesely, A. Ghoshal, L. Burget, and D. Povey, “Sequence-discriminative training of deep neural networks,” in _Proc. Interspeech_ , 2013.
  * [32] D. Povey, V. Peddinti, D. Galvez, P. Ghahrmani, V. Manohar, X. Na, Y. Wang, and S. Khudanpur, “Purely sequence-trained neural networks for ASR based on lattice-free MMI,” in _Proc. Interspeech_ , 2016.
  * [33] H. Hadian, H. Sameti, D. Povey, and S. Khudanpur, “End-to-end speech recognition using lattice-free MMI,” in _Proc. Interspeech_ , 2018.
  * [34] D. B. Paul and J. M. Baker, “The design for the wall street journal-based CSR corpus,” in _Proc. DARPA Speech and Natural Language Workshop_ , 1992.
  * [35] N. Parihar and J. Picone, “Aurora working group: DSR front end LVCSR evaluation AU/384/02,” Inst. for Signal & Inform. Process., Mississippi State Univ., Tech. Rep., 2002. [Online]. Available: http://aurora.hsnr.de/aurora-4.html
  * [36] D. Povey, G. Cheng, Y. Wang, K. Li, H. Xu, M. Yarmohamadi, and S. Khudanpur, “Semi-orthogonal low-rank matrix factorization for deep neural networks,” in _Proc. Interspeech_ , 2018.
  * [37] N. Dehak, P. J. Kenny, R. Dehak, P. Dumouchel, and P. Ouellet, “Front-end factor analysis for speaker verification,” _IEEE Trans. Audio, Speech, Lang. Process._ , vol. 19, no. 4, pp. 788–798, 2011.
  * [38] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _Proc. MICCAI_ , 2015.
  * [39] S. Watanabe, T. Hori, S. Karita, T. Hayashi, J. Nishitoba, Y. Unno, N. E. Y. Soplin, J. Heymann, M. Wiesner, N. Chen, A. Renduchintala, and T. Ochiai, “ESPNet: End-to-end speech processing toolkit,” in _Proc. Interspeech_ , 2018.
  * [40] A. Gulati, J. Qin, C.-C. Chiu, N. Parmar, Y. Zhang, J. Yu, W. Han, S. Wang, Z. Zhang, Y. Wu, and R. Pang, “Conformer: Convolution-augmented transformer for speech recognition,” in _Proc. Interspeech_ , 2020.
  * [41] E. Loweimi, P. Bell, and S. Renals, “On the robustness and training dynamics of raw waveform models,” in _Proc. Interspeech_ , 2020.


11institutetext: Institute for Medical Informatics, Statistics and
Documentation.
Medical University Graz, Austria
<EMAIL_ADDRESS>
22institutetext: MI2.AI, Warsaw University of Technology, Poland
33institutetext: MI2.AI, University of Warsaw, Poland 44institutetext: Human-
Centered AI Lab, University of Natural Resources and Life Sciences, Vienna,
Austria

# Explaining and visualizing black-box models through counterfactual paths

Bastian Pfeifer 11 Mateusz Krzyzinski 22 Hubert Baniecki 33 Anna Saranti 1144
Andreas Holzinger 1144 Przemyslaw Biecek
2233

###### Abstract

Explainable AI (XAI) is an increasingly important area of machine learning
research, which aims to make black-box models transparent and interpretable.
In this paper, we propose a novel approach to XAI that uses the so-called
_counterfactual paths_ generated by conditional permutations of features. The
algorithm measures feature importance by identifying sequential permutations
of features that most influence changes in model predictions. It is
particularly suitable for generating explanations based on counterfactual
paths in knowledge graphs incorporating domain knowledge. _Counterfactual
paths_ introduce an additional graph dimension to current XAI methods in both
explaining and visualizing black-box models. Experiments with synthetic and
medical data demonstrate the practical applicability of our approach.

###### Keywords:

explainable machine learning, knowledge graph, feature importance,
counterfactual explanation

##  1 Introduction

Figure 1: Counterfactual path. Shown is the core principle of our proposed
algorithm. The input data is in tabular form. A certain $k$ number of features
is sampled, forming a candidate counterfactual path. These features are then
subsequently permuted according to the ordering reflected by the candidate
path. The swapped classes are tracked, and in case they exceed a certain
threshold, the permutation procedure terminates and the detected
counterfactual path is stored. According to the edges of that path a
transition matrix ($\mathbf{T}$) is formed, from which feature importances are
derived.

Explainable AI (XAI) methods are becoming a promising solution to the emerging
challenge of providing human-understandable justifications for AI decisions
(Holzinger et al., 2022), e.g. in medicine (Gozzi et al., 2022; Xu et al.,
2022; Krzyziński et al., 2023), biology (Anguita-Ruiz et al., 2020) or finance
(Bücker et al., 2022). One of the key problems in model explanation is
determining the importance of variables, and to this end, many methods have
been proposed either specific to some model families or model agnostic, here
we focus on the latter. One approach to XAI is through counterfactual
explanations, which involve generating alternative scenarios that could have
led to different prediction outcomes. Counterfactual explanations help users
understand how an AI model arrived at a particular decision and what factors
influenced that decision (Chou et al., 2022). Another widely adopted approach
is feature importance methods that express the reliance of model predictive
performance on a specific feature in data (Fisher et al., 2019; Casalicchio et
al., 2019).

However, these methods most often present the importance for variables
separately, without taking into account, often complex relationships between
variables such as interactions or correlations. To challenge this status quo,
we introduce the counterfactual paths (CPATH) algorithm for model-agnostic
global explanations of machine learning predictive models trained on tabular
data. This algorithm is inspired by both counterfactual explanations and
permutational feature importance (see Figure 1 for a graphical illustration).
Compared to classical feature importance methods, CPATH provides additional
graph information about the counterfactual dependence of the black-box model
on particular features, and thus can help to uncover its decision-making
process.

Counterfactual paths represent the relationships between input features and
the model’s output, allowing users to explore how changing individual features
or their combinations affect the model’s predictions. One of the main
advantages of counterfactual paths is their ability to provide insights into
how the model works and to identify potential biases or confounding factors
that may affect its predictions. These insights can then improve the model’s
accuracy and robustness.

Furthermore, counterfactual paths and their visualization could provide a more
intuitive and interpretable explanation of the model’s behaviour than
traditional feature importance methods. For example, rather than simply
providing a list of the top features that contribute to the model’s
predictions, counterfactual paths can show how changes to specific
combinations of features lead to changes in the output of the model. This can
help users to better understand the underlying patterns in the data and the
overall behaviour of the model.

##  2 Related work

####  2.0.1 Model-agnostic feature importance.

Explaining a predictive model on a _global_ level aims to understand how
important is a given feature to its performance. Several model-agnostic, i.e.
explaining any black-box function, feature importance measures have been
proposed. The widely adopted approach is permutation feature importance
(Fisher et al., 2019), which is also extended to local and partial importance
(Casalicchio et al., 2019). Molnar et al. (2021) introduce confidence
intervals for permutation feature importance and Au et al. (2022) propose to
group features to explain their combined importance. In practice, estimating
the marginal importance of a single feature without taking into account
correlation structure in data is a challenge (Molnar et al., 2023). To address
feature dependence, Watson and Wright (2021) proposes to measure the
conditional predictive impact between features and predictions using the
knockoff sampling framework; also for categorical features (Blesch et al.,
2023). Related is work on conditional estimation of Shapley-based feature
attributions Aas et al. (2021). As opposed to considering Shapley-based
feature importance measures (Casalicchio et al., 2019; Covert et al., 2020),
in this paper, we specifically relate permutation importance to counterfactual
explanations.

####  2.0.2 Counterfactual explanations.

A counterfactual explanation of model prediction describes the smallest change
to the feature values that change the predicted class. It is a useful _local_
“what-if” explanation for making actionable decisions Wachter et al. (2017),
Saranti et al. (2022). One can optimize to find counterfactuals using model-
specific optimization methods, e.g. gradients. Our work relates more to model-
agnostic approaches available to any black-box model (Karimi et al., 2020).
Dandl et al. (2020) propose multi-objective counterfactual explanations that
take into account data manifold, i.e. how likely it is that the counterfactual
data point originates from the training data distribution. Mothilal et al.
(2020) focus on finding a diverse set of counterfactual data points for a
given prediction. Most recent work considers the robustness of such
explanations to (potentially adversarial) data perturbations (Pawelczyk et
al., 2023). Contrary to related work, we use the intuition behind
counterfactual explanations to construct a novel global explanation of feature
importance.

####  2.0.3 Domain knowledge in the context of model explanation.

Relating explanations to domain knowledge is an emerging research topic
(Biecek and Burzykowski, 2021). Crucially, interpreting models should be done
with respect to data distribution and its correlation structure (Molnar et
al., 2021; Baniecki et al., 2023). One idea is to build a surrogate model
based on explanations of a black-box function achieving an interpretable
predictive interface consistent with domain knowledge (Alaa et al., 2021;
Gosiewska et al., 2021). Other approaches consider including domain knowledge
directly in the algorithm used to learn a predictive function (Panigutti et
al., 2020; Confalonieri et al., 2021; Pfeifer et al., 2022a), which is
challenging to do in an algorithm-agnostic way.

In this paper, we introduce feature importance explanation and visualization
through counterfactual paths that provide additional information about the
graph structure of features. Domain knowledge can both be derived from
explanations, as well as participate in estimating more accurate explanations.

##  3 Counterfactual paths for explaining black-box models

####  3.0.1 Intuition.

The herein proposed CPATH algorithm randomly selects paths through the feature
space, where each feature represents a node in a fully-connected graph or a
user-specified knowledge graph, masking the sampling scheme. Once a path is
created, CPATH permutes one feature after the other of the sampled path and
terminates as soon as a certain number of class labels swap to the inverse
class. We call these paths counterfactual paths. From these counterfactual
paths, we derive an adjacency matrix weighted by the inferred path lengths.
Based on the weighted adjacency matrix the global importance of a certain
feature is determined.

###  3.1 Mathematical formulation

Let $\mathcal{M}\colon\mathbb{R}^{p}\to\\{1,\ldots,g\\}$, where $g$ is the
number of classes, denote the model of interest, which is a classifier. Given
an observation $\mathbf{x}=(x_{1},x_{2},\ldots,x_{p})\in\mathbb{R}^{p}$, the
model’s prediction is denoted as $\mathcal{M}(\mathbf{x})$.

We describe the problem using weighted directed graphs, (weighted digraphs). A
graph $G_{\mathcal{M}}=(V,E,w_{\mathcal{M}})$ is defined as an ordered triple,
where $V={1,\ldots,p}$ represents the set of predictors (explanatory
variables) that form the vertices of the graph. The set of directed edges
(arcs), denoted as $E$, contains ordered pairs of vertices $(i,j)$ such that
there exists a directed edge from vertex $i$ to vertex $j$. The function
$w_{\mathcal{M}}\colon E\to\mathbb{R}$ maps the edges to their weights and can
also be represented as an adjacency matrix $W$.

The variables used in the $\mathcal{M}$ model which constitute the set of
vertices $V$ are predetermined, i.e., these are the variables used in the
model fitting process. Similarly, the set $E$ is known – its form can be
derived from a domain knowledge graph describing the process of generating the
data. However, if the domain knowledge graph is unknown, it can be assumed
that $G_{\mathcal{M}}$ is a complete digraph, meaning that each pair of
vertices is connected by a symmetric pair of directed arcs.

The objective is to estimate the function $w_{\mathcal{M}}$ (related to the
adjacency matrix $W$) by aggregating _counterfactual paths_ obtained from
sampling trajectories of finite-length random walks on the unweighted version
of the graph $\widetilde{G_{\mathcal{M}}}$.

######  Definition 1 (Counterfactual path)

Given a model’s predictions
$\mathcal{M}(\mathbf{X})=[\mathcal{M}(\mathbf{x}_{1}),\ldots,\mathcal{M}(\mathbf{x}_{n})]\in\\{1,\ldots,g\\}^{n}$
for a dataset $\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})^{T}$, we call
a sequence $\mathbf{v}=(v_{1},\ldots,v_{k})$, $k\leq p$, a counterfactual path
under counterfactual policy $\Psi$ if perturbing (e.g., by permuting) values
in the features $v_{1},\ldots,v_{k}$ leads to a substantial change in the
model’s predictions according to the indicator given by $\Psi$, i.e, if

$\Psi(\mathcal{M}(\mathbf{X}),\mathcal{M}(\mathbf{X}^{\prime}_{\mathbf{v}}))=1.$

In the above definition, $\mathbf{X}^{\prime}_{\mathbf{v}}$ represents a
dataset with perturbed values of features from the sequence $\mathbf{v}$ (note
that perturbations are performed sequentially, for individual variables as
they are added to the path). The counterfactual policy $\Psi$ is used to
determine whether a perturbation of features leads to a significant change in
predictions.

Various choices for the counterfactual policy $\Psi$ are possible. One
approach is to define $\Psi$ based on the fraction of changed predictions $p$,
where

$p=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}(\mathcal{M}(\mathbf{x}_{i})\neq\mathcal{M}(\mathbf{x}^{\prime}_{i})).$

With that, we define $\Psi$ in a stochastic way as follows:

$\Psi_{\mathrm{s}}(\mathcal{M}(\mathbf{X}),\mathcal{M}(\mathbf{X}^{\prime}_{\mathbf{v}}))\coloneqq
X\sim Bern(p).$

Alternatively, $\Psi$ can be defined as an indicator of whether the fraction
of changed predictions exceeds the selected threshold $\kappa$ chosen by a
user, i.e.,

$\Psi_{\mathrm{d}}(\mathcal{M}(\mathbf{X}),\mathcal{M}(\mathbf{X}^{\prime}_{\mathbf{v}}))\coloneqq\mathbbm{1}(p>\kappa).$

Another option is to penalize path lengths by considering their length $k$ in
the definition of $\Psi$.

The selection of an appropriate counterfactual policy $\Psi$ is a crucial step
in the proposed methodology. Once the counterfactual policy is defined,
Algorithm 1 is employed to identify and collect counterfactual paths. This
algorithm takes several inputs, including the model of interest $\mathcal{M}$,
the dataset $\mathbf{X}$, the counterfactual policy $\Psi$, the unweighted
graph $\widetilde{G_{\mathcal{M}}}$, the number of iterations $n_{iter}$, and
the maximal path length $k$. The CPATH algorithm initializes a set to store
counterfactual paths and iterates a specified number of times. In each
iteration, a new path is generated by sampling vertices from the unweighted
graph and checking for significant changes in the model’s predictions based on
the counterfactual policy.

Data: $\mathcal{M}$, $\mathbf{X}$, $\Psi$, $\widetilde{G_{\mathcal{M}}}$,
$n\\_iter\in\mathbf{Z_{+}}$, $k\in\mathbf{Z_{+}}$

1 $CPATHS\leftarrow\\{\\}$;

2 for  _$i=1$ to $n\\_iter$_ do

3 $\mathbf{v}\leftarrow[\,]$;

4 for  _$j=1$ to $k$_ do

5 $v_{j}\leftarrow\mathrm{sample\\_vertex}(\widetilde{G_{\mathcal{M}}})$;

6 $\mathbf{v}\leftarrow\mathbf{v}+[v_{j}]$;

7 if
_$\Psi(\mathcal{M}(\mathbf{X}),\mathcal{M}(\mathbf{X}^{\prime}_{\mathbf{v}}))=1$_
then

8 $CPATHS\leftarrow CPATHS\cup\\{\mathbf{v}\\}$ ;

9 break;

10

11 end if

12

13 end for

14

15 end for

Algorithm 1 CPATH: counterfactual paths generation

Subsequently, the observed counterfactual paths are used in Algorithm 2 to
estimate the function $w_{\mathcal{M}}$, i.e., derive the corresponding
adjacency matrix $W$. Specifically, the transition matrix $\mathbf{T}$ is
calculated by iterating through the generated paths. For each path, the length
$l$ of the path is computed. Then, for each consecutive pair of vertices
$(v_{i},v_{i+1})$ in the path, the corresponding entry in $\mathbf{T}$ is
incremented by $k-l+1$. This process accounts for the penalization of longer
paths, as paths with a shorter length contribute more to the edge weights. The
resulting adjacency matrix should capture the underlying dependencies and
interactions between variables based on the observed counterfactual paths.

Data: $CPATHS$, $k\in\mathbf{Z_{+}}$

1 $\mathbf{T}\leftarrow\mathbf{0}_{p\times p}$;

2 for  _$\mathbf{v}$ $\mathrm{\mathbf{in}}$ $CPATHS$_ do

3 $l\leftarrow\mathrm{length}(\mathbf{v})$;

4 if  _$l=1$_ then

5 $\mathbf{T}[v_{1},v_{1}]\leftarrow\mathbf{T}[v_{1},v_{1}]+k-l+1$;

6

7 end if

8 for  _ $i=1$ to $l-1$_ do

9 $\mathbf{T}[v_{i},v_{i+1}]\leftarrow\mathbf{T}[v_{i},v_{i+1}]+k-l+1$;

10

11 end for

12

13 end for

Algorithm 2 CPATH: transition matrix generation

Based on the determined adjacency matrix $W$, we can calculate the importance
of individual variables. A straightforward approach is to compute the fraction
of weights of edges adjacent to the node corresponding to variable
$\mathbf{X}_{j}$. This can be expressed as:

$Imp(\mathbf{X}_{j})=\frac{\sum_{i\in[p]}\mathbf{T}[v_{i},v_{j}]}{\sum_{i\in[p],k\in[p]}\mathbf{T}[v_{i},v_{k}]},$
(1)

where $\mathbf{T}$ is the transition matrix obtained from Algorithm 2 and
$[p]$ denotes the set of indices from $1$ to $p$.

Another way to estimate the importance values is by considering the stationary
distribution of a Markov chain based on transition matrix $\mathbf{P}$ being
the row-wise normalized matrix $\mathbf{T}$. The stationary distribution
represents the long-term probability distribution of being at each node in the
graph and it is a vector $\pi$ that satisfies the condition
$\pi\mathbf{P}=\pi$. In this setting, the importance of variable
$\mathbf{X}_{j}$ is given by the corresponding element $\pi_{j}$ in the
stationary distribution vector.

###  3.2 Counterfactual paths for causal modeling

In the previous section the transition matrix $\mathbf{T}$ was used to derive
feature importances of single features. However, it also contains information
about important interactions of features in form of a weighted directed graph.
In fact, the generated counterfactual paths in $CPATHS$ can be used for causal
modeling and inference, for instance with Bayesian networks. From the absence
and presence of features within the detected counterfactual paths one can
estimate the conditional probability distribution between the important
features, which may help to discovery causal effects. We illustrate this on
two concrete examples in Section 5.

##  4 Evaluation and experimental set-up

####  4.0.1 General approach.

We have generated synthetic data under three conditions: Conditional
dependency, correlation, and conditional independence (see Appendix section
A.1). For each of these datasets two features were simulated fulfilling the
above-mentioned conditions; the rest of the features were added as noise. A
Random Forest was trained on the generated synthetic data and predictions were
made based on the training set. The model was treated as a black box and
explainable methods were exploited to verify the relevant features.

###  4.1 Explaining the model

Our evaluation strategy was to compare the feature importance scores computed
by the explainers (e.g. SHAP (Aas et al., 2021) and LIME (Ribeiro et al.,
2016)) with the model-intern Gini impurity scores generated by the random
forest, which we consider here as ground truth. A high correlation between the
explainers’ feature importance scores and the model’s Gini impurity values
indicates that the explainer can detect the underlying patterns in the model’s
behaviour. Here, we report on the correlation with the model-specific Gini
impurity scores associated with each simulated feature. We compare our
approach with SHAP, LIME, and Permutation Feature Importance (PFI, Fisher et
al., 2019). We report and analyze the performance of PFI, because it relies on
a permutation scheme to derive feature importances, and thus is related to our
method. Unlike CPATH, however, PFI requires the ground-truth labels to
determine feature importance.

###  4.2 Explaining the data

We are aware that the Gini impurity scores themselves are biased due to
several shortcomings (Nembrini et al., 2018). As a consequence, we also
studied the capability of the explainers to detect the most important features
within the data. We could show that the model-specific Gini impurity scores
efficiently reflect the simulated ground truth, which essentially supports the
validity of the first evaluation set-up as described in Section 4.1. In this
particular investigation, however, the focus was more on the simulated
features within the data and not on the features the model actually preferred,
which is addressed by our first evaluation (see Section 4.1).

In the case of SHAP and LIME, feature importance scores for global
explanations were computed as the mean absolute values of the local
explanation scores. For this data-centric experiment we compared our method
also to Conditional Predictive Impact (CPI, Watson and Wright, 2021). CPI
provides variable importance measures taking into account the association
between one or several features and a given outcome.

###  4.3 Evaluation of the explanations’ quality

There are several methods that are used to measure the quality of
explanations; in a survey of surveys (Schwalbe and Finzel, 2023) describe a
taxonomy of Explainable AI methods and mention the fundamental differences and
commonalities between several quality of explanation metrics. Expressing the
expectation of the fact that a substantial change in the model decision’s
logic will also be reflected in its explanations is made by the sensitivity
metric. Furthermore, fidelity (also referred to as faithfulness) is mainly
used in conjecture with a surrogate model (Ancona et al., 2017). It describes
the agreement between the original and surrogate model that is used to provide
the explanations. Typically, models strive to have low infidelity, which is
computed by the following function:

$\textrm{INFD}(\mathrm{\Phi},\mathbf{f},\mathbf{x})=\mathbb{E}_{\mathbf{I}\sim\mu_{\mathbf{I}}}\biggl{[}\mathbf{I}^{T}\mathrm{\Phi}(\mathbf{f},\mathbf{x})-\Bigl{(}\mathbf{f}(\mathbf{x})-\mathbf{f}(\mathbf{x}-\mathbf{I})\Bigl{)}^{2}\biggr{]}$
(2)

Equation 2 provides the infidelity for a model computing the function
$\mathbf{f}$, having $\mathbf{x}$ as input. The explanation method is
represented by $\mathrm{\Phi}$, called “explanation functional”. This metric
is based on the definition of suitable perturbations, tailored to the needs of
the task, the used model, input and the xAI method. They are represented by
the difference $\mathbf{I}=\mathbf{x}-\mathbf{x}_{0}\in\mathbb{R}^{d}$ between
the input $\mathbf{x}$ and $\mathbf{x}_{0}$, where $\mathbf{x}_{0}$ can be a
baseline value, a noisy baseline value or even a Random Variable (RV) (in the
equation 2, $\mathbf{I}\sim\mu_{\mathbf{I}}$ is a RV having the probability
measure $\mu_{\mathbf{I}}$).

Another metric that is related to infidelity to a certain extent, since it is
also based on the principle of using perturbations is sensitivity. In this
case, it is not the input data $\mathbf{x}$ that is perturbed but the input
features; the perturbation strategy is not defined solely by the data
scientist, but it is driven by the relevance/attribution values that each
feature has - according to the xAI method of interest. This metric expresses
an expectation that the computed attribution values for each feature have some
relationship with the computed output of the model. If a feature that is said
to be highly important or decisive for a particular prediction is removed or
replaced by a less informative one, then this action has to have some impact
on the model’s prediction. Even more, this change has to be in some way
analogous to the assigned relevance/attribution of the feature, as computed by
the xAI method. The extension of sensitivity to subsets of features (and not
just one) is called sensitivity-$n$, where $n$ is the cardinality of the
selected feature subset $S$. The corresponding equation is:

$\text{Sens}_{n}=r\Big{(}\sum_{s\in
S_{i}}\mathbf{e}_{s},\mathbf{f}(\mathbf{x})_{c}-\mathbf{f}(\mathbf{x}_{S_{i}})_{c}\Big{)}$
(3)

where $c$ is the predicted class, $\mathbf{x}$ is the original input,
$\mathbf{x}_{S}$ is $\mathbf{x}$ with all features in subset $S$ removed. The
sum of attributions in the subset $S$ is denoted by $\sum_{s\in
S_{i}}\mathbf{e}_{s}$ and its Pearson correlation $r$ with the aforementioned
difference comprises the value of sensitivity-$n$. What is questionable is the
use of this correlation without considering Spearman’s rank correlation
coefficient or Mutual Information (MI) (MacKay, 2003) for capturing non-linear
correlations.

A detailed description of the relationship between sensitivity and infidelity
for different scenarios is described in (Yeh et al., 2019); the researchers
have shown that a small decrease in sensitivity can be achieved when the
explanation is multiplied by a smooth kernel (in this case a Gaussian) without
increasing the infidelity - and under particular circumstances also decreasing
it. This provided a general methodology that can start from any explanation
towards one that has more beneficial quality properties - at least as far as
those metrics are concerned. A detailed analysis of several related metrics
with different representative experiments testing various datasets in the
image processing domain is presented in (Gevaert et al., 2022).

###  4.4 Explanations based on domain knowledge graphs

In the above-described experiments, the counterfactual paths are computed on a
fully-connected graph. The domain expert, however, might have some prior
knowledge about the functional relationships between the studied features
(nodes in the graph), which could be reflected as edges in a knowledge graph.
In these cases, the explanatory factors induced by the counterfactual paths
and their visiting nodes are restricted and guided by domain knowledge.
In this specific evaluation, we compute the network, features and classes in
the following way. First, we generate Barabasi networks (Barabási and Albert,
1999) of varying size and structure. Second, we compute the features
associated with a node using a normal distribution with $N(\mu=0,\sigma=1)$.
In a next step, we randomly select two connected nodes $v_{1}$ and $v_{2}$ and
apply the following formula

$z=5*v_{1}+3*v_{2}.$

We apply the sigmoid function to transform $z$ into a range between 0 and 1

$S(z)=\frac{1}{1+e^{-z}}.$

After we get the transformed values from the sigmoid function, we underlay a
Bernoulli distribution and sample from it to obtain the outcome class.

In the described experiment the feature sampling in CPATH is guided by the
network through random walks. In the remainder of this paper, we refer to it
as CPATHknow.

Finally, we exemplary show the potential of our approaches (CPATH and
CPATHknow) in a bioinformatics application for biomarker discovery. We used
the gene expression data of human breast cancer patient samples for an
experimental evaluation of the herein proposed methodologies. The data was
retrieved from The Cancer Genome Atlas (TCGA) and was preprocessed as
described in (Chereda et al., 2021). We masked the data by the topology of the
Human Protein Reference Database (HPRD) protein-protein interaction (PPI)
network (Keshava Prasad et al., 2009). The resulting dataset comprised 981
patients and 8469 genes. The binary prediction task was to classify the
samples into a group of patients with the luminal A subtype (499 samples) and
patients with other breast cancer subtypes (482 samples). Knowledge-guided
explanations were generated using CPATHknow for the detection of potential
breast cancer-specific biomarkers.

##  5 Results

###  5.1 Synthetic data: correlation with ground-truth

The correlation with the Gini impurities (herein assumed as the model-specific
ground-truth) is depicted in Figure 2. The results indicate that path-based
explanations (CPATH) and permutation-based feature importances (PFI) are more
accurate compared to SHAP and LIME. Especially in the case of conditional
feature dependency they efficiently reflect the model’s internal behaviour
(see Figure 2 & Table 1). The explanations based on LIME are not informative.
The same observation can be made when interpreting the ability of the methods
to detect the relevant features within the data (see Figure 3 and Table 2).
Here, the simulated features are defined as the data-specific ground-truth. In
the case of conditional feature dependency, CPATH is clearly more accurate
than SHAP and LIME. CPATH almost perfectly aligns with the ground truth. CPI
outperforms all methods when the signal-to-noise ratio is low. However, there
is a substantial performance degradation when this ratio increases. It should
be noted that CPI is specially designed for feature selection purposes and not
for explaining the internals of the model. For instance, specialized feature
selection methods often use a random forest, or any other classifier as a
wrapper to detect the most relevant features within the data (Pfeifer et al.,
2022b).

In the case of correlated features, CPATH performs slightly worse than the
competing SHAP method (see Figure 2). Also Figure 3 suggests that SHAP is the
more appropriate explainer in the presence of correlated features. In this
data-focused experiment (Figure 3, Table 2) also LIME is closer to the ground-
truth. Notably, when the signal-to-noise ratio is high all methods perform
almost identically.

Figure 2: Correlation with Gini importance values based on 50 simulations. Table 1: Correlation to Gini importances (min/mean/max) Experiment |  | CPATH | SHAP | LIME | PFI
---|---|---|---|---|---
Conditional dependency (1) |  | 0.8/0.94/1 | 0.4/0.90/1 | 0.2/0.72/1 | 0.8/0.95/1
Conditional dependency (2) |  | 0.8/0.94/1 | 0.2/0.84/1 | 0.2/0.74/1 | 0.6/0.94/1
Correlation |  | 0.4/0.87/1 | 0.4/0.92/1 | 0/0.84/1 | 0.4/0.90/1
Conditional in-dependency |  | 0.6/0.91/1 | 0.6/0.90/1 | 0.6/0.92/1 | 0.6/0.95/1
Table 2: Mean coverage of the important features within the data Experiment | signal/noise | GINI | CPATH | SHAP | LIME | CPI
---|---|---|---|---|---|---
Conditional dependency (1) | 2/2 | 0.990 | 0.960 | 0.855 | 0.720 | 1
Conditional dependency (2) | 2/2 | 0.990 | 0.950 | 0.905 | 0.730 | 0.995
Correlation | 2/2 | 0.995 | 0.960 | 0.905 | 0.725 | 1
Conditional in-dependency | 2/2 | 1 | 1 | 1 | 1 | 1
Conditional dependency (1) | 2/4 | 0.940 | 0.910 | 0.700 | 0.655 | 0.955
Conditional dependency (2) | 2/4 | 0.940 | 0.880 | 0.710 | 0.665 | 0.925
Correlation | 2/4 | 0.760 | 0.700 | 0.705 | 0.770 | 0.805
Conditional in-dependency | 2/4 | 1 | 1 | 1 | 1 | 1
Conditional dependency (1) | 2/6 | 0.875 | 0.830 | 0.635 | 0.595 | 0.860
Conditional dependency (2) | 2/6 | 0.870 | 0.820 | 0.670 | 0.655 | 0.895
Correlation | 2/6 | 0.730 | 0.745 | 0.675 | 0.770 | 0.775
Conditional in-dependency | 2/6 | 1 | 1 | 1 | 1 | 1
Conditional dependency (1) | 2/8 | 0.860 | 0.850 | 0.640 | 0.600 | 0.505
Conditional dependency (2) | 2/8 | 0.845 | 0.805 | 0.615 | 0.600 | 0.400
Correlation | 2/8 | 0.685 | 0.710 | 0.660 | 0.755 | 0.415
Conditional in-dependency | 2/8 | 1 | 1 | 1 | 1 | 0.500
Conditional dependency (1) | 2/10 | 0.815 | 0.715 | 0.615 | 0.600 | 0.520
Conditional dependency (2) | 2/10 | 0.820 | 0.735 | 0.585 | 0.575 | 0.325
Correlation | 2/10 | 0.645 | 0.670 | 0.640 | 0.655 | 0.390
Conditional in-dependency | 2/10 | 0.995 | 0.900 | 0.995 | 0.995 | 0.495
Figure 3: Mean coverage of important features within the simulated data based
on 50 simulations and along the signal-to-noise ratio.

###  5.2 Sensitivity, Infidelity, and Interpretability

We have evaluated the sensitivity and infidelity of the explanations based on
four different datasets, namely Ionosphere, Breast Cancer, Diabetes, and Iris.
The datasets were retrieved using the R package mlbench (Leisch and
Dimitriadou, 2021). We split the data into a train (80%) and test set (20%)
and generated explanations based on the test set. In terms of sensitivity,
CPATH performs well on all datasets and outperforms LIME and SHAP on two out
of four cases (see  Appendix Figure 1). LIME in this experiment is the worst-
performing method. The fidelity of the explanations based on CPATH is not as
good as SHAP, but it is competitive with LIME (see  Appendix Figure 2).

To showcase the interpretability of our proposed approach we analyzed the
Diabetes dataset in more detail (see Figure 4). The accuracy of the random
forest classifier was $AUC=0.97$. A graphical summary of the generated
counterfactual is shown in Figure 4. We see that the glucose variable causes
the largest number of swapped labels when it is used as the first node in a
counterfactual path. More than 20% labels on average swap when this feature is
permuted. Once permuted, insulin and mass increase the swapped fraction up to
30%. However, the largest fraction is obtained by the path starting with
pedigree and further going through glucose, mass, and age. This particular
path caused $>$ 40% of swapped labels on average. These observations suggest
that a combination of the mentioned variables might be a marker for a detailed
medical risk assessment. The feature importances derived from the
counterfactual paths can be obtained from Figure 5. We further learned a
Bayesian network from the inferred counterfactual paths (Figure 5), using the
R-package bnlearn (Scutari, 2017). While the glucose variable was inferred as
the most important feature, the Directed Acyclic Graph in Figure 5b indicates
that it depends on several features, including age, triceps and insulin,
causing the overall importance. From the learned conditional distributions we
see that when these features are present within the given path, it increases
the probability of glucose up to $0.72$ to be also part of that path, causing
the label swap.

Figure 4: Diabetes dataset. Graphical summary of the counterfactual paths.

(a) Feature importance

(b) Conditional dependencies

Figure 5:  Diabetes dataset. Feature importances and conditional dependencies
inferred by Bayesian network learning and here represented as a Directed
Acyclic Graph.

###  5.3 Knowledge-guided explanations and application on Protein-Protein
Interaction Network (PPI)

####  5.3.1 Syntethic Barabasi networks.

The experiments on synthetic Barabasi networks indicate that incorporated
domain knowledge generates more parsimony and interpretable explanations (see
Figure 6). CPATH with incorporated knowledge is more accurate in detecting the
relevant features when the path length is low. CPATH without guided domain
knowledge requires about three times larger paths to capture the ground truth.
Furthermore, knowledge-guided counterfactual explanations converge faster, as
indicated by the higher performance when the number of sampled paths is low
(Figure 6a). Interestingly, CPATH without domain knowledge significantly
outperforms CPATHknow when the size of the sampled paths is high. In that
case, CPATH is able to explore the feature space more efficiently, while the
incorporated knowledge graph seems to hinder optimal convergence. The reason
for that behavior is due to the edge degrees of the graph which leads to
repeated visits of the same node in a random walk. Thus, we can derive the
conclusion that for incorporated knowledge, in the form of a graph, a high
number of generated paths is essential to ensure a sufficient number of
starting nodes so that the whole graph can be explored by the random walks.

(a) Varying number of paths

(b) Varying length of paths

Figure 6:  Mean coverage in detecting the important features within the
simulated data, masked by Barabasi networks. The results are based on 100
simulations, where data and network topology vary in each iteration. (a) The
path length is set to $k=4$. (b) The number of paths is set to $100$.

####  5.3.2 Application on PPI networks.

The retrieved breast cancer gene expression data was split into a train (80%)
and test set (20%). We trained a random forest comprising $1000$ trees and we
evaluated the accuracy based on the independent hold-out test set. The trained
classifier achieved an Area Under the ROC Curve (AUC) of $0.93$. Knowledge-
guided explanations were generated using CPATHknow using a PPI network. We
have generated $10000$ paths using a path length of $k=10$. Overall, $118$
paths of the generated paths were counterfactuals. The feature importances
derived from these paths can be obtained from Figure 7. The top-3 genes were
LAT, RANBP1, and KDM5B. Single-feature test set accuracy for these genes were
$\text{AUC(LAT)}=0.49$, $\text{AUC(RANBP1)}=0.73$, and
$\text{AUC(KDM5B)}=0.54$. These results suggest that LAT and KDM5B on their
own are not sufficient as predictive markers.

LAT is a protein that plays a critical role in the immune system and is
primarily involved in T cell activation and signaling. It is expressed in T
cells and other immune cells, facilitating the transmission of signals from
the T cell receptor to the downstream signalling molecules.

The KDM5B gene, also known as JARID1B or PLU-1, is a gene that encodes a
protein belonging to the family of histone demethylases. Histone demethylases
are enzymes involved in the regulation of gene expression by modifying histone
proteins, which are involved in packaging DNA within the nucleus. Studies have
implicated KDM5B in various biological processes, including development,
differentiation, and cancer. While KDM5B’s role in breast cancer is still
being actively investigated, emerging research suggests its involvement in
breast cancer progression and metastasis (Di Nisio et al., 2023).

The RANBP1 gene, also known as Ran-binding protein 1, encodes a protein
involved in the regulation of the Ran GTPase cycle. The Ran GTPase is
essential for nucleocytoplasmic transport, which controls the movement of
molecules between the nucleus and the cytoplasm. While RANBP1 itself is not
directly associated with breast cancer, aberrant expression or dysregulation
of proteins involved in the Ran GTPase cycle, including RANBP1, have been
implicated in cancer, including breast cancer (Yuen et al., 2016)(Bamodu et
al., 2016). Some studies have suggested that RANBP1 may have tumor suppressor
properties. Decreased expression of RANBP1 has been associated with poor
prognosis and aggressive features in breast cancer patients. Loss or
downregulation of RANBP1 expression may contribute to tumor development and
progression.

In a further investigation, we retrieved the first-order neighbourhood (N) of
the detected and aforementioned genes. The results were AUC(N(LAT))= 0.78,
AUC(N(RANBP1))= 0.82, AUC(N(KDM5B))= 0.67. The neighborhood of the RANBP1 gene
(N(RANBP1)) included four additional genes, namely RAN, CO93, RCC1, and
RANGRF. For this subset, we repeated the generation of CPATH explanations, but
this time without incorporating domain-knowledge. The number of paths was set
to $1000$, and the length of paths was $k=5$. We obtained 864 counterfactual
paths. The generated counterfactual paths can be obtained from Figure 8. The
shortest counterfactual path goes through RANBB1, RCC1 and RAN, leading to an
average of 50% label swaps.

For the first-order neighborhood of RANBP1 we applied bayesian network
learning based on the detected counterfactual paths. The RANBP1 strongly
depends on the CD93 gene. When CD93 is part of a counterfactual path the
probability is $0.74$ that RANBP1 also is included.

(a) Top-10 relevant genes

(b) Conditional dependencies

Figure 7:  Breast cancer dataset. (a) Top-10 relevant genes inferred by CPATH
with incorporated PPI domain knowledge. (b) Conditional dependencies of the
neighborhood features of RANBP1 inferred by Bayesian network learning and here
represented as a Directed Acyclic Graph.  Figure 8: Breast cancer dataset.
Graphical summary of the counterfactual paths.

##  6 Discussion

Apart from the proposed technique, various other approaches can be used to
obtain importance values based on sampling trajectories of finite-length
random walks on the $\widetilde{G_{\mathcal{M}}}$ graph (with corresponding
changes in model predictions). In particular, other approaches can use any
sampled trajectory to compute importance scores without relying directly on
the definition of counterfactual paths using the counterfactual policy $\Psi$.
Therefore, in our initial experiments, we also tested solutions based on
reinforcement learning (RL) methods. The main difference with the technique
presented earlier is that in reinforcement learning-based approaches,
successive sampled trajectories can be used to update the transition matrix
(or directly the importance vector) iteration by iteration, after each episode
of the learning process.

When RL techniques are used to estimate feature importance, the agent learns a
policy that maximizes the cumulative reward over time. The agent explores the
graph by traversing different paths and observing the resulting changes in the
model predictions. Technically, we can formulate the problem as a Markov
decision process (MDP), which is a tuple $\langle S,A,P,R\rangle$, where

  * $\bullet$

$S$ represents the state space corresponding to the set of vertices in the
graph $G_{\mathcal{M}}$ (explanatory variables).

  * $\bullet$

$A$ is the action space, representing the available actions that an RL agent
can take in each state. In this case, the actions correspond to traversing
from one vertex to another in the graph.

  * $\bullet$

$P$ denotes the transition probabilities to other states, given the current
state and the chosen action. In this particular environment, the state to
which an action leads is deterministic (since the action is to select the
appropriate edge of the graph).

  * $\bullet$

$R$ represents the reward function that provides feedback to the RL agent
based on its actions. In the context of feature importance, the reward can be
defined to reflect the changes in model predictions induced by traversing
different paths.

The terminal state is then determined by a fixed maximum number of variables
on the path or is associated with exceeding the threshold of the reward
(significant change in the model predictions).

In our implementation, a solution using the model-free Q-learning algorithm is
available. It estimates the matrix $Q_{p\times p}$ of state-action quality,
from which the importance of variables can be extracted by aggregation,
analogous to the main method. Other available techniques allow the direct
estimation of the state value function (which translates into the importance
of the variables) – this can be done using the temporal difference learning
TD(0) or the first-visit Monte Carlo algorithm. Although these methods give
promising results, they suffer from sensitivity to hyperparameters. They
require further research, investigation and validation, which is beyond the
scope of this paper.

Furthermore, our proposed approach was herein evaluated based on global
explanations. A single test-instance, however, could be explained by first
learning the transition matrix $\mathbf{T}$ on the training set. From this
transition matrix candidate paths could be generated by a Markov process. The
generated path features and their corresponding values could be exchanged with
those in the training set to test for counterfactuals.

##  7 Conclusion

We have developed a novel explainable AI method: counterfactual paths. Unlike
classical feature importance methods, the generated explanations are
efficiently visualized through graphs, which could help detect causal effects
and important interactions between features.

## Code availability

The herein presented method is implemented within the cpath package available
for R and Python at https://github.com/pievos101/cpath.

## Acknowledgements

Parts of this work have been funded by the Austrian Science Fund (FWF),
Project: P-32554 “explainable Artificial Intelligence” (Grantholder AH). Parts
of this work have been funded by the Polish National Science Centre (NCN)
grant 2019/34/E/ST6/00052 (Grantholder PB).

## References

  * Aas et al. (2021) Aas K, Jullum M, Løland A (2021) Explaining individual predictions when features are dependent: More accurate approximations to shapley values. Artificial Intelligence 298:103502
  * Alaa et al. (2021) Alaa AM, Gurdasani D, Harris AL, Rashbass J, van der Schaar M (2021) Machine learning to guide the use of adjuvant therapies for breast cancer. Nature Machine Intelligence 3(8):716–726
  * Ancona et al. (2017) Ancona M, Ceolini E, Öztireli C, Gross M (2017) Towards better understanding of gradient-based attribution methods for deep neural networks. arXiv preprint arXiv:171106104
  * Anguita-Ruiz et al. (2020) Anguita-Ruiz A, Segura-Delgado A, Alcalá R, Aguilera CM, Alcalá-Fdez J (2020) eXplainable Artificial Intelligence (XAI) for the identification of biologically relevant gene expression patterns in longitudinal human studies, insights from obesity research. PLOS Computational Biology 16(4):e1007792
  * Au et al. (2022) Au Q, Herbinger J, Stachl C, Bischl B, Casalicchio G (2022) Grouped feature importance and combined features effect plot. Data Mining and Knowledge Discovery 36(4):1401–1450
  * Bamodu et al. (2016) Bamodu OA, Huang WC, Lee WH, Wu A, Wang LS, Hsiao M, Yeh CT, Chao TY (2016) Aberrant kdm5b expression promotes aggressive breast cancer through malat1 overexpression and downregulation of hsa-mir-448. BMC cancer 16:1–18
  * Baniecki et al. (2023) Baniecki H, Parzych D, Biecek P (2023) The grammar of interactive explanatory model analysis. Data Mining and Knowledge Discovery pp 1–37
  * Barabási and Albert (1999) Barabási AL, Albert R (1999) Emergence of scaling in random networks. science 286(5439):509–512
  * Biecek and Burzykowski (2021) Biecek P, Burzykowski T (2021) Explanatory Model Analysis. Chapman and Hall/CRC, New York
  * Blesch et al. (2023) Blesch K, Watson DS, Wright MN (2023) Conditional feature importance for mixed data. AStA Advances in Statistical Analysis pp 1–20
  * Bücker et al. (2022) Bücker M, Szepannek G, Gosiewska A, Biecek P (2022) Transparency, auditability, and explainability of machine learning models in credit scoring. Journal of the Operational Research Society 73(1):70–90
  * Casalicchio et al. (2019) Casalicchio G, Molnar C, Bischl B (2019) Visualizing the Feature Importance for Black Box Models. In: Machine Learning and Knowledge Discovery in Databases, pp 655–670
  * Chereda et al. (2021) Chereda H, Bleckmann A, Menck K, Perera-Bel J, Stegmaier P, Auer F, Kramer F, Leha A, Beißbarth T (2021) Explaining decisions of graph convolutional neural networks: patient-specific molecular subnetworks responsible for metastasis prediction in breast cancer. Genome medicine 13:1–16
  * Chou et al. (2022) Chou YL, Moreira C, Bruza P, Ouyang C, Jorge J (2022) Counterfactuals and causability in explainable artificial intelligence: Theory, algorithms, and applications. Information Fusion 81:59–83
  * Confalonieri et al. (2021) Confalonieri R, Weyde T, Besold TR, Moscoso del Prado Martín F (2021) Using ontologies to enhance human understandability of global post-hoc explanations of black-box models. Artificial Intelligence 296:103471
  * Covert et al. (2020) Covert I, Lundberg SM, Lee SI (2020) Understanding global feature contributions with additive importance measures. Advances in Neural Information Processing Systems 33:17212–17223
  * Dandl et al. (2020) Dandl S, Molnar C, Binder M, Bischl B (2020) Multi-Objective Counterfactual Explanations. In: Parallel Problem Solving from Nature – PPSN XVI, pp 448–469
  * Di Nisio et al. (2023) Di Nisio E, Licursi V, Mannironi C, Buglioni V, Paiardini A, Robusti G, Noberini R, Bonaldi T, Negri R (2023) A truncated and catalytically inactive isoform of kdm5b histone demethylase accumulates in breast cancer cells and regulates h3k4 tri-methylation and gene expression. Cancer Gene Therapy pp 1–11
  * Fisher et al. (2019) Fisher A, Rudin C, Dominici F (2019) All Models are Wrong, but Many are Useful: Learning a Variable’s Importance by Studying an Entire Class of Prediction Models Simultaneously. Journal of Machine Learning Research 20(177):1–81
  * Gevaert et al. (2022) Gevaert A, Rousseau AJ, Becker T, Valkenborg D, De Bie T, Saeys Y (2022) Evaluating feature attribution methods in the image domain. arXiv preprint arXiv:220212270
  * Gosiewska et al. (2021) Gosiewska A, Kozak A, Biecek P (2021) Simpler is better: Lifting interpretability-performance trade-off via automated feature engineering. Decision Support Systems 150:113556
  * Gozzi et al. (2022) Gozzi N, Malandri L, Mercorio F, Pedrocchi A (2022) XAI for myo-controlled prosthesis: Explaining EMG data for hand gesture classification. Knowledge-Based Systems 240:108053
  * Holzinger et al. (2022) Holzinger A, Saranti A, Molnar C, Biecek P, Samek W (2022) xxAI - Beyond Explainable AI, Springer, chap Explainable AI Methods – A Brief Overview, pp 13–38
  * Karimi et al. (2020) Karimi AH, Barthe G, Balle B, Valera I (2020) Model-agnostic counterfactual explanations for consequential decisions. In: International Conference on Artificial Intelligence and Statistics, pp 895–905
  * Keshava Prasad et al. (2009) Keshava Prasad T, Goel R, Kandasamy K, Keerthikumar S, Kumar S, Mathivanan S, Telikicherla D, Raju R, Shafreen B, Venugopal A, et al. (2009) Human protein reference database—2009 update. Nucleic acids research 37(suppl_1):D767–D772
  * Krzyziński et al. (2023) Krzyziński M, Spytek M, Baniecki H, Biecek P (2023) Survshap(t): Time-dependent explanations of machine learning survival models. Knowledge-Based Systems 262:110234
  * Leisch and Dimitriadou (2021) Leisch F, Dimitriadou E (2021) mlbench: Machine Learning Benchmark Problems. R package version 2.1-3.1
  * MacKay (2003) MacKay DJ (2003) Information theory, inference and learning algorithms. Cambridge university press
  * Molnar et al. (2021) Molnar C, Freiesleben T, König G, Casalicchio G, Wright MN, Bischl B (2021) Relating the Partial Dependence Plot and Permutation Feature Importance to the Data Generating Process. arXiv preprint arXiv:210901433
  * Molnar et al. (2023) Molnar C, König G, Bischl B, Casalicchio G (2023) Model-agnostic feature importance and effects with dependent features: a conditional subgroup approach. Data Mining and Knowledge Discovery pp 1–39
  * Mothilal et al. (2020) Mothilal RK, Sharma A, Tan C (2020) Explaining machine learning classifiers through diverse counterfactual explanations. In: Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency, pp 607–617
  * Nembrini et al. (2018) Nembrini S, König IR, Wright MN (2018) The revival of the gini importance? Bioinformatics 34(21):3711–3718
  * Panigutti et al. (2020) Panigutti C, Perotti A, Pedreschi D (2020) Doctor XAI: An Ontology-Based Approach to Black-Box Sequential Data Classification Explanations. In: Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency, pp 629–639
  * Pawelczyk et al. (2023) Pawelczyk M, Datta T, den Heuvel JV, Kasneci G, Lakkaraju H (2023) Probabilistically Robust Recourse: Navigating the Trade-offs between Costs and Robustness in Algorithmic Recourse. In: The Eleventh International Conference on Learning Representations
  * Pfeifer et al. (2022a) Pfeifer B, Baniecki H, Saranti A, Biecek P, Holzinger A (2022a) Multi-omics disease module detection with an explainable greedy decision forest. Scientific Reports 12(1):16857
  * Pfeifer et al. (2022b) Pfeifer B, Holzinger A, Schimek MG (2022b) Robust random forest-based all-relevant feature ranks for trustworthy ai. Studies in Health Technology and Informatics 294:137–138
  * Ribeiro et al. (2016) Ribeiro MT, Singh S, Guestrin C (2016) “Why Should I Trust You?”: Explaining the Predictions of Any Classifier. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp 1135–1144
  * Saranti et al. (2022) Saranti A, Hudec M, Minarikova E, Takac Z, Großschedl U, Koch C, Pfeifer B, Angerschmid A, Holzinger A (2022) Actionable explainable ai (axai): a practical example with aggregation functions for adaptive classification and textual explanations for interpretable machine learning. Machine Learning and Knowledge Extraction 4(4):924–953
  * Schwalbe and Finzel (2023) Schwalbe G, Finzel B (2023) A comprehensive taxonomy for explainable artificial intelligence: a systematic survey of surveys on methods and concepts. Data Mining and Knowledge Discovery pp 1–59
  * Scutari (2017) Scutari M (2017) Bayesian network constraint-based structure learning algorithms: Parallel and optimized implementations in the bnlearn R package. Journal of Statistical Software 77(2):1–20
  * Wachter et al. (2017) Wachter S, Mittelstadt B, Russell C (2017) Counterfactual explanations without opening the black box: Automated decisions and the GDPR. Harv JL & Tech 31:841
  * Watson and Wright (2021) Watson DS, Wright MN (2021) Testing conditional independence in supervised learning algorithms. Machine Learning 110(8):2107–2129
  * Xu et al. (2022) Xu C, Liu W, Chen Y, Ding X (2022) A supervised case-based reasoning approach for explainable thyroid nodule diagnosis. Knowledge-Based Systems 251:109200
  * Yeh et al. (2019) Yeh CK, Hsieh CY, Suggala A, Inouye DI, Ravikumar PK (2019) On the (in) fidelity and sensitivity of explanations. Advances in Neural Information Processing Systems 32
  * Yuen et al. (2016) Yuen HF, Chan KK, Platt-Higgins A, Dakir EH, Matchett KB, Haggag YA, Jithesh PV, Habib T, Faheem A, Dean FA, et al. (2016) Ran gtpase promotes cancer progression via met receptor-mediated downstream signaling. Oncotarget 7(46):75854

##  Appendix 0.A Appendix

###  0.A.1 Simulation Algorithms

1 $n=100,p=4,D^{n\times p}\sim N(0,2),y^{n\times 1}\sim Bern(0.5)$;

2 $a=0,b=0$;

3 for  _$i=1$ to $n$ _ do

4

5 if  _$D[i,1]\geq a$_ then

6

7 if  _$D[i,2]\geq b$_ then

8 $y[i]=1$;

9

10 end if

11

12 else

13 if  _$D[i,2]\geq b$_ then

14 $y[i]=0$;

15

16 end if

17

18 end if

19

20 end for

Algorithm 3 Simulation 1. Conditional dependency (1)

1 $n=100,p=4,D^{n\times p}\sim N(0,2),y^{n\times 1}\sim Bern(0.5)$;

2 $a=0,b=0$;

3 for  _$i=1$ to $n$ _ do

4

5 if  _$D[i,1]\geq a$_ then

6

7 if  _$D[i,2]\leq b$_ then

8 $y[i]=1$;

9

10 else

11 $y[i]=0$;

12

13 end if

14

15 end if

16

17 end for

Algorithm 4 Simulation 2. Conditional dependency (2)

1 $n=100,p=4,D^{n\times p}\sim N(0,2),y^{n\times 1}\sim Bern(0.5)$;

2 $a=0$;

3

4for  _$i=1$ to $n$ _ do

5

6 if  _$D[i,1]\geq a\And D[i,2]\geq a$_ then

7 $y[i]=1$;

8

9 end if

10

11 end for

Algorithm 5 Simulation 3. Correlation

1 $n=100,p=4,D^{n\times p}\sim N(0,2),y^{n\times 1}\sim Bern(0.5)$;

2 $a=0$;

3 for  _$i=1$ to $n$ _ do

4

5 if  _$D[i,1]\geq a\And D[i,2]\geq a$_ then

6 $next$;

7

8 end if

9

10 if  _$D[i,1]\geq a$_ then

11 $y[i]=1$;

12

13 end if

14

15 if  _$D[i,2]\geq a$_ then

16 $y[i]=0$;

17

18 end if

19

20 end for

21

Algorithm 6 Simulation 4. Conditional in-dependency

###  0.A.2 Supporting Figures

(a) Ionosphere

(b) Breast cancer

(c) Diabetes

(d) Iris

Figure Appendix Figure 1: Sensitivity

(a) Ionosphere

(b) Breast cancer

(c) Diabetes

(d) Iris

Figure Appendix Figure 2: Infidelity