doc/source/methods/TreeSHAP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[source]](../api/alibi.explainers.shap_wrappers.rst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tree SHAP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tree SHAP (**SH**apley **A**dditive ex**P**lanations) algorithm is based on the paper [From local explanations to global understanding with explainable AI for trees](https://www.nature.com/articles/s42256-019-0138-9) by Lundberg et al. and builds on the open source  [shap library](https://github.com/slundberg/shap) from the paper's first author.\n",
    "\n",
    "\n",
    "The algorithm provides human interpretable explanations suitable for regression and classification of models with tree structure applied to tabular data. This method is a member of the *additive feature attribution methods* class; feature attribution refers to the fact that the change of an outcome to be explained (e.g., a class probability in a classification problem) with respect to a *baseline* (e.g., average prediction probability for that class in the training set) can be attributed in different proportions to the model input features. \n",
    "\n",
    "A simple illustration of the explanation process is shown in Figure 1. Here we see depicted a tree-based model which takes as an input features such as `Age`, `BMI` or `Blood pressure` and outputs `Mortality risk score`, a continuous value. Let's assume that we aim to explain the difference between and observed outcome and no risk, corresponding to a base value of `0.0`. Using the Tree SHAP algorithm, we attribute the `4.0` difference to the input features. Because the sum of the attribute values equals `output - base value`, this method is _additive_. We can see for example that the `Sex` feature contributes negatively to this prediction whereas the remainder of the features have a positive contribution (i.e., increase the mortality risk). For explaining this particular data point, the `Blood Pressure` feature seems to have the largest effect, and corresponds to an increase in the mortality risk. See our example on how to perform explanations with this algorithm and visualise the results using the `shap` library visualisations [here](../examples/interventional_tree_shap_adult_xgb.ipynb) and [here](../examples/path_dependent_tree_shap_adult_xgb.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![overview](tree_exp_overview_img.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 1: Cartoon ilustration of explanation models with Tree SHAP.\n",
    "\n",
    "Image Credit: Scott Lundberg (see source [here](https://www.nature.com/articles/s42256-019-0138-9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to compute the shap values , the following arguments can optionally be set when calling the `explain` method:\n",
    "\n",
    "- `interactions`: set to `True` to decompose the shap value of every feature for every example into a main effect and interaction effects \n",
    "\n",
    "- `approximate`: set to `True` to calculate an approximation to shap values (see our [example](../examples/path_dependent_tree_shap_adult_xgb.ipynb))\n",
    "\n",
    "- `check_additivity`: if the explainer is initialised with `model_output = raw` and this option is `True` the explainer checks that the sum of the shap values is equal to model output - expected value\n",
    "\n",
    "- `tree_limit`: it an `int` is passed, an ensemble formed of only `tree_limit` trees is explained\n",
    "\n",
    "If the dataset contains categorical variables that have been encoded before being passed to the explainer and a single shap value is desired for each categorical variable, the the following options should be specified:\n",
    "\n",
    "\n",
    "- `summarise_result`: set to `True`\n",
    "\n",
    "- `cat_var_start_idx`: a sequence of integers containing the column indices where categorical variables start. If the feature matrix contains a categorical feature starting at index 0 and one at index 10, then `cat_var_start_idx=[0, 10]`\n",
    "\n",
    "- `cat_vars_enc_dim`: a list containing the dimension of the encoded categorical variables. The number of columns specified in this list is summed for each categorical variable starting with the corresponding index in `cat_var_start_idx`. So if `cat_var_start_idx=[0, 10]` and `cat_vars_enc_dim=[3, 5]`, then the columns with indices `0, 1` and `2` and `10, 11, 12, 13` and `14` will be combined to return one shap value for each categorical variable, as opposed to `3` and `5`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Path-dependent feature perturbation algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Initialiastion and fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The explainer is initialised with the following agruments:\n",
    "\n",
    "- a model, which could be an `sklearn`, `xgboost`, `catboost` or `lightgbm` model. Note that some of the models in these packages or models trained with specific objectives may not be supported. In particular, passing raw strings as categorical levels for `catboost` and `lightgbm` is not supported\n",
    "\n",
    "- `model_output` should always default to `raw` for this algorithm\n",
    "\n",
    "- optionally, set `task` to `'classification'` or `'regression'` to indicate the type of prediction the model makes. If set to `regression` the `prediction` field of the response is empty\n",
    "\n",
    "- optionally, a list of feature names via `feature_names`. This is used to provide information about feature importances in the response\n",
    "\n",
    "- optionally, a dictionary, `category_names`, that maps the columns of the categorical variables to a list of strings representing the names of the categories. This may be used for visualisation in the future.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python \n",
    "from alibi.explainers import TreeShap\n",
    "\n",
    "explainer = TreeShap(\n",
    "    model, \n",
    "    feature_names=['size', 'age'], \n",
    "    categorical_names={0: ['S', 'M', 'L', 'XL', 'XXL']}\n",
    ")\n",
    "```\n",
    "\n",
    "For this algorithm, fit is called with no arguments:\n",
    "\n",
    "```python\n",
    "explainer.fit()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Explanation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To explain an instance `X`, we simply pass it to the explain method:\n",
    "\n",
    "```python \n",
    "explanation = explainer.explain(X)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The returned explanation object has the following fields:\n",
    "\n",
    "* `explanation.meta`: \n",
    "\n",
    "```python\n",
    "{'name': 'TreeShap',\n",
    " 'type': ['whitebox'],\n",
    " 'task': 'classification',\n",
    " 'explanations': ['local', 'global'],\n",
    " 'params': {'summarise_background': False, 'algorithm': 'tree_path_dependent' ,'kwargs': {}}\n",
    "}\n",
    "```\n",
    "\n",
    "This field contains metadata such as the explainer name and type as well as the type of explanations this method can generate. In this case, the `params` attribute shows the Tree SHAP variant that will be used to explain the model in the `algorithm` attribute.\n",
    "\n",
    "* `explanation.data`:\n",
    "\n",
    "```python\n",
    "data={'shap_values': [\n",
    "        array([[ 5.0661433e-01,  2.7620478e-02],\n",
    "       [-4.1725192e+00,  4.4859368e-03],\n",
    "       [ 4.1338313e-01, -5.5618007e-02]],\n",
    "      dtype=float32)\n",
    "    ],\n",
    "    'shap_interaction_values': [array([], dtype=float64)],\n",
    "    'expected_value': array([-0.06472124]),\n",
    "    'model_output': 'raw',\n",
    "    'categorical_names': {0: ['S', 'M', 'L', 'XL', 'XXL']},\n",
    "    'feature_names': ['size', 'age'],\n",
    "    'raw': {\n",
    "        'raw_prediction': array([-0.73818872, -8.8434663 , -3.24204564]),\n",
    "        'loss': [],\n",
    "        'prediction': array([0, 0, 0]),\n",
    "        'instances': array([[0, 23],\n",
    "       [4,  55],\n",
    "       [2,  43]]),\n",
    "        'labels': array([], dtype=float64),\n",
    "        'importances': {\n",
    "            '0': {\n",
    "                'ranked_effect': array([1.6975055 , 1.3598266], dtype=float32),\n",
    "                'names': [\n",
    "                    'size',\n",
    "                    'age',\n",
    "                ]\n",
    "            },\n",
    "            'aggregated': {\n",
    "                'ranked_effect': array([1.6975055 , 1.3598266], dtype=float32),\n",
    "                'names': [\n",
    "                    'size',\n",
    "                    'age',\n",
    "                ]\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This field contains:\n",
    "\n",
    "* `shap_values`: a list of length equal to the number of model outputs, where each entry is an array of dimension samples x features of shap values. For the example above , 3 instances with 2 features has been explained so the shap values for each class are of dimension 3 x 2\n",
    "\n",
    "* `shap_interaction_values`: an empty list since this `interactions` was set to `False` in the explain call\n",
    "\n",
    "* `expected_value`: an array containing expected value for each model output \n",
    "\n",
    "* `model_output`: `raw` indicates that the model raw output was explained, the only option for the path dependent algorithm\n",
    "\n",
    "* `feature_names`: a list with the feature names\n",
    "\n",
    "* `categorical_names`: a mapping of the categorical variables (represented by indices in the shap_values columns) to the description of the category\n",
    "\n",
    "* `raw`: this field contains:\n",
    "\n",
    "    * `raw_prediction`: a samples x n_outputs array of predictions for each instance to be explained. \n",
    "\n",
    "    * `prediction`: an array containing the index of the maximum value in the `raw_prediction` array\n",
    "\n",
    "    * `instances`: a samples x n_features array of instances which have been explained\n",
    "    \n",
    "    * `labels`: an array containing the labels for the instances to be explained \n",
    "\n",
    "    * `importances`: a dictionary where each entry is a dictionary containing the sorted average magnitude of the shap value (ranked_effect) along with a list of feature names corresponding to the re-ordered shap values (names). There are n_outputs + 1 keys, corresponding to n_outputs and the aggregated output (obtained by summing all the arrays in shap_values)\n",
    "\n",
    "     \n",
    "Please see our examples on how to visualise these outputs using the shap library visualisations library visualisations [here](../examples/interventional_tree_shap_adult_xgb.ipynb) and [here](../examples/path_dependent_tree_shap_adult_xgb.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Shapley interaction values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Initialisation and fit "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shapley interaction values can only be calculated using the path-dependent feature perturbation algorithm in this release, so no arguments are passed to the `fit` method:\n",
    "\n",
    "```python \n",
    "explainer = TreeShap(\n",
    "    model, \n",
    "    model_output='raw',\n",
    ")\n",
    "\n",
    "explainer.fit()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Explanation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain the Shapley interaction values, the `explain` method is called with the option `interactions=True`:\n",
    "\n",
    "```python\n",
    "explanation = explainer.explain(X, interactions=True)\n",
    "```\n",
    "\n",
    "The explanation contains a list with the shap interaction values for each model output in the `shap_interaction_values` field of the `data` property."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interventional feature perturbation algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Explaining model output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Initialiastion and fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "explainer = TreeShap(\n",
    "    model, \n",
    "    model_output='raw',\n",
    ")\n",
    "\n",
    "explainer.fit(X_reference)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model output can be set to `model_output='probability'` to explain models which return probabilities. Note that this requires the model to be trained with specific objectives. Please the footnote to our path-dependent feature perturbation [example](../examples/path_dependent_tree_shap_adult_xgb.ipynb) for an example of how to set the model training objective in order to explain probability outputs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Explanation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To explain instances in `X`, the explainer is called as follows:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```pyhton\n",
    "explanation = explainer.explain(X)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Explaining loss functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Initialisation and fit "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To explain loss function, the following configuration and fit steps are necessary:\n",
    "\n",
    "```python \n",
    "explainer = TreeShap(\n",
    "    model, \n",
    "    model_output='log_loss',\n",
    ")\n",
    "\n",
    "explainer.fit(X_reference)\n",
    "```\n",
    "\n",
    "Only square loss regression objectives and cross-entropy classification objectives are supported in this release."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Explanation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the labels need to be passed to the `explain` method in order to obtain the explanation:\n",
    "\n",
    "```python\n",
    "explanation = explainer.explain(X, y)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Miscellaneous"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Runtime considerations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Adjusting the size of the reference dataset "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm automatically warns the user if a background dataset size of more than `1000` samples is passed. If the runtime of an explanation with the original dataset is too large, then the algorithm can automatically subsample the background dataset during the `fit` step. This can be achieve by specifying the fit step as \n",
    "\n",
    "```python\n",
    "explainer.fit(\n",
    "    X_reference,\n",
    "    summarise_background=True,\n",
    "    n_background_samples=300,\n",
    ")\n",
    "```\n",
    "\n",
    "or \n",
    "```python\n",
    "explainer.fit(\n",
    "    X_reference,\n",
    "    summarise_background='auto'\n",
    ")\n",
    "```\n",
    "\n",
    "The `auto` option will select `1000` examples, whereas using the boolean argument allows the user to directly control the size of the reference set. If categorical variables are specified, the algorithm uses subsampling of the data. Otherwise, a kmeans clustering algorithm is used to select the background dataset.\n",
    "\n",
    "As describe above, the explanations are performed with respect to the expected output over this dataset so the shap values will be affected by the dataset selection. We recommend experimenting with various ways to choose the background dataset before deploying explanations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theoretical overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Recall that, for a model $f$, the Kernel SHAP algorithm [[1]](#References) explains a certain outcome with respect to a chosen reference (or an expected value) by estimating the shap values of each feature $i$  from $\\{1, ..., M\\}$, as follows:\n",
    "\n",
    "- enumerate all subsets $S$ of the set $F   \\setminus \\{i\\}$ \n",
    "\n",
    "- for each $S \\subseteq F \\setminus \\{i\\}$, compute the contribution of feature $i$ as $C(i|S) = f(S \\cup \\{i\\}) - f(S)$\n",
    "\n",
    "- compute the shap value according to\n",
    "\n",
    "$$\n",
    "\\phi_i := \\frac{1}{M} \\sum \\limits_{{S \\subseteq F \\setminus \\{i\\}}} \\frac{1}{ M - 1 \\choose |S|} C(i|S).         \n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "\n",
    "Since most models do not accept arbitrary patterns of missing values at inference time, $f(S)$ needs to be approximated. The original formulation of the Kernel Shap algorithm [[1]](#References) proposes to compute $f(S)$ as the _observational conditional expectation_\n",
    "\n",
    "$$\n",
    "f(S) := \\mathbb{E}\\left[f(\\mathbf{x}_{S}, \\mathbf{X}_{\\bar{S}} | \\mathbf{X}_S = \\mathbf{x}_S) \\right]\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "where the expectation is taken over a *background dataset*, $\\mathcal{D}$, after conditioning. Computing this expectation involves drawing sufficiently many samples from $\\mathbf{X}_{\\bar{S}}$ for every sample from $\\mathbf{X}_S$, which is expensive. Instead, $(2)$ is approximated by \n",
    "\n",
    "$$\n",
    "f(S) := \\mathbb{E} \\left[f(\\mathbf{x}_{S}, \\mathbf{X}_{\\bar{S}})\\right]\n",
    "$$\n",
    "\n",
    "where features in a subset $S$ are fixed and features in $\\bar{S}$ are sampled from the background dataset. This quantity is referred to as _marginal_ or *interventional conditional expectation*, to emphasise that setting features in $S$ to the values $\\mathbf{x}_{S}$ can be viewed as an intervention on the instance to be explained.\n",
    "\n",
    "As described in [[2]](#References), if estimating impact of a feature $i$ on the function value by $\\mathbb{E} \\left[ f | X_i = x_i \\right]$, one should bear in mind that observing $X_i = x_i$ changes the distribution of the features $X_{j \\neq i}$ if these variables are correlated. Hence, if the conditional expectation if used to estimate $f(S)$, the Shapley values might not be accurate since they also depend on the remaining variables, effect which becomes important if there are strong correlations amongst the independent variables. Furthermore, the authors show that estimating $f(S)$ using the conditional expectation violates the *sensitivity principle*, according to which the Shapley value of a redundant variable should be 0. On the other hand, the intervention breaks the dependencies, ensuring that the sensitivity holds. One potential drawback of this method is that setting a subset of values to certain values without regard to the values of the features in the complement (i.e., $\\bar{S}$) can generate instances that are outside the training data distribution, which will affect the model prediction and hence the contributions.\n",
    "\n",
    "The following sections detail how these methods work and how, unlike Kernel SHAP, compute the exact shap values in polynomial time. The algorithm estimating contributions using interventional expectations is presented, with the remaining sections being dedicated to presenting an approximate algorithm for evaluating the interventional expectation that does not require a background dataset and Shapley interaction values.\n",
    "\n",
    "<a id='source_1'></a>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interventional feature perturbation\n",
    "<a id='interventional'></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interventional feature perturbation algorithm provides an efficient way to calculate the expectation $f(S) := \\mathbb{E} \\left[f(\\mathbf{x}_{S}, \\mathbf{X}_{\\bar{S}})\\right]$ for all possible subsets $S$, and to combine these values according to equation $(1)$ in order to obtain the Shapley value. Intuitively, one can proceed as follows:\n",
    "\n",
    "- choose a background sample $r \\in \\mathcal{D}$\n",
    "\n",
    "- for each feature $i$, enumerate all subsets $S \\subseteq F \\setminus \\{i\\}$\n",
    "\n",
    "- for each such subset, $S$, compute $f(S)$ by traversing the tree with a _hybrid sample_ where the features in $\\bar{S}$ are replaced by their corresponding values in $r$\n",
    "\n",
    "- combine results according to equation $(1)$\n",
    "\n",
    "If $R$ samples from the background distribution are used, then the complexity of this algorithm is $O(RM2^M)$ since we perform $2^M$ enumerations for each of the $M$ features, $R$ times. The key insight into this algorithm is that multiple hybrid samples will end up traversing identical paths and that this can be avoided if the shap values' calculation is reformulated as a summation over the paths in the tree (see [[4]](#References) for a proof):\n",
    "\n",
    "$$\n",
    "\\phi_i = \\sum_{P}\\phi_{i}^P\n",
    "$$\n",
    "\n",
    "where the summation is over paths $P$ in the tree descending from $i$. The value and sign of the contribution of each path descending through a node depends on whether the split from the node is due to a foreground or a background feature, as explained in the practical example below.\n",
    "\n",
    "<a id='source_4'></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing contributions with interventional Tree SHAP: a practical example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![img](interventional_tree_shap_f.png)\n",
    "Figure 2: Ilustration of the feature contribution and expected value estimation process using interventional perturbation Tree SHAP. The positive and the negative contributions of a node are represented in <span style=\"color:green\">green</span> and <span style=\"color:red\">red</span>, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure above, the paths followed due the instance to be explained $x$ are coloured in red, paths followed due to the background sample in red, and common paths in yellow.\n",
    "\n",
    "The instance to be explained is perturbed using a reference sample by the values of the features $F1$, $F3$ and $F5$ in $x$ with the corresponding values in $r$. This process gives the name of the algorithm since following the paths indicated by the background sample is akin to intervening on the instance to be explained with features from the background sample. Therefore, one defines the set $F$ in the previous section as $F = \\{ j: x_{j} \\neq r_{j}\\}$ for this case. Note that these are the only features for which one can estimate a contribution given this background sample; the same path is followed for features $F2$ and $F4$ for both the original and the perturbed sample, so these features do not contribute to explain the difference between the observed outcome ($v_6$) and the outcome that would have been observed if the tree had been traversed according to the reference $(v_{10})$.\n",
    "\n",
    "\n",
    "Considering the structure of the tree for the given $x$ and $r$ together with equation $(1)$ reveals that the left subtree can be traversed to compute the negative terms in the summation whereas the right subtree will provide positive terms. This is because the nodes in the left subtree can only be reached if $F1$ takes the value from the background sample, that is, only $F1$ is missing. Because $F2$ and $F4$ do not contribute to explaining $f(x) - f(r)$, the negative contribution of the left subtree will be equal to the negative contribution of node $8$. This node sums two negative components: one when the downstream feature $F5$ is also missing (corresponding to evaluating $f$ at $S = \\varnothing$) and one when $F5$ is present (corresponding to evaluating $f$ at $S=\\{F5\\}$). These negative values are weighted according to the combinatorial factor in equation $(1)$.  By a similar reasoning, the nodes in the right subtree are reached only if $F1$ is present and they provide the positive terms for the shap value computation. Note that the combinatorial factor in $(1)$ should be evaluated with $|S| \\gets |S| - 1$ for positive contributions since $|S|$ is increased by $1$ because of the feature whose contribution is calculated is present in the right subtree. \n",
    "\n",
    "A similar reasoning is applied to compute the contributions of the downstream nodes. For example, to estimate the contribution of $F5$, one considers a set $S = \\varnothing$ and observes the value of node $10$, and weighs that with the combinatorial factor from equation  $(1)$ where $M-1 = 1$ and $|S|=0$ (because there are no features present on the path) and a positive contribution from node $9$ weighted by the same combinatorial factor (because $S = \\{F5\\}$ so $|S| - 1 = 0$). \n",
    "\n",
    "To summarise, the efficient algorithm relies on the following key ideas:\n",
    "\n",
    "- each node in the tree is assigned a positive contribution reflecting membership of the splitting feature in a subset $S$ and a negative contribution to indicate the feature is missing ($i\\in \\bar{S}$)\n",
    "\n",
    "- the positive and negative contributions of a node can be computed by summing the positive and negative contributions  of the children nodes, in keeping with the fact that the Shapley value can be computed by summing a contribution from each path the feature is on \n",
    "\n",
    "- to compute the contribution of a feature at a node, one adds a positive contribution from the node reached by splitting on the feature from the instance to be explained and a negative contribution from the node reached by splitting on the feature in the background sample\n",
    "\n",
    "- features for which the instance to be explained and the reference follow the same path are assigned $0$ contribution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Explaining loss functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One advantage of the interventional approach is that it allows to approximately transform the shap values to account for nonlinear transformation of outputs, such as the loss function. Recall that given $\\phi_i, ..., \\phi_M$ the local accuracy property guarantees that given $\\phi_0 = \\mathbb{E}[f(x)]$\n",
    "\n",
    "$$\n",
    "f(x) = \\phi_0 + \\sum \\limits_{i=1}^M \\phi_i.\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "Hence, in order to account for the effect of the nonlinear transformation $h$, one has to find the functions $g_0, ..., g_M$ such that\n",
    "\n",
    "$$\n",
    "h(f(x)) = g(\\phi_0) + \\sum \\limits_{i=1}^M g_i(\\phi_i)\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "For simplicity, let $y=h(x)$. Then using a first-order Taylor series expansion around $\\mathbb{E}[y]$ one obtains \n",
    "\n",
    "$$\n",
    "h(y) \\approx h(\\mathbb{E}[y]) + \\frac{\\partial h(y) }{\\partial y} \\Bigr|_{y=\\mathbb{E}[y]}(y - \\mathbb{E}[y]).\n",
    "\\tag{5}\n",
    "$$\n",
    "\n",
    "Substituting $(3)$ in $(5)$ and comparing coefficients with $(4)$ yields\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "g_0 & \\approx &  \\; h(\\mathbb{E}[y]) \\\\\n",
    "g_i &\\approx  &  \\; \\phi_i \\frac{\\partial h(y) }{\\partial y} \\Bigr|_{y=\\mathbb{E}[y]} .\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Hence, an approximate correction is given by simply scaling the shap values using the gradient of the nonlinear function. Note that in practice one may take the Taylor series expansion at a reference point $r$ from the background dataset and average over the entire background dataset to compute the scaling factor. This introduces an additional source of noise since $h(\\mathbb{E}[y]) = \\mathbb{E}[h(y)]$ only when $h$ is linear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computational complexity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a single foreground and background sample and a single tree, the algorithm runs in $O(LD)$ time. Thus, using $R$ background samples and a model containing $T$ trees, yields a complexity of $O(TRLD)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Path dependent feature perturbation\n",
    "<a id='path_dependent'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to approximate equation $(2)$ to compute $f(S)$ given an instance $x$ and a set of missing features $\\bar{S}$ is to recursively follow the decision path through the tree and:\n",
    "\n",
    "- return the node value if a split on a feature $i \\in S$ is performed\n",
    "\n",
    "- take a weighted average of the values returned by children if $i \\in \\bar{S}$, where the weighing factor is equal to the proportion of training examples flowing down each branch. This proportion is a property of each node, sometimes referred to as _weight_ or _cover_ and measures how important is that node with regard to classifying the training data.\n",
    "\n",
    "Therefore, in the path-dependent perturbation method, we compute the expectations with respect to the training data distribution by weighting the leaf values according to the proportion of the training examples that flow to that leaf.\n",
    "\n",
    "To avoid repeating the above recursion $M2^M$ times, one first notices that for a single decision tree, applying a perturbation would result in the sample ending up in a different leaf. Therefore, following each path from the root to a leaf in the tree is equivalent to perturbing subsets of features of varying cardinalities. Consequently, each leaf will contain a certain proportion of all possible subsets $S \\subseteq F$. Therefore, to compute the shap values, the following quantities are computed at each leaf, *for every feature $i$ on the path leading to that leaf*:\n",
    "\n",
    "- the proportion of subsets $S$ at the leaf that contain $i$ and the proportion of subsets $S$ that do not contain $i$\n",
    "\n",
    "- for each cardinality, the proportion of the sets of that cardinality contained at the leaf. Tracking each cardinality  as opposed to a single count of subsets falling into a given leaf is necessary since it allows to apply the weighting factor in equation (1), which depends on the subset size, $|S|$.\n",
    "\n",
    "This intuition can be summarised as follows:\n",
    "$$\n",
    "\\phi_i := \\sum \\limits_{j=1}^L \\sum \\limits_{P \\in {S_j}} \\frac {w(|P|, j)}{ M_j {M_j - 1 \\choose |P|}} (p_o^{i,j} - p_z^{i, j}) v_j\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "where $S_j$ is the set of present feature subsets at leaf $j$, $M_j$ is the length of the path and $w(|P|, j)$ is the proportion of all subsets of cardinality $P$ at leaf $j$,  $p_o^{i, j}$  and $p_z^{i, j}$ represent the fractions of subsets that contain or do not contain feature $i$ respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computational complexity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the above quantities, one can compute the _contribution_ of each leaf to the Shapley value of every feature. This algorithm has complexity $O(TLD^2)$ for an ensemble of trees where $L$ is the number of leaves, $T$ the number of trees in the ensemble and $D$ the maximum tree depth. If the tree is balanced, then $D=\\log L$ and the complexity of our algorithm is $O(TL\\log^2L)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Expected value for the path-dependent perturbation algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that although a background dataset is not provided, the expected value is computed using the node cover information, stored at each node. The computation proceeds recursively, starting at the root. The contribution of a node to the expected value of the tree is a function of the expected values of the children and is computed as follows:\n",
    "\n",
    "$$\n",
    "c_j = \\frac{c_{r(j)}r_{r(j)} + c_{l(j)}r_{l(j)}}{r_j}\n",
    "$$\n",
    "\n",
    "where $j$ denotes the node index, $c_j$ denotes the node expected value, $r_j$ is the cover of the $j$th node and $r(j)$ and $l(j)$ represent the indices of the right and left children, respectively. The expected value used by the tree is simply $c_{root}$. Note that for tree ensembles, the expected values of the ensemble members is weighted according to the tree weight and the weighted expected values of all trees are summed to obtain a single value.\n",
    "\n",
    "The cover depends on the objective function and the model chosen. For example, in a gradient boosted tree trained with squared loss objective, $r_j$ is simply the number of training examples flowing through $j$. For an arbitrary objective, this is the sum of the Hessian of the loss function evaluated at each point flowing through $j$, as explained [here](../examples/xgboost_model_fitting_adult.ipynb). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shapley interaction values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the Shapley values provide a solution to the problem of allocating a function variation  to the input features, in practice it might be of interest to understand how the importance of a feature depends on the other features. The Shapley interaction values can solve this problem, by allocating the change in the function amongst the individual features (*main effects*) and all pairs of features (*interaction effects*). Thus, they are defined as \n",
    "\n",
    "$$\n",
    "\\Phi_{i, j}(f, x) = \\sum_{S \\subseteq {F \\setminus \\{i, j\\}}} \\frac{1}{2|S| { M-1 \\choose |S| - 1}} \\nabla_{ij}(f, x, S), \\; i \\neq j\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\nabla_{ij}(f, x, S) = \\underbrace{f_{x}(S \\cup \\{i, j\\}) - f_x(S \\cup \\{j\\})}_{j \\; present} - \\underbrace{[f_x(S \\cup \\{i\\}) - f_x(S)]}_{j \\; not \\; present}. \n",
    "\\tag{8}\n",
    "$$\n",
    "\n",
    "Therefore, the interaction of features $i$ and $j$ can be computed by taking the difference between the shap values of $i$ when $j$ is present and when $j$ is not present. The main effects are defined as \n",
    "\n",
    "$$\n",
    "\\Phi_{i,i}(f, x) = \\phi_i(f, x) - \\sum_{i \\neq j} \\Phi_{i, j}(f, x),\n",
    "$$\n",
    "\n",
    "Setting $\\Phi_{0, 0} = f_x(\\varnothing)$ yields the local accuracy property for Shapley interaction values:\n",
    "\n",
    "$$f(x) = \\sum \\limits_{i=0}^M \\sum \\limits_{j=0}^M \\Phi_{i, j}.(f, x) $$.\n",
    "\n",
    "The interaction is split equally between feature $i$ and $j$, which is why the division by two appears in equation $(7)$. The total interaction effect is defined as $\\Phi_{i, j}(f, x) + \\Phi_{j, i}(f,x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computational complexity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to equation $(8)$, the interaction values can be computed by applying either the interventional or path-dependent feature perturbation algorithm twice: once by fixing the value of feature $j$ to $x_j$ and computing the shapley value for feature $i$ in this configuration, and once by fixing $x_j$ to a \"missing\" value and performing the same computation. Thus, the interaction values can be computed in $O(TMLD^2)$ with the path-dependent perturbation algorithm and $O(TMLDR)$ with the interventional feature perturbation algorithm.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison to other methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tree-based models are widely used in areas where model interpretability is of interest because node-level statistics gathered from the training data can be used to provide insights into the behaviour of the model across the training dataset, providing a _global explanation_ technique. As shown in our [example](../examples/path_dependent_tree_shap_adult_xgb.ipynb), considering different statistics gives rise to different importance rankings. As discussed in [[1]](#References) and [[3]](#References), depending on the statistic chosen, feature importances derived from trees are not *consistent*, meaning that a model where a feature is known to have a bigger impact might fail to have a larger importance. As such, feature importances cannot be compared across models.  In contrast, both the path-dependent and interventional perturbation algorithms tackle this limitation.\n",
    "\n",
    "In contrast to feature importances derived from tree statistics, the Tree SHAP algorithms can also provide local explanations, allowing the identification of features that are globally \"not important\", but can affect specific outcomes significantly, as might be the case in healthcare applications. Additionally, it provides a means to succinctly summarise the effect magnitude and direction (positive or negative) across potentially large samples. Finally, as shown in [[1]](#References) (see [here](https://static-content.springer.com/esm/art%3A10.1038%2Fs42256-019-0138-9/MediaObjects/42256_2019_138_MOESM1_ESM.pdf), p. 26),  averaging the instance-level shap values importance to derive a global score for each feature can result in improvements in feature selection tasks. \n",
    "\n",
    "Another method to derive instance-level explanations for tree-based model has been proposed by Sabaas [here](https://github.com/andosa/treeinterpreter). This feature attribution method is similar in spirit to Shapley value, but does not account for the effect of variable order as explained [here](https://static-content.springer.com/esm/art%3A10.1038%2Fs42256-019-0138-9/MediaObjects/42256_2019_138_MOESM1_ESM.pdf) (pp. 10-11) as well as not satisfying consistency ([[3]](#References)).\n",
    "\n",
    "Finally, both Tree SHAP algorithms exploit model structure to provide exact Shapley values computation albeit using different estimates for the effect of missing features, achieving explanations in low-order polynomial time. The KernelShap method relies on post-hoc (black-box) function modelling and approximations to approximate the same quantities and given enough samples has been shown to to the exact values (see experiments [here](https://static-content.springer.com/esm/art%3A10.1038%2Fs42256-019-0138-9/MediaObjects/42256_2019_138_MOESM1_ESM.pdf) and our [example](../examples/interventional_tree_shap_adult_xgb.ipynb)). Our Kernel SHAP [documentation](KernelSHAP.ipynb) provides comparisons of feature attribution methods based on Shapley values with other algorithms such as LIME and [anchors](Anchors.ipynb). \n",
    "\n",
    "<a id='source_3'></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='References'></a>\n",
    "\n",
    "[[1]](#source_1) Lundberg, S.M. and Lee, S.I., 2017. A unified approach to interpreting model predictions. In Advances in neural information processing systems (pp. 4765-4774).\n",
    "\n",
    "[[2]](#source_2) Janzing, D., Minorics, L. and Blöbaum, P., 2019. Feature relevance quantification in explainable AI: A causality problem. arXiv preprint arXiv:1910.13413.\n",
    "\n",
    "[[3]](#source_3) Lundberg, S.M., Erion, G.G. and Lee, S.I., 2018. Consistent individualized feature attribution for tree ensembles. arXiv preprint arXiv:1802.03888.\n",
    "\n",
    "[[4]](#source_4) Chen, H., Lundberg, S.M. and Lee, S.I., 2018. Understanding Shapley value explanation algorithms for trees. Under review for publication in Distill, draft available [here](https://hughchen.github.io/its_blog/index.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Path-dependent Feature Perturbation Tree SHAP "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Explaing tree models with path-dependent feature perturbation Tree SHAP](../examples/path_dependent_tree_shap_adult_xgb.nblink)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interventional Feature Perturbation Tree SHAP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Explaing tree models with path-dependent feature perturbation Tree SHAP](../examples/interventional_tree_shap_adult_xgb.nblink)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}