# Explaining Predictive Models with Mixed Features
Using Shapley Values and Conditional Inference Trees

Annabelle Redelmeier, Martin Jullum, Kjersti Aas

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

## To cite this version:

Annabelle Redelmeier, Martin Jullum, Kjersti Aas. Explaining Predictive Models with Mixed Features Using Shapley Values and Conditional Inference Trees. 4th International Cross-Domain Conference for Machine Learning and Knowledge Extraction (CD-MAKE), Aug 2020, Dublin, Ireland. pp.117-137, 10.1007/978-3-030-57321-8_7. hal-03414718

HAL Id: hal-03414718

https://inria.hal.science/hal-03414718

Submitted on 4 Nov 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Explaining predictive models with mixed features using Shapley values and conditional inference trees

Annabelle Redelmeier, Martin Jullum, and Kjersti Aas

Norwegian Computing Center, P.O. Box 114, Blindern, N-0314 Oslo, Norway
{annabelle.redelmeier,jullum,kjersti}@nr.no

**Abstract.** It is becoming increasingly important to explain complex, black-box machine learning models. Although there is an expanding literature on this topic, Shapley values stand out as a sound method to explain predictions from any type of machine learning model. The original development of Shapley values for prediction explanation relied on the assumption that the features being described were independent. This methodology was then extended to explain dependent features with an underlying continuous distribution. In this paper, we propose a method to explain mixed (i.e. continuous, discrete, ordinal, and categorical) dependent features by modeling the dependence structure of the features using conditional inference trees. We demonstrate our proposed method against the current industry standards in various simulation studies and find that our method often outperforms the other approaches. Finally, we apply our method to a real financial data set used in the 2018 FICO Explainable Machine Learning Challenge and show how our explanations compare to the FICO challenge Recognition Award winning team.

**Keywords:** Explainable AI · Shapley values · conditional inference trees · feature dependence · prediction explanation

## 1 Introduction

Due to the ongoing data and artificial intelligence (AI) revolution, an increasing number of crucial decisions are being made with complex automated systems. It is therefore becoming ever more important to understand how these systems make decisions. Such systems often consist of ‘black-box’ machine learning models which are trained to predict an outcome/decision based on various input data (i.e. features). Consider, for instance, a model that predicts the price of car insurance based on the features age and gender of the individual, type of car, time since the car was registered, and number of accidents in the last five years. For such a system to work in practice, the expert making the model, the insurance brokers communicating the model, and the policyholders vetting the model should know which features drive the price of insurance up or down.

Although there are numerous ways to explain complex models, one way is to show how the individual features contribute to the overall predicted value
for a given individual<fnref n="1" />. The Shapley value framework is recent methodology to calculate these contributions [15, 20, 21]. In the framework, a Shapley value is derived for each feature given a prediction (or ‘black-box’) model and the set of feature values for the given individual. The methodology is such that the sum of the Shapley values for the individual equals their prediction value so that the features with the largest (absolute) Shapley values are the most important.

The Shapley value concept is based on economic game theory. The original setting is as follows: Imagine a game where *N* players cooperate in order to maximize the *total gains* of the game. Suppose further that each player is to be given a certain payout for his or her efforts. Lloyd Shapley [19] discovered a way to distribute the total gains of the game among the players according to certain desirable axioms. For example, players that do not contribute anything get a payout of 0; two players that contribute the same regardless of other players get the same payout; and the sum of the payouts equals the total gains of the game. A player's payout is known as his or her *Shapley value*.

[15, 20, 21] translate Shapley values from the game theory setting to a machine learning setting. The cooperative game becomes the individual, the total gains of the game become the prediction value, and the players become the feature values. Then, analogous to game theory, the Shapley value of one of the features (called the *Shapley value explanation*) is how the feature contributes to the overall prediction value.

Figure 1 shows how such Shapley value explanations can be visualized for two examples of the aforementioned car insurance scenario. For the individual on the left, ‘number of accidents’ pulls the predicted insurance price up (its Shapley value is positive) whereas ‘gender’ and ‘age’ pull it down. The features ‘type of car’ and ‘time since registration’ only minimally affect the prediction. For the individual on the right, ‘gender’ and ‘type of car’ pull the predicted insurance price up whereas ‘age’, ‘time since registration’, and ‘number of accidents’ pull it marginally down. The sum of the Shapley values of each individual gives the predicted price of insurance (123.5 and 229.9 USD/month, respectively). Note that ‘none’ is a fixed average prediction contribution not due to any of the features in the model.

As we demonstrate in Section 2, calculating Shapley values is not necessarily straightforward or computationally simple. To simplify the estimation problem, [15, 20, 21] assume the features are independent. However, [1] shows that this may lead to severely inaccurate Shapley value estimates and, therefore, incorrect explanations when the features are not independent. [1] extends [15]’s methodology to handle dependent *continuously* distributed features by modeling/estimating the dependence between the features. However, as exemplified by the aforementioned car insurance example, practical modeling scenarios often involve a mixture of different feature types: continuous (age and time since the car was registered), discrete (number of accidents in the last five years), and categorical (gender and type of car). Thus, there is a clear need to extend the

<footnote n="1">¹ Here, ‘individual’ could be an individual person or an individual non-training observation - not necessarily a person.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** An example of using Shapley values to show how the predicted price of car insurance can be broken down into the respective features.</figcaption></figure>

Shapley value methodology to handle dependent mixed (i.e. continuous, discrete, ordinal, categorical) features.

While it is, in principle, possible to naively apply some of the methods proposed by [1] to discrete or categorical features, it is unlikely that they will function well. This will typically require encoding categorical features with *L* different categories into *L* − 1 new indicator features using one-hot encoding. The main drawback of this approach is that the feature dimension increases substantially unless there are very few categories. Computational power is already a non-trivial issue with the Shapley value framework (see [1]), so this is not a feasible approach unless the number of categories or features is very small.

The aim of this paper is to show how we can extend the Shapley value frame-
work to handle mixed features without assuming the features are independent.
We propose to use a special type of tree model, known as conditional inference
trees [12], to model the dependence between the features. This is similar to [1]'s
extension of [15]'s work but for mixed features. We use tree models since they
are inherently good at modeling both simple and complex dependence structures
in mixed data types [8]. The conditional inference tree model has the additional
advantage of naturally extending to multivariate responses, which is required
in this setting. Since conditional inference trees handle categorical data, this
approach does not require one-hot encoding any features resulting in a much
shorter computation time. In addition, we do not run the risk of estimating
one-hot encoded features using a method not designed for the sort.

The rest of the paper is organized as follows. We begin by explaining the fun-
damentals of the Shapley value framework in an explanation setting in Section 2
and then outline how to extend the method to mixed features using conditional
inference trees in Section 3. In Section 4, we present various simulation studies
for both continuous and categorical features that demonstrate that our method
works in a variety of settings. Finally, in Section 5, we apply our method to the
2018 FICO Explainable Machine Learning Challenge data set and show how the estimated Shapley values differ when calculated using various feature distribution assumptions. We also compare the feature importance rankings calculated using Shapley values with the rankings calculated by the 2018 FICO challenge Recognition Award winning team. In Section 6, we conclude.

The Shapley methodology from [1] is implemented in the software package `shapr` [18] in the R programming language [17]. Our new approach is implemented as an extension to the `shapr` package. To construct the conditional inference trees in R, we use the packages `party` and `partykit` [12, 13].

