A Visual Introduction to Shapley Values

UCLouvain Applied Statistics Workshop

Author

Kris Sankaran

Published

June 19, 2026

Readings: 1, 2, 3, Code

Setup

Goal. Given a model $f$ and a sample $x \in \reals^{D}$, return a local feature attribution $\varphi_d\left(x\right)$ that quantifies the contribution of feature $d$ to the prediction $f\left(x\right)$.

Requirements.

Local feature attributions. Unlike global variable importances, $\varphi_d\left(x\right)$ are specific to sample $x$. This matters a stakeholder cares specifically about particular prediction $f(x)$ and wants an explanation for it.
Model agnostic. Attributions should be computable by querying $f$ alone, without any assumptions about what kind of model it is – it could be a black box.
Principled. The attribution measure should be derivable from a clear set of mathematical axioms.

Approach. SHAP values satisfy all three requirements. This handout focuses on their conceptual foundation (Section 1) and practical computation (Section 2).

Conceptual Foundation

Local Feature Attributions

High-stakes decisions. When an individual has a medical diagnosis made or an insurance claim denied, knowing the globally most important features isn’t enough. They deserve an explanation specific to their case.
- A related use case is algorithmic recourse. What could a stakeholder change to reverse a decision? (e.g., which change to their resume would have gotten them a job interview?)
Model debugging. A model might classify $x$ as a husky because it had snow in the background (large $\varphi_d(x)$ on pixels $d$ in the snow region) rather than the dog itself. This means the model has learned a “shortcut” and won’t generalize well – a wolf in the snow might get misclassified as a husky (Ribeiro et al. 2016).
Scientific discovery. In heterogeneous populations (e.g., different disease subtypes), a model might rely on different sets of features for each subpopulation. Local attributions can highlight these differences – e.g., identifying which features drive drug effectiveness in one subgroup vs. another.

Exercise: Give one example (hypothetical, or from your own experience) where local feature attribution would be useful. How would it differ from global variable importance?

Game Theory Definitions

SHAP is motivated by the credit assignment problem from game theory. Imagine agents $\mathcal{D} = \{1, \dots, D\}$ where any subset $S$ earns profit $v(S)$. How much of the total profit $v(\mathcal{D})$ should agent $d$ receive? This share is the Shapley value $\varphi_{d}(v)$.
Intuitively, agent $d$’s contribution depends on how much they add to each team $S$. Define the marginal contributions of agent $d$ to team $S$ as, \[C\left(d \vert S\right) = v(S \cup \{d\}) - v(S)\]
The Shapley value is a weighted average of these marginal contributions across all subsets $S \subseteq{\mathcal{D} - \{d\}}$ excluding agent $d$,

\[ \varphi_d(v) = \sum_{S \subseteq \mathcal{D} - \{d\}} \frac{1}{D {D - 1 \choose \left|S\right|}} C(d \vert S) \tag{1}\] The summation is over all subsets that don’t include agent $\{d\}$ (if it had included agent $d$, then the definition of the contribution $C(d \vert S)$ of $d$ to $S$ wouldn’t make sense).

Exercise: Express the following in terms of $v$. Which of these could appear in the definition of the Shapley value?
- $C(1 \vert 2)$
- $C(1 \vert \emptyset)$
- $C(1 \vert 1)$
- $C(3 \vert \{1, 2\})$
- $C(\{1, 2\} \vert 3)$
The weights $1/(D {D - 1 \choose \left|S\right|})$ ensure that the weights sum to 1, making $\varphi_{d}(v)$ a proper weighted average of marginal contributions. For any coalition of size $s$, there are ${D - 1 \choose s}$ subsets with that size, so the total weight is,

\[ \sum_{s=0}^{D-1} \binom{D-1}{s} \frac{1}{D\binom{D-1}{s}} = \sum_{s=0}^{D-1} \frac{1}{D} = 1. \]
The Shapley value is the unique solution satisfying the axioms below, giving this approach a principled justification,
- Efficiency. Shapley values sum to the total profit: \[ \sum_{d = 1}^{D} \varphi_{d}(v) = v\left(\mathcal{D}\right) - v\left(\emptyset\right) \] All profit is distributed and nothing is left over.
- Monotonicity. If agent $d$ contributes at least as much in game $v_1$ as in $v_2$ for every coalition¹, then $\varphi_d(v_1) \geq \varphi_d(v_2)$.
- Symmetry. Equal contributors receive equal credit. That is, if $v(S \cup \{d\}) = v(S \cup \{d'\})$ for every $S$, then $\varphi_d(v) = \varphi_{d'}(v)$.
- Missingness. An agent that never contributes receives nothing. That is, if $v\left(S \cup \{d\}\right) = v\left(S\right)$ for every $S$, then $\varphi_d(v) = 0$.

Machine Learning Analogy

The key insight is that we can map this game theoretic setup to the local feature attribution problem:
- “profit $v(\mathcal{D})$” $\to$ “prediction $f(x)$”
- “agent $d$” $\to$ “feature $d$.”
- “team $S$” $\to$ “subset of features $S \subset \mathcal{D}$.”
Instead of distributing profit across agents, we attribute a prediction $f(x)$ across features. We denote this attribution $\varphi_d(f, x)$
For this to work, we need to define $C(d \vert S)$ as the change in prediction feature $d$ is included vs. removed. Different answers to this feature removal question lead to different Shap values. But once $C$ is defined, we can substitute it into Equation 3 to get a local feature attribution for sample $x$.

Exercise: Pick one of the four axioms for game theoretic SHAP. What does it imply about $\varphi_d(f, x)$?

Deterministic Feature Removal

There are three common approaches to feature removal: baseline, marginal, and conditional. We’ll review each in turn.
Let $x'$ denote a baseline value. For example, $\mathbf{0}\in \reals^{D}$, or the sample mean $\bar{x} \in \reals^{D}$. Define \[ v(S) = f(x_{S}, x'_{\bar{S}}) \\ \] where $x_{S}$ and $x_{\bar{S}}$ index the coordinates of $x$ included and excluded from $S$. This uses the real feature values from sample $x$ for coordinates in $S$ and substitutes the baseline $x'$ elsewhere.

Exercise: Suppose that $x$ is an image and that $x'$ is the all zeros image. What would $\left(x_{S}, x'_{\bar{S}}\right)$ look like?
The marginal contributions become,

\[ \begin{align*} C(d \vert S) &= v\left(S \cup \{d\}\right) - v(S)\\n &= f\left(x_{S \cup \{d\}}, x'_{\overline{S \cup \{d\}}}\right) - f\left(x_{S}, x'_{\bar{S}}\right) \end{align*} \] This is the change in prediction when feature $d$ is included (top) vs. replaced by the baseline (bottom).

The figure below represents $C\left(2 \vert 1\right)$ under baseline removal.
The downside of this approach is that it depends on the choice of baseline $x'$, and there is no obvious principled way to choose it.

Sampling-based Feature Removal