# 2 Shapley values

## 2.1 Shapley values in game theory

Suppose we are in a cooperative game setting with $M$ players, $j = 1, \dots, M$, trying to maximize a payoff. Let $\mathcal{M}$ be the set of all players and $\mathcal{S}$ any subset of $\mathcal{M}$. Then the Shapley value [19] for the $j$th player is defined as

$$ \phi_j = \sum_{S \subseteq \mathcal{M} \setminus \{j\}} \frac{|S|!(M - |S| - 1)!}{M!} (v(S \cup \{j\}) - v(S)). \quad (1) $$

$v(\mathcal{S})$ is the contribution function which maps subsets of players to real numbers representing the worth or contribution of the group $\mathcal{S}$ and $|\mathcal{S}|$ is the number of players in subset $\mathcal{S}$.

In the game theory sense, each player receives $\phi_j$ as their payout. From the formula, we see that this payout is just a weighted sum of the player's marginal contributions to each group $\mathcal{S}$. Lloyd Shapley [19] proved that distributing the total gains of the game in this way is 'fair' in the sense that it obeys certain important axioms.

## 2.2 Shapley values for explainability

In a machine learning setting, imagine a scenario where we fit $M$ features, $\boldsymbol{x} = (x_1, \dots, x_M)$, to a univariate response $y$ with the model $f(\boldsymbol{x})$ and want to explain the prediction $f(\boldsymbol{x})$ for a specific feature vector $\boldsymbol{x} = \boldsymbol{x}^*$. [15, 20, 21] suggest doing this with Shapley values where the predictive model replaces the cooperative game and the features replace the players. To use (1), [15] defines the contribution function $v(\mathcal{S})$ as the following expected prediction

$$ v(\mathcal{S}) = \mathbb{E}[f(\boldsymbol{x})|\boldsymbol{x}_{\mathcal{S}} = \boldsymbol{x}_{\mathcal{S}}^*]. \quad (2) $$

Here $\boldsymbol{x}_{\mathcal{S}}$ denotes the features in subset $\mathcal{S}$ and $\boldsymbol{x}_{\mathcal{S}}^*$ is the subset $\mathcal{S}$ of the feature vector $\boldsymbol{x}^*$ that we want to explain. Thus, $v(\mathcal{S})$ denotes the expected prediction given that the features in subset $\mathcal{S}$ take the value $\boldsymbol{x}_{\mathcal{S}}^*$.

Calculating the Shapley value for a given feature $x_j$ thus becomes the arduous task of computing (1) but replacing $v(\mathcal{S})$ with the conditional expectation (2).
It is clear that the sum in (1) grows exponentially as the number of features, $M$, increases. [15] cleverly approximates this weighted sum in a method they call Kernel SHAP. Specifically, they define the Shapley values as the optimal solution to a certain weighted least squares problem. They prove that the Shapley values can explicitly be written as

$$ \phi = (Z^T W Z)^{-1} Z^T W v, \qquad (3) $$

where $Z$ is the $2^M \times (M+1)$ binary matrix representing all possible combinations of the $M$ features, $W$ is the $2^M \times 2^M$ diagonal matrix containing Shapley weights, and $v$ is the vector containing $v(S)$ for every $S$. The full derivation is described in [1].

To calculate (3), we still need to compute the contribution function $v(S)$ for different subsets of features, $S$. When the features are continuous, we can write the conditional expectation (2) as

$$ \mathbb{E}[f(\mathbf{x})|\mathbf{x}_S = \mathbf{x}_S^*] = \mathbb{E}[f(\mathbf{x}_{\bar{S}}, \mathbf{x}_S)|\mathbf{x}_S = \mathbf{x}_S^*] = \int f(\mathbf{x}_{\bar{S}}, \mathbf{x}_S^*)p(\mathbf{x}_{\bar{S}}|\mathbf{x}_S = \mathbf{x}_S^*) d\mathbf{x}_S^*, \qquad (4) $$

where $\mathbf{x}_{\bar{S}}$ is the vector of features not in $S$ and $p(\mathbf{x}_{\bar{S}}|\mathbf{x}_S = \mathbf{x}_S^*)$ is the conditional distribution of $\mathbf{x}_{\bar{S}}$ given $\mathbf{x}_S = \mathbf{x}_S^*$. Note that in the rest of the paper we use $p(\cdot)$ to refer to both probability mass functions and density functions (made clear by the context). We also use lower case $x$-s for both random variables and realizations to keep the notation concise.

Since the conditional probability function is rarely known, [15] replaces it with the simple (unconditional) probability function

$$ p(\mathbf{x}_{\bar{S}}|\mathbf{x}_S = \mathbf{x}_S^*) = p(\mathbf{x}_{\bar{S}}). \qquad (5) $$

The integral then becomes

$$ \mathbb{E}[f(\mathbf{x})|\mathbf{x}_S = \mathbf{x}_S^*] = \int f(\mathbf{x}_{\bar{S}}, \mathbf{x}_S^*)p(\mathbf{x}_{\bar{S}}) d\mathbf{x}_{\bar{S}}^*, \qquad (6) $$

which is estimated by randomly drawing $K$ times from the full training data set and calculating

$$ v_{\text{KerSHAP}}(S) = \frac{1}{K} \sum_{k=1}^{K} f(\mathbf{x}_{\bar{S}}^k, \mathbf{x}_S^*), \qquad (7) $$

where $\mathbf{x}_{\bar{S}}^k$, $k = 1, \dots, K$ are the samples from the training set and $f(\cdot)$ is the estimated prediction model.

Unfortunately, when the features are not independent, [1] demonstrates that naively replacing the conditional probability function with the unconditional one leads to very inaccurate Shapley values. [1] then proposes multiple methods for estimating $p(\mathbf{x}_S|\mathbf{x}_S = \mathbf{x}_S^*)$ without relying on the naive assumption in (5). However, these methods are only constructed for continuous features. In the next section we demonstrate how we can use conditional inference trees to extend the current Shapley framework to handle mixed features.
## 3 Extending the Shapley framework with conditional inference trees

Conditional inference trees (ctree) [12] is a type of recursive partitioning algorithm like CART (classification and regression trees) [4] and C4.5 [16]. Just like these algorithms, ctree builds trees recursively, making binary splits on the feature space until a given stopping criterion is fulfilled. The difference between ctree and CART/C4.5 is how the feature/split point and stopping criterion are chosen. CART and C4.5 solve for the feature and split point simultaneously: each feature and split point is tried together and the best pair is the combination that results in the smallest error (often based on the squared error loss or binary cross-entropy loss depending on the response). Ctree, on the other hand, proceeds sequentially: the splitting feature is chosen using statistical significance tests and then the split point is chosen using any type of splitting criterion [12]. According to [12], choosing the splitting feature without first checking for the potential split points avoids being biased towards features with many split points. In addition, unlike CART and C4.5, ctree is defined independently of the dimension of the response variable. This is advantageous since proper handling of multivariate responses is crucial for our problem.

### 3.1 Conditional inference tree algorithm

Suppose that we have a training data set with $p$ features, a $q$ dimensional response, and $n$ observations: $\{y_i, x_i\}_{i=1,...,n}$ with $y_i = (y_{i1},...,y_{iq})$ and $x_i = (x_{i1},...,x_{ip})$. Suppose further that the responses come from a sample space $\mathcal{Y} = \mathcal{Y}_1 \times \cdots \times \mathcal{Y}_q$ and the features come from a sample space $\mathcal{X} = \mathcal{X}_1 \times \cdots \times \mathcal{X}_p$. Then conditional inference trees are built using the following algorithm:

1. For a given node in the tree, test the global null hypothesis of independence between all of the $p$ features and the response $y$. If the global hypothesis cannot be rejected, do not split the node. Otherwise, select the feature $x_j$ that is the least independent of $y$.

2. Choose a splitting point in order to split $\mathcal{X}_j$ into two disjoint groups.

Steps 1 and 2 are repeated until no nodes are split.

The global null hypothesis can be written as

$$H_0: \cap_{j=1}^{p} H_0^j,$$

where the $p$ partial hypotheses are

$$H_0^j: F(\mathbf{Y}|X_j) = F(\mathbf{Y}),$$

and $F(\cdot)$ is the distribution of $\mathbf{Y}$.

Specifically, we calculate the $p$ $P$-values for the partial hypotheses and combine them to form the global null hypothesis $P$-value. If the $P$-value for the global null hypothesis is smaller than some predetermined level $\alpha$, we reject the
global null hypothesis and assume that there is some dependence between the features and the response. The feature that is the least independent of the response (i.e. has the smallest partial P-value) becomes the splitting feature. If the global null hypothesis is not rejected, we do not split the node. The size of the tree is controlled using the parameter α. As α increases, we are more likely to reject the global null hypothesis and therefore split the node. This results in deeper trees. However, if α is too large, we risk that the tree overfits the data.

Step 2 can be done using any type of splitting criterion, specifically it can be done with the permutation test framework devised by [12]. Note that this method is not tied to a specific feature type and can be used with mixtures of continuous, discrete, ordinal, and categorical features. We refer the reader to the original paper [12] for more details on how to form the test statistic, associated distribution, and P-values.

## 3.2 Extending the Shapley value framework with conditional inference trees

As already mentioned, one of the main limitations with the Shapley value framework is estimating the contribution function (2) when the conditional distribution of the features is unknown but the features are assumed dependent. [1] estimates (2) by modeling the conditional probability density function

$$p(x_S|\mathbf{x}_S = \mathbf{x}_S^*) \qquad (8)$$

using various approaches. Then, [1] samples K times from this modeled conditional distribution function and uses these samples to estimate the integral (4) using (7).

We extend this approach to mixed features by modeling the conditional distribution function (8) using conditional inference trees. We fit a tree to our training data where the features are $\mathbf{x}_S$ and the response is $\mathbf{x}_S$ with the algorithm described in Section 3.1. Then for a given $\mathbf{x}_S^*$, we find its leaf in the tree and sample K times from the $\mathbf{x}_S$ part of the training observations in that node to obtain $\mathbf{x}_S^k$, $k = 1, \dots, K$. Finally, we use these samples to estimate (4) using the approximation (7).

We fit a new tree to every combination of features $\mathbf{x}_S$ and response $\mathbf{x}_S$. Once $v(S)$ is estimated for every $S$, we follow [14]'s steps and estimate the Shapley value of this feature with (3). Since conditional inference trees handle continuous, discrete, ordinal, and categorical features; univariate and multivariate responses; and any type of dependence structure, using conditional inference trees to estimate (8) is a natural extension to [1]'s work. Below, we use the term ctree to refer to estimating Shapley value explanations using conditional inference trees.

## 4 Simulation studies

In this section, we discuss two simulation studies designed to compare different ways to estimate Shapley values. Specifically, we compare our ctree estimation
approach with [15]'s independence estimation approach (below called *independence*) and [1]'s empirical and Gaussian estimation approaches. A short description of each approach is in Table 1.

<table><caption>**Table 1.** A short description of the approaches used to estimate (8) in the simulation studies.</caption><thead><tr><th>Method</th><th>Citation</th><th>Description</th></tr></thead><tbody><tr><td>independence</td><td>[14]</td><td>Assume the features are independent. Assume (2) is (6) and estimate it with (7) where $x_S^k$ are sub-samples from the training data set.</td></tr><tr><td>empirical</td><td>[1]</td><td>Calculate the distance between the set of features being explained and every training instance. Use this distance to calculate a weight for each training instance. Approximate (2) using a function of these weights.</td></tr><tr><td>Gaussian (100)</td><td>[1]</td><td>Assume the features are jointly Gaussian. Estimate the mean/covariance of this conditional distribution and then sample 100 times from this distribution. Estimate (2) with (7) using this sample.</td></tr><tr><td>Gaussian (1000)</td><td>[1]</td><td>The same as Gaussian (100), but we sample 1000 times.</td></tr><tr><td>ctree</td><td></td><td>See Section 3.2. Set $\alpha = 0.5$.</td></tr><tr><td>ctree-onehot</td><td></td><td>Convert the categorical features into one-hot encoded features and then apply the algorithm in Section 3.2 to these binary features. This approach is used only as a reference.</td></tr></tbody></table>

The independence, empirical, and Gaussian approaches are all implemented in the R package `shapr` [18]. We implement the ctree method in the `shapr` package as an additional method. Building the conditional inference trees for each combination of features is done using either the `party` package or `partykit` package in R [12, 13]. Although `party` is faster than `partykit`, it sometimes runs into a convergence error related to the underlying linear algebra library in R (error code 1 from Lapack routine 'dgesdd'). We therefore fall back to `partykit` when this error occurs. Both packages typically give identical results.

In the first simulation study, we simulate only categorical features and in the second, we simulate both categorical and continuous features. Then, we estimate the Shapley values of each test observation with the methods in Table 1 and compare them against the truth using a mean absolute error type of performance measure.

For simplicity, in both situations we restrict ourselves to a linear predictive function of the form

$$f(\boldsymbol{x}) = \alpha + \sum_{\{j:j \in C_{\text{cat}}\}} \sum_{l=2}^{L} \beta_{jl} \mathbf{1}(x_j = l) + \sum_{\{j:j \in C_{\text{cont}}\}} \gamma_j x_j, \quad (9)$$
where $C_{cat}$ and $C_{cont}$ denote, respectively, the set of categorical and continuous features, $L$ is the number of categories for each of the categorical features, and $1(x_j = l)$ is the indicator function taking the value 1 if $x_j = l$ and 0 otherwise. $\alpha, \beta_{jl}$ for $j \in C_{cat}$, $l = 2, \dots, L$, and $\gamma_j$ for $j \in C_{cont}$ are the parameters in the linear model. We define $M = |C_{cat}| + |C_{cont}|$, where $|C_{cat}|$ and $|C_{cont}|$ denote the number of categorical and continuous features, respectively.

The empirical and Gaussian methods cannot handle categorical features. For these methods, we transform the categorical features into one-hot encoded features. If the categorical feature originally has $L$ categories, the one-hot encoded transformation creates $L-1$ binary features representing the second, third (etc) categories. The first category is represented by the intercept. The Shapley value of each categorical feature is then the sum of the Shapley values of the corresponding one-hot encoded features.

## 4.1 Evaluation method

We measure the performance of each method based on the mean absolute error (MAE), across both the features and sample space. This is defined as

$$ \text{MAE}(\text{method } q) = \frac{1}{M} \sum_{j=1}^{M} \sum_{i=1}^{T} p(\boldsymbol{x}_i) |\phi_{j,\text{true}}(\boldsymbol{x}_i) - \phi_{j,q}(\boldsymbol{x}_i)|, \quad (10) $$