Both the marginal and conditional approaches replace the deterministic baseline with an expectation over randomly sampled coordinates. Let $X_{\bar{S}}$ be a random vector of the features not in $S$, drawn from the training distribution. The marginal approach defines, \[ v(S) = \Esubarg{p(X_{\bar{S}})}{f(x_{S}, X_{\bar{S}})} \tag{2}\]
In practice, this expectation can be approximated by the training data, $x_1, \dots, x_{N}$, \[ v(S) = \frac{1}{N} \sum_{i = 1}^{N} f(x_S, x_{i,\bar{S}}). \]
To understand this definition geometrically, consider the figure below, which represents $C(2 | 1)$ under marginal removal.
- Left panel. The curved purple lines are contours of $f$ — think of the surface coming out of the page. The large dot at the intersection of the dashed lines is the observation $x = (x_1, x_2)$ whose prediction we want to explain. The smaller dots are training observations $x_i = (x_{i1}, x_{i2})$, and the blue dashed lines connect each $x_i$ to the point $(x_1, x_{i2})$, where the first coordinate has been replaced by $x_1$. These quantities appear directly in $C(2 \mid 1)$:
\[ C(2 \mid 1) = v(\{1,2\}) - v(\{1\}) = f(x) - \mathbb{E}[f(x_1, X_2)] \approx f(x) - \frac{1}{N} \sum_{i} f(x_1, x_{i2}) \]
- Right panel. The decreasing curve is a cross-section of $f$ along that vertical dashed line (i.e., holding $x_1$ fixed). The pink dashed line marks $x_2$ for the point we want to explain, and the solid horizontal line has $y$-value equal to $f(x)$. The blue dashed lines show the $x_{i2}$ from training samples, and the blue histogram summarizes the resulting values $f(x_1, x_{i2})$. The gap between $f(x)$ and the center of the histogram is $C(2 \mid 1)$.
Exercise: What would $C\left(1 \vert 2\right)$ look like in this figure? What about $C\left(1 \vert \emptyset\right)$?

Exercise: Would you expect $\varphi_1(f, x)$ to be larger or smaller than $\varphi_2(f, x)$ for the $f$ and $x$ shown? Explain your reasoning.
A potential downside of this approach is that, if features are correlated, the marginal approach may evaluate $f$ at unrealistic input combinations. This is because $X_{\bar{S}}$ is sampled independently from the observed $x_{S}$.

Exercise: Modify the previous figure to illustrate this extrapolation issue.
The conditional approach addresses this by drawing $X_{\bar{S}}$ from the distribution conditioned on the observed $x_{S}$,

\[ v(S) = \Esubarg{p\left(X_{\bar{S}} \vert X_{S} = x_{S}\right)}{f(x_{S}, X_{\bar{S}})}. \]
Unlike the marginal approach, there is no generic estimator for this conditional expectation. Under some assumptions (e.g., multivariate gaussianity), it might be available in closed form. Alternatively, one approach is to train a surrogate model to approximate and sample from the required conditional distributions.

Causal Intervention Perspective

Some researchers (Janzing et al. 2020) have argued that point (4) in the previous section is not actually a problem, and that, from a causality perspective, the marginal approach computes the more meaningful expectation.
To see this, we distinguish between observed data $\tilde{x}$ and model inputs $x$. An intervention conditional distribution sets the input coordinates ${S}$ to $x_{S}$ while leaving $X_{\bar{S}}$ undisturbed, written as \[\Earg{f(X_{S}, X_{\bar{S}}) \vert \text{do}(X_{S} = x_{S})}\] This differs from ordinary conditioning, which asks “given $X_{S} = x_{S}$, what do we expect $X_{\bar{S}}$ to look like?” In contrast, intervening asks “What if we force $X_S = x_S$ and leave everything else untouched?”

Exercise: Compare and contrast interventional conditional distributions with CP profiles.
Conditioning on $X_S = x_S$ can change the distribution of $X_{\bar{S}}$ in a way that creates nonzero attributions for any feature $d$ that is correlated with $X_S$, even even if they aren’t in the model. Interventions avoid this because they leave $X_{\bar{S}}$ unchanged.
To make this concrete, suppose $X_1$ is a thermometer measurement, $X_2$ is the actual temperature, and $f(X_1, X_2) = X_2$ only depends on the actual temperature. Since the two variables are correlated, $\Esubarg{p(X_2 \vert X_1 = x_1)}{X_2} \neq \Earg{X_2}$, so $C(1 \vert \emptyset) \neq 0$ and we get a nonzero attribution for $X_1$, even though it isn’t in the model².
It turns out that \[ \Esubarg{p(X_{\bar{S}})}{f(X_{S}, X_{\bar{S}}) \vert \text{do}(X_{S} = x_{S})} = \Esubarg{p(X_{\bar{S}})}{f(x_{S}, X_{\bar{S}})} \] The right hand side is exactly the marginal expectation from Equation 2. So, the marginal approach matches a causally meaningful quantity that reflects what it means to remove a feature’s influence.

Computation

Challenge

Let’s first see why naive SHAP computation is slow. Recall the key quantities,
- Feature $d$’s attribution, \[ \varphi_d(f, x) = \sum_{S \subseteq \mathcal{D} - \{d\}} \frac{1}{D {D - 1 \choose \left|S\right|}} C(d \mid S) \tag{3}\] where $C(d \mid S) = v(S \cup \{d\}) - v(S)$ is the marginal contribution of $d$ to coalition $S$. Both $C$ and $v$ depend on $f$ and $x$, but we suppress this from our notation. We also write $\varphi_d(x)$ instead of $\varphi_d(f, x)$, since $f$ is fixed in these notes.
- The value function $v(S)$ evaluates $f$ on $x$ after “removing” features outside of $S$ removed. We covered three removal strategies:
  - Baseline $v(S) = f(x_S, x'_{\bar{S}})$
  - Marginal³: $v(S) = \mathbb{E}_{p(X_{\bar{S}})}[f(x_S, X_{\bar{S}})] \approx \frac{1}{B}\sum_{i=1}^{B} f(x_S, x_{i,\bar{S}})$
  - Conditional: $v(S) = \mathbb{E}_{p(X_{\bar{S}} \mid X_S = x_S)}[f(x_S, X_{\bar{S}})]$
To compute all $D$ attributions for one sample, we need $v(S)$ for every $S \subseteq \mathcal{D}$. Using marginal feature removal, each $v(S)$ costs $B$ evaluations of $f$. Let’s try a Fermi calculation. Take $N = 100, D = 10, B = 100$. Attribution across all samples and features requires \[10^2 \text{ samples} \times 2^{10} \text{ subsets} \times 10^2 \text{ calls per subset} \approx 10^7 \text{ evaluations}.\] Supposing 1ms per call, this takes $\approx 3$ hours (and this is only 100 samples with 10 features). More generally, marginal removal costs $NB\,2^D$ evaluations.

Exercise: What are the relative costs of marginal vs. baseline feature removal?

Sampling Perspective

The weights in Equation 3 define a probability distribution, \[ \mathbb{P}(S) = \frac{1}{D \binom{D-1}{|S|}}, \] because they are nonnegative and sum to one. This lets us rewrite $\varphi_d(x)$ as an expectation, \[\begin{align} \label{eq-distribution} \varphi_d(x) &= \sum_{S \subseteq \mathcal{D} \setminus \{d\}} \mathbb{P}(S)\, C(d \mid S) = \mathbb{E}_{\mathbb{P}(S)}[C(d \mid S)]. \end{align}\] i.e., the Shapley value is the expected contribution of feature $d$ to a random coalition.

Instead of enumerating all $2^{D-1}$ subsets to get the exact expectation, we sample $S_1, \dots, S_M \sim \mathbb{P}(S)$ and estimate \[ \hat{\varphi}_d(x) = \frac{1}{M}\sum_{m=1}^{M} C(d \mid S_m). \tag{4}\] This is unbiased ($\mathbb{E}[\hat{\varphi}_d(x)] = \varphi_d(x)$) and increasing $M$ improves the approximation.

Exercise: Suppose $\mathcal{D} = \{1, \dots, 4\}$ and $d = 1$. What are the $\mathbb{P}(S)$ for the sets in Equation 3?
To draw $S_{m} \sim \mathbb{P}(S)$ in practice, we use random permutations. Draw $\pi : \{1, \dots, D\} \to \{1, \dots, D\}$ uniformly at random. For example, for 10 features we could use,
```
pi <- sample(1:10)
pi
```
```
 [1]  9  1  6  5  3  4  7 10  8  2
```
Let $\text{Pre}^{d}\left(\pi\right)$ the features appearing before $d$ in $\pi$. E.g., for $d = 4$,
```
 pi[seq_len(match(4, pi) - 1)]
```
```
[1] 9 1 6 5 3
```
It’s not obvious, but $\text{Pre}^d\left(\pi\right)$ gives a sample from $\mathbb{P}\left(S\right)$.
To summarize, the procedure is,
```
for m = 1, ..., M:
  sample a random permutation π_m over D features
  set S_m = Pre^d(π_m)
  compute C(d | S_m)
return average of C(d | S_m) over m
```
This is often called IME (“Interactions-based Method for Explanation”). Larger $M$ gives better approximations but at higher computing cost.
Adaptive sampling. The terms $C\left(d \vert S_{m}\right)$ are random variables with mean $\varphi_{d}\left(x\right)$.

For features that are irrelevant, the contributions are tightly clustered near zero, so even small values of $M$ are enough. Inspired by this observation, adaptive sampling tracks the variance of $C(d \mid S_m)$ across features and allocates more samples to high-variance features.
Antithetic sampling. Recall that \[ \text{Var}(Z_1 + Z_2) = \text{Var}(Z_1) + \text{Var}(Z_2) + 2\text{Cov}(Z_1, Z_2). \] so weakly correlated variables have smaller variances when added.

Exercise: If $\text{Cor}(Z_1, Z_2) = -1$ and $Z_1, Z_2$ have the same variance, what is $\text{Var}(Z_1 + Z_2)$?

More generally, \[\text{Var}\!\left(\sum_{m=1}^M Z_m\right) = \sum_m \text{Var}(Z_m) + \sum_{m \neq m'} \text{Cov}(Z_m, Z_{m'}),\] so again smaller pairwise correlations reduce variance.
Applying this to Equation 4, we can reduce the correlation between $C\left(d \vert S_m\right)$ and $C\left(d \vert S_{m'}\right)$ by pairing each $S_m$ with its complement,
- Sample $S_m \sim \mathbb{P}(S)$
- Set $S_{m'} = \bar{S}_{m}$ the complement of $S_m$ in $\{1, \dots, D\}$.
Despite its simplicity, this yields state-of-the-art variance reductions (Mitchell et al. 2022).

Weighted Least Squares

A classical result in game theory (Charnes et al. 1988) shows that SHAP values can be characterized through a weighted least squares problem, \[ \varphi\left(x\right) = \arg \min_{\beta \in \reals^{D + 1}} \sum_{S \subset \mathcal{D}} w(S) \left(v(S) - \beta_0 - \sum_{d \in S}\beta_{d}\right)^2 \tag{5}\] where $w(S) = \left[\binom{D}{|S|} |S|(D - |S|)\right]^{-1}$ and $\varphi(x) = (\varphi_1(x), \dots, \varphi_D(x))^\top$.
This characterization still enumerates $2^D$ subsets, but it suggests an approximation: sample random subsets $S_1, \dots, S_M$ and solve the weighted regression on the reduced design. This is “KernelSHAP.”
Let $z_{md} = \mathbf{1}[d \in s_m]$ be the associated design matrix. The regression suggests $v(\emptyset) \approx \beta_0$ and \[ v(S_m) - v(\emptyset) \approx \sum_{d = 1}^{D} z_{md}\beta_d. \tag{6}\] The efficiency constraint becomes, \[ v(\mathcal{D}) - v(\emptyset) \approx \sum_{d = 1}^{D} \beta_d =: \Delta\\ \] We can improve the quality of our approximation by enforcing these constraints, substituting $v(\emptyset) = \Esubarg{p(x)}{f(X)} \approx \frac{1}{N}\sum_{i = 1}^{N} f(x_i)$ and $v(\mathcal{D}) = f(x)$.

Exercise: In fact, the regression objective ensures that $v\left(\emptyset\right) = \beta_0$ and $v(\mathcal{D}) - v(\emptyset) = \Delta$ hold exactly. Can you see why? (Hint: Consider the weights)
To see this concretely, let’s try a code implementation (validated against the shap package below). First, initialize the binary mask matrix $z_{md}$, kernel weights, and feature-removed inputs:
```
mask_matrix = np.zeros((M, D))
kernel_weights = np.zeros(M)
synth_data = np.tile(data, (M, 1))
```

Next we draw random subsets $S_m$ and the corresponding feature-removed data $(x_{S_m}, x_{i,\bar{S}_m})$:

for i in range(M):
     mask = np.random.choice([0, 1], size=D)
     mask_matrix[i] = mask
     synth_data[i * N:(i + 1) * N, mask == 1] = x[0, mask == 1] # explaining x[0]
     kernel_weights[i] = 1 / (comb(self.D, mask.sum()) * mask.sum() * (self.D - mask.sum()))

We evaluate $f$ on the feature-removed data. Here f_bar estimates $v(\emptyset)$.

 f_bar = np.mean(model(data))
 f = model(synth_data)
 v = np.mean(f.reshape(M, N, -1), axis=1)

Exercise: In the line synth_data[i, ...], how many entries of synth_data get modified? What is the interpretation?

Finally, we run the weighted regression. We incorporate efficiency constraints from point (3). Since $\beta_D = \Delta - \sum_{d= 1}^{D - 1} \beta_d$, we can substitute into Equation 6, \[ v(S_m) - v(\emptyset) - z_{mD}\Delta \approx \sum_{d = 1}^{D - 1}\left(z_{md}- z_{mD}\right)\beta_{d} \] The left hand side is y below and the right hand side is X.
```
y = (v - f_bar).flatten() - mask_matrix[:, -1] * (model(x) - f_bar)
X = mask_matrix - mask_matrix[:, -1][:, None]
lm = LinearRegression(fit_intercept=False).fit(X, y, sample_weight=kernel_weights)

# final shapley values
phi = np.zeros(D)
phi[:-1] = lm.coef_[:-1]
phi[-1] = (model(x) - f_bar) - np.sum(lm.coef_)
```

Amortization

So far, we’ve only explained a single observation $x$. In practice, we want attributions for many $\{x_i\}_{i =1}^{N}$. An important insight is that attributions across observations are often related, so we can pool our estimates across all samples. This lets us tolerate rougher estimates for each individual sample.

Notice that many samples have similar attributions.
Amortization frames parameter estimation as supervised learning. For concreteness, consider noisy KernelSHAP estimates using small $M$, \[ \hat{\varphi}_d^{M}(x_i) = \arg \min_{\beta \in \reals^{D + 1}} \sum_{m = 1}^{M} w(S_m) \left(v(S_m) - \beta_0 - \sum_{d \in S_m}\beta_{d}\right)^{2} \] Then train a model $a_\theta : \reals^D \to \reals^D$ (e.g. linear or neural network) to predict these estimates, \[ \hat{\theta} = \arg\min_{\theta \in \Theta} \sum_{i=1}^{N} \left\|\hat{\varphi}^M(x_i) - a_\theta(x_i)\right\|^2. \] This is now like a supervised learning problem where the “labels” are $\hat{\varphi}^{M}(x_i)$ and the input features are $x_i$.
Once $\hat{\theta}$ has been found, we can very quickly get attributions for a new sample $x^\ast$ by just calling $a_{\hat{\theta}}(x^\ast)$. We don’t even need to sample any subsets $S_m$ or evaluate any $v(S_m)$. Moreover, the learned model can smooth over the noise in the per-sample estimates, leading to $a_{\hat{\theta}}(x^*)$ that approximates $\varphi(x^*)$ more accurately than $\hat{\varphi}^M(x^*)$.
Some caveats. Training $a_{\theta}$ can itself be computationally intensive. Also, if the labels $\hat{\varphi}^M(x_i)$ are too noisy, the learned mapping may be inaccurate. Nonetheless, early research suggests that even moderate $M$ and practical $a_{\theta}$ can learn $\varphi(x)$ with high fidelity.

Model-Specific Versions

Within specific model classes, we can derive closed-form SHAP formulas that bypass enumeration over subsets. For example, it is straightforward to show that linear models $f(x) = \beta_0 + \sum_{d} \beta_{d}x_{d}$ with baseline removal satisfy, \[ \varphi_{d}(x) = \beta_{d}\left(x_{d} - x_{d}'\right) \] The attribution only depends on feature $d$’s deviation from the baseline and coefficient $\beta_d$. A similar derivation shows that for marginal feature removal, \[ \varphi_{d}(x) = \beta_{d}\left(x_{d} - \Esubarg{p(x)}{X_{d}}\right) \] so using the marginal removal is the same as using the sample mean as a baseline.
There are two other important model-specific methods,
- TreeSHAP: This gives exact SHAP values tree-based model (CART, bagging, boosting, random forest). It’s available as TreeExplainer in the shap python package and the treeshap R package.
- DeepSHAP: This is an approximation for deep learning models. It’s not exact, but has lower variance than model-agnostic approaches. It is available through the DeepExplainer class in shap and the DeepSHAP class in the innsight R package.

Case Studies

KernelSHAP Implementation

This example shows that the KernelSHAP sketch above indeed agrees with the official implementation. We’ve copied our code into the explain method in the class below. shap_values simply loops explain over many samples.

import numpy as np
from math import comb
from sklearn.linear_model import LinearRegression

class SimpleKernelExplainer:
    def __init__(self, model, data):
        self.model = model
        self.data = data
        self.N, self.D = data.shape
        self.f_bar = np.mean(model(data))

    def shap_values(self, X, M):
        explanations = []
        for i in range(X.shape[0]):
            explanations.append(self.explain(X[i:(i + 1)], M))
        return np.array(explanations)

    def explain(self, x, M):
        # initialize data structures
        mask_matrix = np.zeros((M, self.D))
        kernel_weights = np.zeros(M)
        synth_data = np.tile(self.data, (M, 1))

        # sample random subsets
        for i in range(M):
            mask = np.random.choice([0, 1], size=self.D)
            mask_matrix[i] = mask
            synth_data[i * self.N:(i + 1) * self.N, mask == 1] = x[0, mask == 1]

            # edge 0 and self.D handled by constraints
            if 0 < mask.sum() < self.D:
                kernel_weights[i] = 1 / (comb(self.D, mask.sum()) * mask.sum() * (self.D - mask.sum()))

        # model predictions on the synthetic data
        model_out = self.model(synth_data)
        v = np.mean(model_out.reshape(M, self.N, -1), axis=1)

        # weighted linear regression on adjusted f
        X = mask_matrix - mask_matrix[:, -1][:, None]
        y = (v - self.f_bar).flatten() - mask_matrix[:, -1] * (self.model(x) - self.f_bar)
        lm = LinearRegression(fit_intercept=False).fit(X, y, sample_weight=kernel_weights)

        # return shapley approximation
        phi = np.zeros(self.D)
        phi[:-1] = lm.coef_[:-1]
        phi[-1] = (self.model(x) - self.f_bar) - np.sum(lm.coef_)
        return phi

We’ll test this implementation using a random forest model trained on the make_friedman1 dataset (Friedman 1991).

from sklearn.ensemble import RandomForestRegressor
from sklearn.datasets import make_friedman1

X, y = make_friedman1(n_samples=100, n_features=5, noise=0.1)
model = RandomForestRegressor(n_estimators=100)
model.fit(X, y)

RandomForestRegressor()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

We run both our simple and the official explainer in the block below.

import shap

explainer = SimpleKernelExplainer(model.predict, X)
official_explainer = shap.KernelExplainer(model.predict, X)

shap_values = explainer.shap_values(X, M=1000)
official_shap_values = official_explainer.shap_values(X, nsamples=1000, silent=True)

There’s randomness in the estimates, but they’re close.

print("Package: ", np.round(official_shap_values[:5], 4))

Package:  [[ 2.2775  1.1336 -0.0345  2.5457  0.5317]
 [ 2.9126  1.2404  0.7779 -2.7753 -0.1276]
 [ 2.207   1.2553  0.1723  5.6079  0.2038]
 [ 1.2856  1.9748  0.5054  2.6508 -0.3298]
 [ 1.5955  0.9755 -0.0149 -2.1233  0.7739]]

print("Ours: ", np.round(shap_values[:5], 4))

Ours:  [[ 2.2861e+00  1.1283e+00 -2.7200e-02  2.5376e+00  5.2930e-01]
 [ 2.9760e+00  1.1748e+00  7.5550e-01 -2.8075e+00 -7.0800e-02]
 [ 2.2006e+00  1.2681e+00  1.8490e-01  5.5471e+00  2.4550e-01]
 [ 1.2249e+00  1.9926e+00  4.8680e-01  2.6801e+00 -2.9760e-01]
 [ 1.5935e+00  9.7440e-01  3.2000e-03 -2.1365e+00  7.7210e-01]]

Income Data

Let’s apply the shap python package to identify important variables in the adult dataset. Each $x_i$ is a survey response for person $i$, and $y_i$ indicates whether they make more than $50K/year.

import shap
import matplotlib.pyplot as plt
plt.rcParams['figure.autolayout'] = True

X, y = shap.datasets.adult()
X, y = X.iloc[:2000], y[:2000]   # subsample for speed
X

       Age  Workclass  Education-Num  ...  Capital Loss  Hours per week  Country
0     39.0          7           13.0  ...           0.0            40.0       39
1     50.0          6           13.0  ...           0.0            13.0       39
2     38.0          4            9.0  ...           0.0            40.0       39
3     53.0          4            7.0  ...           0.0            40.0       39
4     28.0          4           13.0  ...           0.0            40.0        5
...    ...        ...            ...  ...           ...             ...      ...
1995  44.0          4           10.0  ...           0.0            50.0       39
1996  49.0          4            9.0  ...           0.0            60.0       39
1997  75.0          6           14.0  ...           0.0            50.0       39
1998  37.0          4           13.0  ...           0.0            55.0       39
1999  51.0          7            9.0  ...           0.0            38.0       39

[2000 rows x 12 columns]

We train a random forest model to these data using the sklearn package.

from sklearn.ensemble import RandomForestClassifier
model = RandomForestClassifier(n_estimators=100)
model.fit(X, y)

RandomForestClassifier()

Here, we’re explaining the first 50 samples $x_1, \dots, x_{50}$ using the marginal feature removal approach, as implemented by KernelExplainer. The $N = 2000$ rows in X are used in estimating the marginal expectation Equation 2. Notice that the explainer only needs access to the anonymous function predictor – it doesn’t require any knowledge of the type of model implemented within it.
```
X_explain = X.iloc[:50]
predictor = lambda X: model.predict_proba(X)[:, 1]

#explainer = shap.KernelExplainer(predictor, X) # if you want to run accurate version
explainer = shap.KernelExplainer(predictor, X.sample(100)) # if you want to run the fast version
sv = explainer.shap_values(X_explain)
```
Exercise: What is the dimension of sv?

A waterfall plot shows $\varphi_d(f, x)$ for all $d$ features. By the efficiency axiom, the sum of these bars is equal to the predicted response.

exp_single = shap.Explanation(
 values = sv[0],
 base_values = explainer.expected_value,
 data = X_explain.iloc[0].values,
 feature_names = X_explain.columns.tolist(),
)
shap.plots.waterfall(exp_single)

Here is an alternative visualization well-suited for visualizing attributions for multiple samples (lines) simultaneously. This visualization is helpful for identifying clusters of samples with similar feature attributions.
```
shap.decision_plot(explainer.expected_value, sv, X_explain)
```
Exercise: Imagine explaining this visualization to a non-data scientist. Describe each component without jargon and summarize the main takeaways within the salary prediction context.

Drosophila Development

Problem Formulation

Enhancers are DNA sequences that coordinate downstream gene expression. Their activity depends on genome state, like transcription factor binding and chromatin state.
Understanding enhancers means understanding how genotype becomes phenotype, a fundamental problem in biology. To make the problem tractable, we usually study it in simple organisms, like the fruit fly Drosophila melanogaster, often during embryo development.
(Basu et al. 2018) asked a precise question: given the genomic state of a fruit fly embryo at a fixed developmental stage, can we predict whether a given enhancer is active? They used genome-wide measurements from blastoderm (stage 5) embryos:
- DNA occupancy for 23 transcription factors
- Activity for 13 chromatin markers
- Binary enhancer activity labels
This is a binary classification problem. Each observation is a genomic region with associated regulatory features.
- Enhancer status: $y \in \{0,1\}$
- Predictors: $x = (x_1, \ldots, x_D)$ transcription factor binding intensities and chromatin signals
- Model: $f(x_i) = \mathbb{P}(y_i = 1 \mid x_i)$
SHAP quantifies each feature’s contribution to $f(x_i)$.

Setup

The first three packages below train the boosting model. The rest are for visualization and SHAP explanation.

source("https://github.com/krisrs1128/stat479_notes/raw/refs/heads/master/activities/05-helpers.R")
library(mlr3)
library(mlr3learners)
library(mlr3tuning)
library(patchwork)
library(shapviz)
library(tidyverse)
theme_set(theme_classic())
set.seed(20251227)
options(mlr3.store_backends = TRUE)

This block downloads the data. Enhancer is stored in the first column y.

f <- tempfile()
download.file("https://zenodo.org/records/18371236/files/enhancer.Rdata?download=1", f)
load(f)

enhancer <- bind_cols(y = as.factor(Y), X)
enhancer

# A tibble: 7,809 × 81
   y     H3_c12 H3_c14a H3_c14c H3_c8 H3K18ac_c12 H3K18ac_c14a H3K18ac_c14c
   <fct>  <dbl>   <dbl>   <dbl> <dbl>       <dbl>        <dbl>        <dbl>
 1 0          6      10      10     8           5            7            8
 2 0          6       7       6     5          19           10           11
 3 0         14       9      14    13          11            9           19
 4 0         13      17      13    13          19           14           25
 5 1          6       6       6    10          13           23           28
 6 1         10       9       9    11          12           23           28
 7 0          6       8       8    10          11           13           24
 8 0          8      10      16     6           7           19           16
 9 0         16      20      11    10          34            8           15
10 0         16      16      16    11          13           20           22
# ℹ 7,799 more rows
# ℹ 73 more variables: H3K18ac_c8 <dbl>, H3K27ac_c12 <dbl>, H3K27ac_c14a <dbl>,
#   H3K27ac_c14c <dbl>, H3K27ac_c8 <dbl>, H3K27me3_c12 <dbl>,
#   H3K27me3_c14a <dbl>, H3K27me3_c14c <dbl>, H3K36me3_c12 <dbl>,
#   H3K36me3_c14a <dbl>, H3K36me3_c14c <dbl>, H3K4me1_c12 <dbl>,
#   H3K4me1_c14a <dbl>, H3K4me1_c14c <dbl>, H3K4me1_c8 <dbl>,
#   H3K4me3_c12 <dbl>, H3K4me3_c14a <dbl>, H3K4me3_c14c <dbl>, …

Model Training

We train a boosting classifier on the enhancer data. The tuning grid considers the number of trees nrounds and learning rate eta, selecting hyperparameters using three fold CV.

task <- as_task_classif(enhancer, target = "y", id = "enhancer", positive = "1")
learner <- lrn("classif.xgboost", nrounds = to_tune(100, 200), eta = to_tune(1e-3, 0.1), predict_type = "prob")

instance <- ti(
    task = task,
    learner = learner,
    resampling = rsmp("cv", folds = 3),
    measures = msr("classif.auc"),
    terminator = trm("none")
)

We run the hyperparameter tuning next.

tuner <- tnr("grid_search", resolution = 5, batch_size = 5)
tuner$optimize(instance)

       eta nrounds learner_param_vals  x_domain classif.auc
     <num>   <int>             <list>    <list>       <num>
1: 0.02575     100          <list[6]> <list[2]>    0.825632

Let’s refit the model with the optimized hyperparameters.

learner$param_set$values <- instance$result_learner_param_vals
task <- as_task_classif(enhancer, target = "y", id = "enhancer", positive = "1")
learner$train(task)

The histogram of predicted probabilities shows that the model behaves as expected: inactive enhancers (class 0) clusters near zero probability, and active enhancers (class 1) is more spread to high probability regions. Though, there is still noticeable overlap.

tibble(p_hat = learner$predict_newdata(enhancer)$prob[, 1], y = enhancer$y) |>
    ggplot() +
    geom_histogram(aes(p_hat, fill = y)) +
    facet_grid(y ~ ., scales = "free_y")

Interpretation Setup

Entry $(i, d)$ of shp$S is the SHAP attribution $\varphi_d(x_i)$, feature $d$’s contribution to $f(x_i)$.

X_mat <- as.matrix(X)
shp <- shapviz(learner$model, X_pred = X_mat)
shp$S[1:5, 1:8]

             H3_c12      H3_c14a      H3_c14c       H3_c8   H3K18ac_c12
VT0002  0.068839401 -0.034920834 -0.008018139 -0.07555990  0.0327704698
VT0003  0.064319402 -0.004164241  0.007724053 -0.03716381 -0.0258688685
VT0004  0.083743654 -0.037574559 -0.007068704 -0.07883391  0.0002055851
VT0005  0.055233277  0.009904596 -0.002665190 -0.03974437 -0.0256465226
VT0006 -0.004646949 -0.001682470 -0.053047363 -0.02233944 -0.0076694093
       H3K18ac_c14a  H3K18ac_c14c  H3K18ac_c8
VT0002  0.041389603  0.0031723969 0.011215613
VT0003  0.016818035  0.0003347002 0.008839345
VT0004  0.041915886 -0.0014176624 0.001084273
VT0005 -0.001247977  0.0006725105 0.016408261
VT0006 -0.009059546 -0.0014770469 0.007148674

The efficiency axiom holds, though on a log-odds scale: $\sum_d \varphi_d(x_i) + \varphi_0 = \text{logit}^{-1}\mathbb{P}(y \vert x_i)$. Most feature importance measures lack this sum-to-prediction property; this is something that helps SHAP stand out.

y_hat <- predict(learner$model, X_mat, outputmargin = TRUE)

tibble(y_hat = y_hat, sum_phi = rowSums(shp$S) + shp$baseline) |>
    ggplot() +
    geom_point(aes(y_hat, sum_phi)) +
    labs(
        x = "Model Prediction (Log-Odds)",
        y = "Sum of SHAP Attributions"
    )

Variable-Level Plots

The beeswarm plot ranks features using $\left|\varphi_d(x_i)\right|$. Across many samples, the wt_ZLD, input_c14c, and sna1 features stand out. For most features, the attributions have consistent signs across samples, but the tails show that the sign can flip, perhaps when the feature has an unusual value/context.

sv_importance(shp, kind = "bee")

The dependence plot for sna1 shows that high sna1 TF binding leads to positive attribution.

sv_dependence(shp, "sna1")

This is consistent with the data, where large sna1 increases enhancer activity. The sna1 gene is generally a transcriptional repressor, but in these data, it seems that when sna1 is present at a location, the enhancer is more likely to be active.

tibble(p_hat = predict(learner$model, X_mat, outputmargin = TRUE), y = enhancer$y, x= enhancer$sna1) |>
    ggplot() +
    geom_point(aes(x, p_hat, col = y)) +
    labs(x = "Value of SNA1 Feature", y = "Log-Odds", col = "True Class")

Sample-Level Plots

Alternatively, waterfall plots break a single prediction down into feature contributions. Each bar is one feature, and tracking the endpoints across bars shows how the model arrives at the output by weighing evidence across many features.

map(seq_len(6), \(i) sv_waterfall(shp, max_display = 8, row = i)) |>
    reduce(`+`)

Waterfall plots don’t scale well to many samples. A stacked bar alternative solves this issue – each vertical bar below is one sample, colored by SHAP attributions for the 11 most improtant features. Samples are sorted so that those with similar attribution profiles are side-by-side. The “x” marks give the actual log-odds predictions, and features above or below indicate positive or negative attributions, respectively.

The block below calls our helper functions (in 05-helpers.R) to perform this ordering.

var_order <- order(abs(colSums(shp$S)), decreasing = TRUE)
important_vars <- colnames(shp$S)[var_order] |>
    head(11)

ix <- sample(1:nrow(enhancer), 500)
shap_bars <- shap_to_long(shp$S[ix, ], y_hat[ix], shp$baseline, important_vars) |>
    add_waterfall_offsets()

This view shows that samples with nearly identical predictions can have different reasons for those predictions. For example, kr1 drives high predictions for a cluster of sequences near the left, but has barely any influence elsewhere.

ggplot(shap_bars) +
    geom_hline(yintercept = 0) +
    geom_rect(aes(xmin = xpos - 0.45, xmax = xpos + 0.45, ymin = ymin, ymax = ymax, fill = name, col = name)) +
    geom_point(
        data = distinct(shap_bars, xpos, pred),
        aes(xpos, pred),
        size = 3,
        shape = "x"
    ) +
    labs(y = "Prediction relative to E[f(x)]", x = "Samples", fill = "Feature", col = "Feature") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_fill_brewer(palette = "Set3", na.value = "#dcdcdc") +
    scale_color_brewer(palette = "Set3", na.value = "#dcdcdc") +
    theme(
        axis.ticks = element_blank(),
        axis.text.x = element_blank(),
    )

References

Basu, Sumanta, Karl Kumbier, James B. Brown, and Bin Yu. 2018. “Iterative Random Forests to Discover Predictive and Stable High-Order Interactions.” Proceedings of the National Academy of Sciences 115 (8): 1943–48. https://doi.org/10.1073/pnas.1711236115.

Charnes, A., B. Golany, M. Keane, and J. Rousseau. 1988. “Extremal Principle Solutions of Games in Characteristic Function Form: Core, Chebychev and Shapley Value Generalizations.” In Econometrics of Planning and Efficiency. Springer Netherlands. https://doi.org/10.1007/978-94-009-3677-5_7.

Friedman, Jerome H. 1991. “Multivariate Adaptive Regression Splines.” The Annals of Statistics 19 (1). https://doi.org/10.1214/aos/1176347963.

Janzing, Dominik, Lenon Minorics, and Patrick Bloebaum. 2020. Feature Relevance Quantification in Explainable AI: A Causal Problem. Proceedings of machine learning research, vol. 108: 2907–16. https://proceedings.mlr.press/v108/janzing20a.html.

Mitchell, Rory, Joshua Cooper, Eibe Frank, and Geoffrey Holmes. 2022. “Sampling Permutations for Shapley Value Estimation.” J. Mach. Learn. Res. 23 (1).

Ribeiro, Marco Tulio, Sameer Singh, and Carlos Guestrin. 2016. “‘Why Should i Trust You?’: Explaining the Predictions of Any Classifier.” Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, August, 1135–44. https://doi.org/10.1145/2939672.2939778.

Footnotes

That is, if $v_1\left(S \cup \{d\}\right) - v_1(S) \geq v_2\left(S \cup \{d\}\right) - v_2(S)$.↩︎
In full detail, $C(1 \vert \emptyset) = v(\{1\}) - v(\emptyset) = \Esubarg{p(X_2 \vert X_1 = x_1)}{X_2} - \Esubarg{p(X_1, X_2)}{X_2} \neq 0$ and $C(1 \vert 2) = v(\{1, 2\}) - v(\{2\}) = f(x_1, x_2) - \Esubarg{p(X_1 \vert X_2 = x_2)}{f(X_1, x_2)} = x_2 - x_2 = 0$ where in the last step we used the fact that $f(X_1, x_2) = x_2$ deterministically. Therefore the terms don’t cancel and we obtain a nonzero Shapley value.↩︎
Compared to the earlier discussion, we’re allowing for background subsamples B < N.↩︎

	n_estimators n_estimators: int, default=100 The number of trees in the forest. .. versionchanged:: 0.22 The default value of ``n_estimators`` changed from 10 to 100 in 0.22.	100
	criterion criterion: {"squared_error", "absolute_error", "friedman_mse", "poisson"}, default="squared_error" The function to measure the quality of a split. Supported criteria are "squared_error" for the mean squared error, which is equal to variance reduction as feature selection criterion and minimizes the L2 loss using the mean of each terminal node, "friedman_mse", which uses mean squared error with Friedman's improvement score for potential splits, "absolute_error" for the mean absolute error, which minimizes the L1 loss using the median of each terminal node, and "poisson" which uses reduction in Poisson deviance to find splits. Training using "absolute_error" is significantly slower than when using "squared_error". .. versionadded:: 0.18 Mean Absolute Error (MAE) criterion. .. versionadded:: 1.0 Poisson criterion.	'squared_error'
	max_depth max_depth: int, default=None The maximum depth of the tree. If None, then nodes are expanded until all leaves are pure or until all leaves contain less than min_samples_split samples.	None
	min_samples_split min_samples_split: int or float, default=2 The minimum number of samples required to split an internal node: - If int, then consider `min_samples_split` as the minimum number. - If float, then `min_samples_split` is a fraction and `ceil(min_samples_split * n_samples)` are the minimum number of samples for each split. .. versionchanged:: 0.18 Added float values for fractions.	2
	min_samples_leaf min_samples_leaf: int or float, default=1 The minimum number of samples required to be at a leaf node. A split point at any depth will only be considered if it leaves at least ``min_samples_leaf`` training samples in each of the left and right branches. This may have the effect of smoothing the model, especially in regression. - If int, then consider `min_samples_leaf` as the minimum number. - If float, then `min_samples_leaf` is a fraction and `ceil(min_samples_leaf * n_samples)` are the minimum number of samples for each node. .. versionchanged:: 0.18 Added float values for fractions.	1
	min_weight_fraction_leaf min_weight_fraction_leaf: float, default=0.0 The minimum weighted fraction of the sum total of weights (of all the input samples) required to be at a leaf node. Samples have equal weight when sample_weight is not provided.	0.0
	max_features max_features: {"sqrt", "log2", None}, int or float, default=1.0 The number of features to consider when looking for the best split: - If int, then consider `max_features` features at each split. - If float, then `max_features` is a fraction and `max(1, int(max_features * n_features_in_))` features are considered at each split. - If "sqrt", then `max_features=sqrt(n_features)`. - If "log2", then `max_features=log2(n_features)`. - If None or 1.0, then `max_features=n_features`. .. note:: The default of 1.0 is equivalent to bagged trees and more randomness can be achieved by setting smaller values, e.g. 0.3. .. versionchanged:: 1.1 The default of `max_features` changed from `"auto"` to 1.0. Note: the search for a split does not stop until at least one valid partition of the node samples is found, even if it requires to effectively inspect more than ``max_features`` features.	1.0
	max_leaf_nodes max_leaf_nodes: int, default=None Grow trees with ``max_leaf_nodes`` in best-first fashion. Best nodes are defined as relative reduction in impurity. If None then unlimited number of leaf nodes.	None
	min_impurity_decrease min_impurity_decrease: float, default=0.0 A node will be split if this split induces a decrease of the impurity greater than or equal to this value. The weighted impurity decrease equation is the following:: N_t / N * (impurity - N_t_R / N_t * right_impurity - N_t_L / N_t * left_impurity) where ``N`` is the total number of samples, ``N_t`` is the number of samples at the current node, ``N_t_L`` is the number of samples in the left child, and ``N_t_R`` is the number of samples in the right child. ``N``, ``N_t``, ``N_t_R`` and ``N_t_L`` all refer to the weighted sum, if ``sample_weight`` is passed. .. versionadded:: 0.19	0.0
	bootstrap bootstrap: bool, default=True Whether bootstrap samples are used when building trees. If False, the whole dataset is used to build each tree.	True
	oob_score oob_score: bool or callable, default=False Whether to use out-of-bag samples to estimate the generalization score. By default, :func:`~sklearn.metrics.r2_score` is used. Provide a callable with signature `metric(y_true, y_pred)` to use a custom metric. Only available if `bootstrap=True`. For an illustration of out-of-bag (OOB) error estimation, see the example :ref:`sphx_glr_auto_examples_ensemble_plot_ensemble_oob.py`.	False
	n_jobs n_jobs: int, default=None The number of jobs to run in parallel. :meth:`fit`, :meth:`predict`, :meth:`decision_path` and :meth:`apply` are all parallelized over the trees. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details.	None
	random_state random_state: int, RandomState instance or None, default=None Controls both the randomness of the bootstrapping of the samples used when building trees (if ``bootstrap=True``) and the sampling of the features to consider when looking for the best split at each node (if ``max_features < n_features``). See :term:`Glossary ` for details.	None
	verbose verbose: int, default=0 Controls the verbosity when fitting and predicting.	0
	warm_start warm_start: bool, default=False When set to ``True``, reuse the solution of the previous call to fit and add more estimators to the ensemble, otherwise, just fit a whole new forest. See :term:`Glossary ` and :ref:`tree_ensemble_warm_start` for details.	False
	ccp_alpha ccp_alpha: non-negative float, default=0.0 Complexity parameter used for Minimal Cost-Complexity Pruning. The subtree with the largest cost complexity that is smaller than ``ccp_alpha`` will be chosen. By default, no pruning is performed. See :ref:`minimal_cost_complexity_pruning` for details. See :ref:`sphx_glr_auto_examples_tree_plot_cost_complexity_pruning.py` for an example of such pruning. .. versionadded:: 0.22	0.0
	max_samples max_samples: int or float, default=None If bootstrap is True, the number of samples to draw from X to train each base estimator. - If None (default), then draw `X.shape[0]` samples. - If int, then draw `max_samples` samples. - If float, then draw `max(round(n_samples * max_samples), 1)` samples. Thus, `max_samples` should be in the interval `(0.0, 1.0]`. .. versionadded:: 0.22	None
	monotonic_cst monotonic_cst: array-like of int of shape (n_features), default=None Indicates the monotonicity constraint to enforce on each feature. - 1: monotonically increasing - 0: no constraint - -1: monotonically decreasing If monotonic_cst is None, no constraints are applied. Monotonicity constraints are not supported for: - multioutput regressions (i.e. when `n_outputs_ > 1`), - regressions trained on data with missing values. Read more in the :ref:`User Guide `. .. versionadded:: 1.4	None

	n_estimators n_estimators: int, default=100 The number of trees in the forest. .. versionchanged:: 0.22 The default value of ``n_estimators`` changed from 10 to 100 in 0.22.	100
	criterion criterion: {"gini", "entropy", "log_loss"}, default="gini" The function to measure the quality of a split. Supported criteria are "gini" for the Gini impurity and "log_loss" and "entropy" both for the Shannon information gain, see :ref:`tree_mathematical_formulation`. Note: This parameter is tree-specific.	'gini'
	max_depth max_depth: int, default=None The maximum depth of the tree. If None, then nodes are expanded until all leaves are pure or until all leaves contain less than min_samples_split samples.	None
	min_samples_split min_samples_split: int or float, default=2 The minimum number of samples required to split an internal node: - If int, then consider `min_samples_split` as the minimum number. - If float, then `min_samples_split` is a fraction and `ceil(min_samples_split * n_samples)` are the minimum number of samples for each split. .. versionchanged:: 0.18 Added float values for fractions.	2
	min_samples_leaf min_samples_leaf: int or float, default=1 The minimum number of samples required to be at a leaf node. A split point at any depth will only be considered if it leaves at least ``min_samples_leaf`` training samples in each of the left and right branches. This may have the effect of smoothing the model, especially in regression. - If int, then consider `min_samples_leaf` as the minimum number. - If float, then `min_samples_leaf` is a fraction and `ceil(min_samples_leaf * n_samples)` are the minimum number of samples for each node. .. versionchanged:: 0.18 Added float values for fractions.	1
	min_weight_fraction_leaf min_weight_fraction_leaf: float, default=0.0 The minimum weighted fraction of the sum total of weights (of all the input samples) required to be at a leaf node. Samples have equal weight when sample_weight is not provided.	0.0
	max_features max_features: {"sqrt", "log2", None}, int or float, default="sqrt" The number of features to consider when looking for the best split: - If int, then consider `max_features` features at each split. - If float, then `max_features` is a fraction and `max(1, int(max_features * n_features_in_))` features are considered at each split. - If "sqrt", then `max_features=sqrt(n_features)`. - If "log2", then `max_features=log2(n_features)`. - If None, then `max_features=n_features`. .. versionchanged:: 1.1 The default of `max_features` changed from `"auto"` to `"sqrt"`. Note: the search for a split does not stop until at least one valid partition of the node samples is found, even if it requires to effectively inspect more than ``max_features`` features.	'sqrt'
	max_leaf_nodes max_leaf_nodes: int, default=None Grow trees with ``max_leaf_nodes`` in best-first fashion. Best nodes are defined as relative reduction in impurity. If None then unlimited number of leaf nodes.	None
	min_impurity_decrease min_impurity_decrease: float, default=0.0 A node will be split if this split induces a decrease of the impurity greater than or equal to this value. The weighted impurity decrease equation is the following:: N_t / N * (impurity - N_t_R / N_t * right_impurity - N_t_L / N_t * left_impurity) where ``N`` is the total number of samples, ``N_t`` is the number of samples at the current node, ``N_t_L`` is the number of samples in the left child, and ``N_t_R`` is the number of samples in the right child. ``N``, ``N_t``, ``N_t_R`` and ``N_t_L`` all refer to the weighted sum, if ``sample_weight`` is passed. .. versionadded:: 0.19	0.0
	bootstrap bootstrap: bool, default=True Whether bootstrap samples are used when building trees. If False, the whole dataset is used to build each tree.	True
	oob_score oob_score: bool or callable, default=False Whether to use out-of-bag samples to estimate the generalization score. By default, :func:`~sklearn.metrics.accuracy_score` is used. Provide a callable with signature `metric(y_true, y_pred)` to use a custom metric. Only available if `bootstrap=True`. For an illustration of out-of-bag (OOB) error estimation, see the example :ref:`sphx_glr_auto_examples_ensemble_plot_ensemble_oob.py`.	False
	n_jobs n_jobs: int, default=None The number of jobs to run in parallel. :meth:`fit`, :meth:`predict`, :meth:`decision_path` and :meth:`apply` are all parallelized over the trees. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details.	None
	random_state random_state: int, RandomState instance or None, default=None Controls both the randomness of the bootstrapping of the samples used when building trees (if ``bootstrap=True``) and the sampling of the features to consider when looking for the best split at each node (if ``max_features < n_features``). See :term:`Glossary ` for details.	None
	verbose verbose: int, default=0 Controls the verbosity when fitting and predicting.	0
	warm_start warm_start: bool, default=False When set to ``True``, reuse the solution of the previous call to fit and add more estimators to the ensemble, otherwise, just fit a whole new forest. See :term:`Glossary ` and :ref:`tree_ensemble_warm_start` for details.	False
	class_weight class_weight: {"balanced", "balanced_subsample"}, dict or list of dicts, default=None Weights associated with classes in the form ``{class_label: weight}``. If not given, all classes are supposed to have weight one. For multi-output problems, a list of dicts can be provided in the same order as the columns of y. Note that for multioutput (including multilabel) weights should be defined for each class of every column in its own dict. For example, for four-class multilabel classification weights should be [{0: 1, 1: 1}, {0: 1, 1: 5}, {0: 1, 1: 1}, {0: 1, 1: 1}] instead of [{1:1}, {2:5}, {3:1}, {4:1}]. The "balanced" mode uses the values of y to automatically adjust weights inversely proportional to class frequencies in the input data as ``n_samples / (n_classes * np.bincount(y))`` The "balanced_subsample" mode is the same as "balanced" except that weights are computed based on the bootstrap sample for every tree grown. For multi-output, the weights of each column of y will be multiplied. Note that these weights will be multiplied with sample_weight (passed through the fit method) if sample_weight is specified.	None
	ccp_alpha ccp_alpha: non-negative float, default=0.0 Complexity parameter used for Minimal Cost-Complexity Pruning. The subtree with the largest cost complexity that is smaller than ``ccp_alpha`` will be chosen. By default, no pruning is performed. See :ref:`minimal_cost_complexity_pruning` for details. See :ref:`sphx_glr_auto_examples_tree_plot_cost_complexity_pruning.py` for an example of such pruning. .. versionadded:: 0.22	0.0
	max_samples max_samples: int or float, default=None If bootstrap is True, the number of samples to draw from X to train each base estimator. - If None (default), then draw `X.shape[0]` samples. - If int, then draw `max_samples` samples. - If float, then draw `max(round(n_samples * max_samples), 1)` samples. Thus, `max_samples` should be in the interval `(0.0, 1.0]`. .. versionadded:: 0.22	None
	monotonic_cst monotonic_cst: array-like of int of shape (n_features), default=None Indicates the monotonicity constraint to enforce on each feature. - 1: monotonic increase - 0: no constraint - -1: monotonic decrease If monotonic_cst is None, no constraints are applied. Monotonicity constraints are not supported for: - multiclass classifications (i.e. when `n_classes > 2`), - multioutput classifications (i.e. when `n_outputs_ > 1`), - classifications trained on data with missing values. The constraints hold over the probability of the positive class. Read more in the :ref:`User Guide `. .. versionadded:: 1.4	None