where $\phi_{j,q}(\boldsymbol{x})$ and $\phi_{j,\text{true}}(\boldsymbol{x})$ denote, respectively, the Shapley value estimated with method $q$ and the corresponding true Shapley value for the prediction $f(\boldsymbol{x})$. In addition, $M$ is the number of features and $T$ is the number of test observations. For the case with only categorical features, the set $\{\boldsymbol{x}_i : i = 1, \dots, T\}$ corresponds to all the unique combinations of features and $p(\boldsymbol{x}_i)$ is the probability mass function of $\boldsymbol{x}$ evaluated at $\boldsymbol{x}_i^2$. In the case where we have both categorical and numerical features, the set $\{\boldsymbol{x}_i : i = 1, \dots, T\}$ is sampled from the distribution of $\boldsymbol{x}$ and $p(\boldsymbol{x}_i)$ is set to $1/T$ for all $i$.

## 4.2 Simulating dependent categorical features

To simulate $M$ dependent categorical features with $L$ categories each, we first simulate an $M$-dimensional Gaussian random variable with a specified mean $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$

$$ (\tilde{x}_1, \dots, \tilde{x}_M) \sim N_M(\boldsymbol{\mu}, \boldsymbol{\Sigma}). \quad (11) $$

We then transform each feature, $\tilde{x}_j$, into a categorical feature, $x_j$, using the following transformation:

$$ x_j = l, \text{ if } v_l < \tilde{x}_j \le v_{l+1}, \text{ for } l = 1, \dots, L \text{ and } j = 1, \dots, M, \quad (12) $$

<footnote n="2">² If there are many categorical features or number of categories, we instead use a subset of the most likely combinations and scale the probabilities such that they sum to 1 over those combinations.</footnote>
where $v_1, \dots, v_{L+1}$ is an increasing, ordered set of cut-off values defining the categories with $v_1 = -\infty$ and $v_{L+1} = +\infty$. We redo this $n_{train}$ times to create a training data set of $M$ dependent categorical features. The strength of the dependencies between the categorical features is controlled by the correlations specified in $\Sigma$. Note that the actual value of $x_j$ is irrelevant – the features are treated as non-ordered categorical features.

For the simulation setting in Section 4.5 where there are both categorical and continuous features, we first sample $M = |C_{cat}| + |C_{cont}|$ features using (11). Then we transform the features $\tilde{x}_j$ where $j \in C_{cat}$ to categorical ones using (12), and leave the remaining features untouched (i.e letting $x_j = \tilde{x}_j$, when $j \in C_{cont}$). This imposes dependence both within and between all feature types.

## 4.3 Calculating the true Shapley values

To evaluate the performance of the different methods with the MAE from Section 4.1, we need to calculate the true Shapley values, $\phi_{j,\text{true}}(x^*)$, $j = 1, \dots, M$, for all feature vectors where $x^* = x_i$, $i = 1, \dots, T$. This requires the true conditional expectation (2) for all feature subsets $S$. We compute these expectations differently depending on whether the features are all categorical or whether there are both categorical and continuous features. The linearity of the predictive function (9) helps to simplify the computations for the latter case. Since there is no need of it in the former case, we present that case more generally.

When the features are all categorical, the desired conditional expectation can be written as

$$ \mathbb{E}[f(\boldsymbol{x})|\boldsymbol{x}_S = \boldsymbol{x}_S^*] = \sum_{\boldsymbol{x}_S \in \mathcal{X}_S} f(\boldsymbol{x}_S^*, \boldsymbol{x}_{\bar{S}})p(\boldsymbol{x}_{\bar{S}}|\boldsymbol{x}_S = \boldsymbol{x}_S^*), $$

where $\mathcal{X}_{\bar{S}}$ denotes the feature space of the feature vector $\boldsymbol{x}_{\bar{S}}$ which contains $|\bar{S}|^L$ unique feature combinations. Thus, all we need is the conditional probability $p(\boldsymbol{x}_{\bar{S}}|\boldsymbol{x}_S = \boldsymbol{x}_S^*)$ for each combination of $\boldsymbol{x}_{\bar{S}} \in \mathcal{X}_{\bar{S}}$. Using standard probability theory, this conditional probability can be written as

$$ p(\boldsymbol{x}_{\bar{S}}|\boldsymbol{x}_S) = \frac{p(\boldsymbol{x}_{\bar{S}}, \boldsymbol{x}_S)}{p(\boldsymbol{x}_S)}, $$

and then evaluated at the desired $\boldsymbol{x}_S^*$. Since all feature combinations correspond to hyperrectangular subspaces of Gaussian features, we can compute all joint probabilities exactly using the cut-offs $v_1, \dots, v_{L+1}$:

$$ p(x_1 = l_1, \dots, x_M = l_M) = P(v_{l_1} < \tilde{x}_1 \le v_{l_1+1}, \dots, v_{l_M} < \tilde{x}_M \le v_{l_M+1}), $$

for $l_j = 1, \dots, L$, $j = 1, \dots, M$. Here $p(\cdot)$ denotes the joint probability mass function of $\boldsymbol{x}$ while $P(\cdot)$ denotes the joint continuous distribution function of $\tilde{x}$. The probability on the right is easy to compute based on the cumulative distribution function of the multivariate Gaussian distribution (we used the R package mvtnorm [11]). The marginal and joint probability functions based on
only a subset of the features are computed analogously based on a subset of the full Gaussian distribution, which is also Gaussian.

For the situation where some of the features are categorical and some are continuous, the computation of the conditional expectation is more arduous. However, due to the linearity of the predictive function (9), the conditional expectation reduces to a linear combination of two types of univariate expectations. Let $S_{cat}$ and $\tilde{S}_{cat}$ refer to, respectively, the $S$ and $\tilde{S}$ part of the categorical features, $C_{cat}$, with analogous sets $S_{cont}$ and $\tilde{S}_{cont}$ for the continuous features. We then write the desired conditional expectation as

$$ 
\begin{aligned}
\mathbb{E}[f(\mathbf{x})|\mathbf{x}_S = \mathbf{x}_S^*] &= \mathbb{E}\left[\alpha + \sum_{j \in C_{\text{cat}}} \sum_{l=2}^{L} \beta_{jl} \mathbf{1}(x_j = l) + \sum_{j \in C_{\text{cont}}} \gamma_k x_j \middle| \mathbf{x}_S = \mathbf{x}_S^*\right] \\
&= \alpha + \sum_{j \in C_{\text{cat}}} \sum_{l=2}^{L} \beta_{jl} \mathbb{E}[\mathbf{1}(x_j = l)|\mathbf{x}_S = \mathbf{x}_S^*] + \sum_{j \in C_{\text{cont}}} \gamma_j \mathbb{E}[x_j|\mathbf{x}_S = \mathbf{x}_S^*] \\
&= \alpha + \sum_{j \in \tilde{S}_{\text{cat}}} \sum_{l=2}^{L} \beta_{jl} \mathbb{E}[\mathbf{1}(x_j = l)|\mathbf{x}_S = \mathbf{x}_S^*] + \sum_{j \in \tilde{S}_{\text{cat}}} \sum_{l=2}^{L} \beta_{jl} \mathbf{1}(x_j^* = l) \\
&\quad + \sum_{j \in \tilde{S}_{\text{cont}}} \gamma_j \mathbb{E}[x_j|\mathbf{x}_S = \mathbf{x}_S^*] + \sum_{j \in \tilde{S}_{\text{cont}}} \gamma_j x_j^*.
\end{aligned}
$$

Then, we just need expressions for the two conditional expectations: $\mathbb{E}[\mathbf{1}(x_j = l)|\mathbf{x}_S = \mathbf{x}_S^*]$ for $j \in \tilde{S}_{cat}$ and $\mathbb{E}[x_j|\mathbf{x}_S = \mathbf{x}_S^*]$ for $j \in \tilde{S}_{cont}$. To calculate them, we use results from [3] on selection (Gaussian) distributions in addition to basic probability theory and numerical integration. Specifically, the conditional expectation for the continuous features takes the form

$$ \mathbb{E}[x_j|\mathbf{x}_S] = \int_{-\infty}^{\infty} xg(x) \frac{p(\mathbf{x}_S|x_j=x)}{p(\mathbf{x}_S)} dx, \quad (13) $$

where $g(x)$ denotes the density of the standard normal (Gaussian) distribution, $p(\mathbf{x}_S|x_j=x)$ is the conditional distribution of $\mathbf{x}_S$ given $x_j$, and $p(\mathbf{x}_S)$ is the marginal distribution of $\mathbf{x}_S$. The latter two are both Gaussian and can be evaluated at the specific vector $\mathbf{x}_S^*$ using the R package mvtnorm [11]. Finally, the integral is solved using numerical integration.

For the second expectation, recall that $x_j = l$ corresponds to the original Gaussian variable $\tilde{x}_j$ falling in the interval $(v_l, v_{l+1}]$. Then, the conditional expectation for the categorical features takes form

$$ 
\begin{aligned}
\mathbb{E}[\mathbf{1}(x_j = l)|\mathbf{x}_S] &= P(v_l < \tilde{x}_j \le v_{l+1}|\mathbf{x}_S) \\
&= \int_{v_l}^{v_{l+1}} g(x) \frac{p(\mathbf{x}_S|\tilde{x}_j = x)}{p(\mathbf{x}_S)} dx,
\end{aligned}
$$

which can be evaluated similarly to (13) and solved with numerical integration. Once we have computed the necessary conditional expectations for each of the
$2^M$ feature subsets $\mathcal{S}$, we compute the Shapley values using (3). This goes for both the pure categorical case and the case with both categorical and continuous features.

## 4.4 Simulation study with only categorical features

We evaluate the performance of the different Shapley value approximation methods in the case of only categorical features with six different experimental setups. Table 2 describes these different experiments.

In each experiment, we sample $n_{train} = 1000$ training observations using the approach from Section 4.2, where the mean $\mu$ is 0 and the covariance matrix is constructed with $\Sigma_{j,j} = 1$, $j = 1, \dots, M$, $\Sigma_{i,j} = \rho$, for $i \neq j$, where $\rho \in \{0, 0.1, 0.3, 0.5, 0.8, 0.9\}$. We set the response to

$$y_i = \alpha + \sum_{j=1}^{M} \sum_{l=2}^{L} \beta_{jl} \mathbf{1}(x_{ij} = l) + \varepsilon_i, \quad (14)$$

where $x_{ij}$ is the jth feature of the ith training observation, $\varepsilon_i$, $i = 1, \dots, n_{train}$, are i.i.d. random variables sampled from the distribution $N(0, 0.01)$, and $\alpha, \beta_{jl}$, $j = 1, \dots, M$, $l = 1, \dots, L$ are parameters sampled from $N(0, 1)$, which are fixed for every experiment. The predictive model, $f(\cdot)$, takes the same form without the noise term (i.e. (9) with $C_{cont} = \emptyset$), where the parameters are fit to the $n_{train}$ training observations using standard linear regression. Then, we estimate the Shapley values using the different methods from Table 1.

Table 2 shows that only ctree and the independence method are used when $M > 4$. This is because for $M > 4$, the methods that require one-hot encoding are too computationally expensive. In the same three cases, the number of unique feature combinations ($M^L$) is so large that a subset of the $T = 2000$ most likely feature combinations are used instead – see the discussion related to (10).

<table><caption>**Table 2.** An outline of the simulation study when using only categorical features. $M$ denotes the number of features, $L$ denotes the number of categories, and $T$ denotes the number of unique test observations used to compute (10).</caption><thead><tr><th>M</th><th>L</th><th>T</th><th>Categorical cut-off values</th><th>Methods used</th></tr></thead><tbody><tr><td>3</td><td>3</td><td>27</td><td>(-∞, 0, 1, ∞)</td><td>all</td></tr><tr><td>3</td><td>4</td><td>81</td><td>(-∞, -0.5, 0, 1, ∞)</td><td>all</td></tr><tr><td>4</td><td>3</td><td>64</td><td>(-∞, 0, 1, ∞)</td><td>all</td></tr><tr><td>5</td><td>6</td><td>2000</td><td>(-∞, -0.5, -0.25, 0, 0.9, 1, ∞)</td><td>ctree, independence</td></tr><tr><td>7</td><td>5</td><td>2000</td><td>(-∞, -0.5, -0.25, 0, 1, ∞)</td><td>ctree, independence</td></tr><tr><td>10</td><td>4</td><td>2000</td><td>(-∞, -0.5, 0, 1, ∞)</td><td>ctree, independence</td></tr></tbody></table>

The results of these experiments are shown in Table 3. When the dependence between the features is small, the performance of each method is almost the same.
Note that when the correlation, $\rho$, is 0 (i.e. the features are independent) and $M \le 4$, the ctree and independence methods perform equally in terms of MAE, and, in fact, give identical Shapley values. This is because when $\rho = 0$, ctree never rejects the hypothesis of independence when fitting any of the trees. As a result, ctree weighs all training observations equally which is analogous to the independence method. This ability to adapt the complexity of the dependence modeling to the actual dependence in the data is a major advantage of the ctree approach. When $M > 4$, the results of the independence and ctree methods for $\rho = 0$ are slightly different. The reason is that when the dimension is large, ctree tests many more hypotheses and therefore is more likely to reject some of the hypotheses. Since the independence method performs better than ctree in these cases, this suggests that the parameter $\alpha$ could be reduced for higher dimensions to improve the performance in low-correlation settings. This remains to be investigated, however.

As expected, the ctree method outperforms the independence method unless the dependence between the features is very small. The ctree approach also always outperforms (albeit marginally) the empirical, Gaussian, and ctree-onehot approaches. In addition, a major advantage of using ctree is that it does not require one-hot encoding. Since the computational complexity of computing Shapley values grows exponentially in the number of features (one-hot encoded or not), the computation time for methods requiring one-hot-encoding grows quickly compared to ctree. In Table 4, we show the average run time (in seconds) per test observation of each method. The average is taken over all correlations since the computation times are almost the same for each correlation.

The empirical method is the fastest amongst the one-hot encoded methods and is still between two and five times slower than the ctree method. This means that if the number of features/categories is large, using one-hot encoding is not suitable. For the Gaussian method we calculate (7) using both 100 and 1000 samples from the conditional distribution. Table 3 shows that the MAE is slightly smaller in the latter case but from Table 4, we see that it is nearly three times slower. Such a small performance increase is probably not worth the extra computation time.

## 4.5 Simulation study with both categorical and continuous features

We also perform a simulation study with both categorical and continuous features. Because we need to use numerical integration to calculate the true Shapley values (see Section 4.3), the computational complexity is large even for lower-dimensional settings. Therefore, we restrict ourselves to an experiment with two categorical features with $L = 4$ categories each and two continuous features. Unless otherwise mentioned, the simulation setup follows that of Section 4.4. As described in Section 4.2, we simulate dependent categorical/continuous data by only transforming two of the four original Gaussian features. The cut-off vector for the categorical features is set to $(-\infty, -0.5, 0, 1, \infty)$. Similarly to (14), the
<table><caption>**Table 3.** The MAE of each method and correlation, $\rho$, for each experiment. The bolded numbers denote the smallest MAE per experiment and $\rho$.</caption><thead><tr><th rowspan="2">M</th><th rowspan="2">L</th><th rowspan="2">Method</th><th colspan="6">ρ</th></tr><tr><th>0</th><th>0.1</th><th>0.3</th><th>0.5</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><td rowspan="6">3</td><td rowspan="6">3</td><td>empirical</td><td>0.0308</td><td>0.0277</td><td>0.0358</td><td>0.0372</td><td>0.0419</td><td>0.0430</td></tr><tr><td>Gaussian (100)</td><td>0.0308</td><td>0.0237</td><td>0.0355</td><td>0.0330</td><td>0.0320</td><td>0.0384</td></tr><tr><td>Gaussian (1000)</td><td>0.0307</td><td>0.0236</td><td>0.0354</td><td>0.0327</td><td>0.0318</td><td>0.0383</td></tr><tr><td>ctree-onehot</td><td>0.0278</td><td>0.0196</td><td>0.0345</td><td>0.0363</td><td>0.0431</td><td>0.0432</td></tr><tr><td>ctree</td><td><b>0.0274</b></td><td><b>0.0191</b></td><td><b>0.0302</b></td><td><b>0.0310</b></td><td><b>0.0244</b></td><td><b>0.0259</b></td></tr><tr><td>independence</td><td><b>0.0274</b></td><td><b>0.0191</b></td><td>0.0482</td><td>0.0777</td><td>0.1546</td><td>0.2062</td></tr><tr><td rowspan="6">3</td><td rowspan="6">4</td><td>empirical</td><td>0.0491</td><td>0.0465</td><td>0.0447</td><td>0.0639</td><td>0.0792</td><td>0.0659</td></tr><tr><td>Gaussian (100)</td><td>0.0402</td><td>0.0350</td><td>0.0358</td><td>0.0620</td><td>0.0762</td><td>0.0724</td></tr><tr><td>Gaussian (1000)</td><td>0.0403</td><td>0.0353</td><td>0.0361</td><td>0.0624</td><td>0.0763</td><td>0.0738</td></tr><tr><td>ctree-onehot</td><td>0.0324</td><td>0.0244</td><td>0.0429</td><td>0.0617</td><td>0.0808</td><td>0.0680</td></tr><tr><td>ctree</td><td><b>0.0318</b></td><td>0.0331</td><td><b>0.0369</b></td><td><b>0.0422</b></td><td><b>0.0416</b></td><td><b>0.0291</b></td></tr><tr><td>independence</td><td><b>0.0318</b></td><td><b>0.0283</b></td><td>0.0774</td><td>0.1244</td><td>0.2060</td><td>0.2519</td></tr><tr><td rowspan="6">4</td><td rowspan="6">3</td><td>empirical</td><td>0.0385</td><td>0.0474</td><td>0.0408</td><td>0.0502</td><td>0.0473</td><td>0.0389</td></tr><tr><td>Gaussian (100)</td><td>0.0312</td><td>0.0381</td><td>0.0327</td><td>0.0459</td><td>0.0475</td><td>0.0409</td></tr><tr><td>Gaussian (1000)</td><td>0.0312</td><td>0.0385</td><td>0.0330</td><td>0.0453</td><td>0.0480</td><td>0.0410</td></tr><tr><td>ctree-onehot</td><td>0.0234</td><td>0.0305</td><td>0.0402</td><td>0.0530</td><td>0.0484</td><td>0.0397</td></tr><tr><td>ctree</td><td><b>0.0223</b></td><td>0.0414</td><td><b>0.0387</b></td><td><b>0.0453</b></td><td><b>0.0329</b></td><td><b>0.0253</b></td></tr><tr><td>independence</td><td><b>0.0223</b></td><td><b>0.0355</b></td><td>0.0961</td><td>0.1515</td><td>0.2460</td><td>0.2848</td></tr><tr><td rowspan="2">5</td><td rowspan="2">6</td><td>ctree</td><td>0.0237</td><td>0.0492</td><td><b>0.0621</b></td><td><b>0.0760</b></td><td><b>0.0767</b></td><td><b>0.0899</b></td></tr><tr><td>independence</td><td><b>0.0222</b></td><td><b>0.0469</b></td><td>0.1231</td><td>0.1803</td><td>0.2835</td><td>0.3039</td></tr><tr><td rowspan="2">7</td><td rowspan="2">5</td><td>ctree</td><td>0.0209</td><td><b>0.0333</b></td><td><b>0.0402</b></td><td><b>0.0542</b></td><td><b>0.0530</b></td><td><b>0.0559</b></td></tr><tr><td>independence</td><td><b>0.0193</b></td><td>0.0345</td><td>0.0794</td><td>0.1294</td><td>0.1908</td><td>0.2397</td></tr><tr><td rowspan="2">10</td><td rowspan="2">4</td><td>ctree</td><td>0.0169</td><td><b>0.0505</b></td><td><b>0.0617</b></td><td><b>0.0607</b></td><td><b>0.0627</b></td><td><b>0.0706</b></td></tr><tr><td>independence</td><td><b>0.0153</b></td><td>0.0544</td><td>0.1593</td><td>0.2180</td><td>0.3017</td><td>0.3412</td></tr></tbody></table>

response is given by

$$y_i = \alpha + \sum_{j=1}^{2} \sum_{l=2}^{4} \beta_{jl} \mathbf{1}(x_{ij} = l) + \sum_{j=3}^{4} \gamma_j x_{ij} + \varepsilon_i,$$

for which we fit a linear regression model of the same form without the error term to act as the predictive model $f(\cdot)$.

Then, we estimate the Shapley values using the methods from Table 1 except for the Gaussian method with 1000 samples. This method is excluded since Section 4.4 showed that its performance was very similar to that of the Gaussian method with 100 samples but significantly more time consuming. To compare the performance of the different methods, we sample $T = 500$ observations from the joint distribution of the features and compute the MAE using (10) as described in Section 4.1.

The results are displayed in Table 5. The Gaussian method is the best performing method when $\rho = 0.1, 0.3, 0.5$ while the empirical and ctree methods are the best performing when $\rho = 0.8$ and $\rho = 0.9$, respectively. The results are not
<table><caption>**Table 4.** The mean run time (in seconds) per test observation, *T*, where the mean is taken over all correlations, *ρ*.</caption><thead><tr><th>M</th><th>L</th><th>T</th><th>Method</th><th>Mean time<br>per test obs</th></tr></thead><tbody><tr><td rowspan="6">3</td><td rowspan="6">3</td><td rowspan="6">27</td><td>empirical</td><td>0.086</td></tr><tr><td>Gaussian (100)</td><td>4.833</td></tr><tr><td>Gaussian (1000)</td><td>13.295</td></tr><tr><td>ctree-onehot</td><td>0.338</td></tr><tr><td>ctree</td><td>0.040</td></tr><tr><td>independence</td><td>0.013</td></tr><tr><td rowspan="6">3</td><td rowspan="6">4</td><td rowspan="6">64</td><td>empirical</td><td>0.553</td></tr><tr><td>Gaussian (100)</td><td>8.041</td></tr><tr><td>Gaussian (1000)</td><td>29.160</td></tr><tr><td>ctree-onehot</td><td>1.807</td></tr><tr><td>ctree</td><td>0.023</td></tr><tr><td>independence</td><td>0.007</td></tr><tr><td rowspan="6">4</td><td rowspan="6">3</td><td rowspan="6">81</td><td>empirical</td><td>0.293</td></tr><tr><td>Gaussian (100)</td><td>3.845</td></tr><tr><td>Gaussian (1000)</td><td>12.983</td></tr><tr><td>ctree-onehot</td><td>0.841</td></tr><tr><td>ctree</td><td>0.052</td></tr><tr><td>independence</td><td>0.012</td></tr><tr><td rowspan="2">5</td><td rowspan="2">6</td><td rowspan="2">2000</td><td>ctree</td><td>0.118</td></tr><tr><td>independence</td><td>0.030</td></tr><tr><td rowspan="2">7</td><td rowspan="2">5</td><td rowspan="2">2000</td><td>ctree</td><td>0.590</td></tr><tr><td>independence</td><td>0.158</td></tr><tr><td rowspan="2">10</td><td rowspan="2">4</td><td rowspan="2">2000</td><td>ctree</td><td>6.718</td></tr><tr><td>independence</td><td>2.066</td></tr></tbody></table>

surprising since we only have two categorical features. With more categorical features or categories, we expect that the ctree method would outperform the other ones when $\rho$ is not small. We also show the run time of each method in Table 6. The one-hot encoded methods are between nine and 75 times slower than the ctree method. This demonstrates, again, the value of using the ctree method when estimating Shapley values with categorical features.

## 5 Real data example

Although there is a growing literature of how to explain black-box models, there are very few studies that focus on quantifying the relevance of these methods [2]. This makes it difficult to compare different explainability methods on a real data set since there is no ground truth. One partial solution is to compare how different explainability models rank the same features for predictions based on specific test observations.
<table><caption>**Table 5.** The MAE of each method and correlation, $\rho$, for the experiment with two continuous and two categorical features ($L = 4$ categories each). The bolded numbers denote the smallest MAE per $\rho$.</caption><thead><tr><th rowspan="2">M</th><th rowspan="2">L</th><th rowspan="2">Method</th><th colspan="6">ρ</th></tr><tr><th>0</th><th>0.1</th><th>0.3</th><th>0.5</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><td rowspan="5">2 cont/2 cat 4</td><td>empirical</td><td></td><td>0.0853</td><td>0.0852</td><td>0.0898</td><td>0.0913</td><td><strong>0.0973</strong></td><td>0.1027</td></tr><tr><td>Gaussian (100)</td><td></td><td>0.0570</td><td><strong>0.0586</strong></td><td><strong>0.0664</strong></td><td><strong>0.0662</strong></td><td>0.1544</td><td>0.2417</td></tr><tr><td>ctree-onehot</td><td></td><td>0.0266</td><td>0.0714</td><td>0.1061</td><td>0.1024</td><td>0.1221</td><td>0.1188</td></tr><tr><td>ctree</td><td></td><td><strong>0.0093</strong></td><td>0.0848</td><td>0.1073</td><td>0.1060</td><td>0.0977</td><td><strong>0.0917</strong></td></tr><tr><td>independence</td><td></td><td><strong>0.0093</strong></td><td>0.0790</td><td>0.2178</td><td>0.3520</td><td>0.5524</td><td>0.6505</td></tr></tbody></table>

<table><caption>**Table 6.** The mean run time (in seconds) per test observation, $T$, where the mean is taken over all correlations, $\rho$. The simulation study has two continuous features and two categorical features ($L = 4$ categories each).</caption><thead><tr><th>M</th><th>L</th><th>T</th><th>Method</th><th>Mean time<br>per test obs</th></tr></thead><tbody><tr><td rowspan="5">4</td><td rowspan="5">4</td><td rowspan="5">500</td><td>empirical</td><td>0.758</td></tr><tr><td>Gaussian (100)</td><td>5.914</td></tr><tr><td>ctree-onehot</td><td>1.514</td></tr><tr><td>ctree</td><td>0.082</td></tr><tr><td>independence</td><td>0.057</td></tr></tbody></table>

In this section, we use a data set from the 2018 FICO Explainable Machine Learning Challenge [9] aimed at motivating the creation of explainable predictive models. The data set is of Home Equity Line of Credit (HELOC) applications made by homeowners. The response is a binary feature called RiskPerformance that takes the value 1 ('Bad') if the customer is more than 90 days late on his or her payment and 0 ('Good') otherwise. 52 percent of customers have the response 'Bad' and 48 percent have the response 'Good'. There are 23 features: 21 continuous and two categorical (with eight and nine categories, respectively) which can be used to model the probability of being a 'Bad' customer. Features with the value -9 are assumed missing. We remove the rows where all features are missing.

We first use this data set to compare the Shapley values calculated using the independence approach with those calculated with our ctree approach. Then, for a few test observations, we see how the Shapley explanations compare with the explanations from the 2018 FICO challenge Recognition Award winning team from Duke University [5] (hereafter referred to as just 'Duke').

After removing the missing data and a test set of 100 observations, we use the remaining 9,765 observations to train a 5-fold cross validation (CV) model using xgboost [6] and then average these to form the final model. Our model achieves an accuracy of 0.737 (compared to Duke's accuracy of 0.74). In our experience, explanation methods often behave differently when there is dependence
between the features. As a measure of dependence, we use the standard Pearson
correlation for all continuous features, Cramer's V correlation measure [7] for
categorical features, and the correlation ratio [10] for continuous and categorical
features. The feature 'MSinceMostRecentInqexcl7days' is the least correlated
with the rest of the features with correlations between -0.109 and 0.07. The
22 other features are strongly correlated with at least one other feature (max
absolute correlations between 0.4 and 0.99).

Turning to the Shapley value comparisons, we estimate the Shapley val-
ues of the features belonging to the 100 test observations using the indepen-
dence approach and the ctree approach. Then, we plot the Shapley value es-
timates against each other for a selection of four features in Figure 2. The
top left panel shows the Shapley values of one of the two categorical features
('MaxDelq2PublicRecLast12M'). Both methods give Shapley values fairly close
to 0 for this feature, but there are some differences. The top right panel shows an
example of a feature ('ExternalRiskEstimate') where the two methods estimate
quite different Shapley values. Since these are some of the largest (absolute)
Shapley values, this is one of the most influential features. We also see that
for most test observations, the independence method estimates more extreme
Shapley values than the ctree method.

The bottom left panel shows a feature ('MSinceMostRecentInqexcl7days')
where the two methods estimate relatively similar Shapley values. As noted
above, this feature is the least correlated with the rest of the features. We believe
the two methods behave similarly because the independence method performs
best when dealing with nearly independent features. Finally, the bottom right
panel shows a feature ('NumTrades90Ever2DerogPubRec') where the indepen-
dence method assigns most test observations a Shapley value very close to 0
while the ctree method does not. Although not plotted, we see this trend for 6
out of the 23 features. We notice that each of these 6 features are highly cor-
related with at least one other feature (max absolute correlation between 0.46
and 0.99). We suspect that the methods behave differently because of the inde-
pendence method's failure to account for dependence between features. We also
colour three random test observations to show that for some test observations
(say 'green'), all Shapley values are estimated very similarly for the two methods,
while for others (say 'blue'), the methods are sometimes quite different.

We also attempt to compare explanations based on Shapley values (estimated
with either the independence or ctree method) with those based on Duke's ap-
proach [5] We reinforce that there is no ground truth when it comes to expla-
nations and that Duke does not use Shapley values in their solution. Therefore,
we only compare how these three methods rank feature importance for specific
test observations.

Duke's explainability approach is based on 10 smaller regression models fit
to 10 partitions (below called 'groups') of the features. They use a combination
of learned weights and risks to calculate the most influential group for each
customer/test observation. Based on our understanding, if the customer has
a large predicted probability of being 'Bad', the largest weight × risk is the
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** The Shapley values of 100 test observations calculated using both the independence and the ctree method for four of the 23 features.</figcaption></figure>

most influential group (given rank 1), while for small predicted probabilities, the largest weight divided by risk is given rank 1. It is not clear how they rank medium-range predictions.

We speculate that Duke does not properly account for feature dependence since they fit 10 independent models that only interact using an overall learned risk for each model. To test this hypothesis, we compare the group rankings of a few<fnref n="3" /> customers/test observations calculated by 1. Duke, 2. The Shapley approach under the independence assumption, and 3. Our new Shapley approach that uses conditional inference trees.

To calculate group importance based on Shapley values for a given test observation, we first estimate the Shapley values of each feature either either the independence or ctree approach. Then, we sum the Shapley values of the features belonging to the same group. This gives 10 new grouped Shapley values. Finally, we rank the grouped Shapley values by giving rank 1 to the group with

<footnote n="3">³ Duke's explanations were not readily available. For a given test observation, we had to manually input 23 feature values into a web application to get an explanation. Therefore, it was too time consuming to compare many test observations.</footnote>
the largest absolute grouped Shapley value. While the prediction models being
compared here are different (we use an xgboost model while Duke does not), the
two models have very similar overall performance and give similar predictions to
specific test observations. We believe this validates the rough comparison below.

We observe that for test observations with a large predicted probability of
being ‘Bad’, there is little pattern among the rankings calculated by the three
explanation methods. However, when Duke and independence give a group the
same ranking out of 10 (and ctree does not), we notice that this group includes
at least one feature that is very correlated with a feature in another group. On
the other hand, for test observations with a small predicted probability of being
‘Bad’, we notice that all three explanation methods rank the groups similarly.
Again, the only time ctree ranks a group differently than Duke and independence
(but Duke and independence rank similarly), the group includes at least one
feature highly correlated with a feature in another group.

# 6 Conclusion

The aim of this paper was to extend [14] and [1]’s Shapley methodology to ex-
plain mixed dependent features using conditional inference trees [12]. We showed
in two simulation studies that when the features are even mildly dependent, it is
advantageous to use our ctree method over the traditional independence method.
Although ctree often has comparable accuracy to some of [1]’s methods, those
methods require transforming the categorical features to one-hot encoded fea-
tures. We demonstrated that such one-hot encoding leads to a substantial in-
crease in computation time, making it infeasible in high dimensions.

We also demonstrated our methodology on a real financial data set. We
first compared the Shapley values of 100 test observations calculated using the
independence and ctree approaches. We noticed that the methods performed
similarly for an almost independent feature but otherwise performed quite dif-
ferently.

Then, we compared explanations based on the independence/ctree Shapley
approaches with those based on Duke's approach [5]. We had to fall back to com-
paring how the different methods ranked features rather than the explanations
themselves because there is no ground truth when it comes to explainability.
It was difficult to argue for one explanatory approach over another; however,
Duke's rankings seemed to agree more with the rankings based on Shapley val-
ues calculated under independence (which we saw was inaccurate in simulation
studies), than with our proposed ctree based Shapley value estimation method.

References

1. Aas, K., Jullum, M., Løland, A.: Explaining individual predictions when features are dependent: More accurate approximations to shapley values. arXiv preprint arXiv:1903.10464 (2019)
2. Adadi, A., Berrada, M.: Peeking inside the black-box: A survey on explainable artificial intelligence (xai). IEEE Access **6**, 52138–52160 (2018)

3. Arellano-Valle, R.B., Branco, M.D., Genton, M.G.: A unified view on skewed distributions arising from selections. Canadian Journal of Statistics **34**(4), 581–601 (2006)

4. Breiman, L., Friedman, J., Olshen, R., Stone, C.: Classification and regression trees. Chapman and Hall (1984)

5. Chen, C., Lin, K., Rudin, C., Shaposhnik, Y., Wang, S., Wang, T.: An interpretable model with globally consistent explanations for credit risk. arXiv preprint arXiv:1811.12615 (2018)

6. Chen, T., Guestrin, C.: XGBoost: A scalable tree boosting system. In: Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD). pp. 785–794. ACM (2016)

7. Cramir, H.: Mathematical methods of statistics. Princeton U. Press, Princeton p. 500 (1946)

8. Elith, J., Leathwick, J.R., Hastie, T.: A working guide to boosted regression trees. Journal of Animal Ecology **77**(4), 802–813 (2008)

9. FICO: Explainable machine learning challenge (2018), https://community.fico.com/s/explainable-machine-learning-challenge

10. Fisher, R.A.: Statistical methods for research workers. In: Breakthroughs in statistics, pp. 66–70. Springer (1992)

11. Genz, A., Bretz, F.: Computation of Multivariate Normal and t Probabilities. Lecture Notes in Statistics, Springer-Verlag, Heidelberg (2009)

12. Hothorn, T., Hornik, K., Zeileis, A.: Unbiased recursive partitioning: A conditional inference framework. Journal of Computational and Graphical statistics **15**(3), 651–674 (2006)

13. Hothorn, T., Zeileis, A.: partykit: A modular toolkit for recursive partytioning in R. Journal of Machine Learning Research **16**, 3905–3909 (2015)

14. Lundberg, S.M., Erion, G.G., Lee, S.I.: Consistent individualized feature attribution for tree ensembles. arXiv preprint arXiv:1802.03888 (2018)

15. Lundberg, S.M., Lee, S.I.: A Unified Approach to Interpreting Model Predictions. In: Advances in Neural Information Processing Systems. pp. 4768–4777. Curram Associates Inc. (2017)

16. Quinlan, J.R.: C4.5: Programs for Machine Learning. Morgan Kaufmann Publishers Inc. (1993)

17. R Core Team: R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (2019)

18. Sellereite, N., Jullum, M.: shapr: An r-package for explaining machine learning models with dependence-aware shapley values. Journal of Open Source Software **5**(46), 2027 (2020)

19. Shapley, L.S.: A Value for N-Person Games. Contributions to the Theory of Games **2**, 307–317 (1953)

20. Štrumbelj, E., Kononenko, I.: An Efficient Explanation of Individual Classifications using Game Theory. Journal of Machine Learning Research **11**, 1–18 (2010)

21. Štrumbelj, E., Kononenko, I.: Explaining prediction models and individual predictions with feature contributions. Knowledge and Information Systems **41**, 647–665 (2014)