Stabilizing Estimates of Shapley Values with Control Variates

Abstract

Shapley values are among the most popular tools for explaining predictions of black-box machine learning models. However, their high computational cost motivates the use of sampling approximations, inducing a considerable degree of uncertainty. To stabilize these model explanations, we propose ControlSHAP, an approach based on the Monte Carlo technique of control variates. Our methodology is applicable to any machine learning model and requires virtually no extra computation or modeling effort. On several high-dimensional datasets, we find it can produce dramatic reductions in the variability of Shapley estimates.

Appendices

Appendix

A Exact Shapley Values

1.1 A.1 Shapley Value, Assuming Feature Independence

Proof

Assume features are sampled independently in the value function. Consider a second-order Taylor approximation to f at x:

$$\begin{aligned} g(x') := f(x) + (x'-x)^T \nabla f(x) + \frac{1}{2}(x'-x)^T \nabla ^2 f(x) (x'-x). \end{aligned}$$

Let $J = \nabla f(x)$, $H = \nabla ^2 f(x)$, $\varSigma =\text {Cov}(X)$, and $\sigma ^2_{k\ell }=\text {Cov}(X_k, X_\ell )$. The Shapley value function is

$$\begin{aligned} v_x(S) &= E[g(X')|X'_S = x_S]\\ &= f(x) + E[(x'-x)|X'_S=x_S]^T \nabla f(x) \\ &\qquad \qquad + \frac{1}{2}E\Big [(x'-x)^T \nabla ^2 f(x) (x'-x)|X'_S=x_S\Big ]\\ &= f(x) + (\mu _{S^C} - x_{S^C})^T J_{S^C}\\ &\qquad \qquad + \frac{1}{2}E\Big [(x_{S^C}'-x_{S^C})^T \nabla ^2 f(x)_{{S^C}{S^C}} (x_{S^C}'-x_{S^C})\Big ], \end{aligned}$$

where the quadratic term is equal to

$$\begin{aligned} &tr(H_{{S^C}{S^C}}\varSigma _{{S^C}{S^C}}) + (\mu _{S^C}- x_{S^C})^T H_{{S^C}{S^C}}(\mu _{S^C}- x_{S^C})\\ = &\sum _{k\in {S^C}}\sum _{\ell \in {S^C}} H_{k\ell }\Big (\sigma ^2_{k\ell } + (\mu _k - x_k)(\mu _\ell - x_\ell )\Big )\\ = &H_{jj}(\sigma ^2_{jj} + (\mu _j-x_j)^2). \end{aligned}$$

Recall $\displaystyle \phi _j(x) := \frac{1}{d} \sum _{S \subseteq [d]\backslash \{j\}} {d-1\atopwithdelims ()|S|}^{-1} \big (v_x(S \cup \{j\}) - v_x(S)\big )$. The difference in value functions is

Define the Shapley weight $w_S = \frac{1}{d}{d-1\atopwithdelims ()|S|}^{-1}$.

(10)

Noting subsets of equal size have the same Shapley weight, we can easily show $\displaystyle \sum _{S \subseteq [d]\backslash \{j\}}w_S = 1$.

$$\begin{aligned} \begin{aligned} \sum _{S \subseteq [d]\backslash \{j\}}w_S &= \sum _{S \subseteq [d]\backslash \{j\}} \frac{1}{d}{d-1\atopwithdelims ()|S|}^{-1}\\ &= \sum _{a=0}^{d-1} {d-1\atopwithdelims ()a} \frac{1}{d}{d-1\atopwithdelims ()a}^{-1} = \sum _{a=0}^{d-1}\frac{1}{d}=1 \end{aligned} \end{aligned}$$

(11)

With a bit more arithmetic, we can show $\displaystyle \sum _{S \subseteq [d]\backslash \{j, k\}}w_S = \frac{1}{2}$.

$$\begin{aligned} \begin{aligned} \sum _{S \subseteq [d]\backslash \{j, k\}}w_S &= \sum _{a=0}^{d-2} {d-2\atopwithdelims ()a} \frac{1}{d}{d-1\atopwithdelims ()a}^{-1}\\ &= \sum _{a=0}^{d-2} \frac{(d-2)!}{a!(d-a-2)!}\frac{a!(d-a-1)!}{d!}\\ &= \sum _{a=0}^{d-2} \frac{d-a-1}{d(d-1)}\\ &= \frac{1}{d(d-1)}\Big [d(d-1) - \sum _{a=0}^{d-2}(a+1) \Big ]\\ &= \frac{1}{d(d-1)}\Big [d(d-1) - \frac{d(d-1)}{2} \Big ] = \frac{1}{2}. \end{aligned} \end{aligned}$$

(12)

Plugging the results of 11 and 12 into 10 yields the final expression for the Shapley value:

   $\square$

1.2 A.2 Shapley Value, Correlated Features

Consider the linear model $g(x') = \beta ^T x' + b$, where $x' \sim \mathcal {N}(\mu , \varSigma )$. (For our problem, $\beta = \nabla f(x)$ and $b = f(x)$.)

Let $P_S$ be the projection matrix selecting set S; define $R_S = P_{\bar{S}}^T P_{\bar{S}} \varSigma P_S^T (P_S \varSigma P_S^T)^{-1} P_S$ and $Q_S = P_S^T P_S$. We consider the set of permutations of [d], each of which indexes a subset S as the features that appear before j. Averaging over all such permutations, the Shapley values of g(x) [1, 7, 23] are

$$\begin{aligned} \phi _j(x') &= \beta \underbrace{\left[ \frac{1}{d!} \sum _m ([Q_{\{S^m \cup j\}^C} - R_{S^m \cup j}]- [Q_{\{S^m\}^C} - R_{S^m}])\right] }_{C_j} \mu \\ &\qquad \qquad + \beta \underbrace{\left[ \frac{1}{d!} \sum _m ([Q_{S^m \cup j} + R_{S^m \cup j}] - [Q_{S^m} + R_{S^m}]) \right] }_{D_j} x'. \end{aligned}$$

Lastly, we observe that $C_j = -D_j$. For each subset S, the R terms cancel out; $Q_{S^C} + Q_S = I_d$, so $Q_{\{S\cup j\}^C} - Q_{S^C} = -(Q_{S\cup j} - Q_S)$. This yields our expression for the dependent Shapley value:

$$\begin{aligned} \phi _j^g(x) = \beta ^T D_j (x - \mu ). \end{aligned}$$

B Comprehensive Results

We display results across all 5 datasets and 3 machine learning predictors. All datasets were binary classification problems, so we used the same models to fit them. For logistic regression and random forest, we used the sklearn implementation with default hyperparameters. For the neural network, we fit a two-layer MLP in Pytorch with 50 neurons in the hidden layer and hyperbolic tangent activation functions.

Figure 6 displays the variance reductions of our ControlSHAP methods in the four settings: Independent vs Dependent Features, and Shapley Sampling vs KernelSHAP. The error bars span the 2^th to 7^th percentiles of the variance reductions for 40 held-out samples.

Figure 7 compares the average number of changes in rankings between the original and ControlSHAP Shapley estimates. Specifically, we look the Shapley estimates obtained via KernelSHAP, assuming correlated features.

Fig. 6.

Variance Reductions

Fig. 7.

Average Number of Rank Changes, with and without ControlSHAP’s adjustment via Control Variates (CV).

C Anticipated Correlation

Figure 8 compares the observed and anticipated variance reductions, across two combinations of dataset and predictor. Recall that the variance reduction of the control variate estimate for $\phi _j(x)$ is $\rho ^2(\hat{\phi }_j(x)^\text {model}, \hat{\phi }_j(x)^\text {approx})$, where $\rho$ is the Pearson’s correlation coefficient. We compute the sample correlation coefficient of the Shapley values on sample x as follows:

$$\begin{aligned} \hat{\rho } := \frac{\widehat{\text {Cov}}(\hat{\phi }_j(x)^\text {model}, \hat{\phi }_j(x)^\text {approx})}{\sqrt{\widehat{\text {Var}}(\hat{\phi }_j(x)^\text {model})\widehat{\text {Var}}(\hat{\phi }_j(x)^\text {approx})}} \end{aligned}$$

We average this across the 50 iterations to obtain a single estimate for the correlation $\hat{\rho }$. The plots display the median and error bars for $\hat{\rho }^2(\hat{\phi }_j(x)^\text {model}, \hat{\phi }_j(x)^\text {approx})$ across 40 samples. Once again, the error bars span the 2^th to 7^th percentiles.

Fig. 8.

Anticipated vs Observed Variance Reduction

D Comparing KernelSHAP Variance Estimators

Section 4.2 of the paper details methods for computing the variance and covariance between KernelSHAP estimates. Bootstrapping and the least squares covariance produce extremely similar estimates. In turn, these produce ControlSHAP estimates that reduce variance by roughly the same amount, as shown in Fig. 9. This indicates that both methods are appropriate choices.

In contrast, we were not able to get the grouped method to reliably work. Its variance estimates were somewhat erratic, as they are drawn from a heavy-tailed $\chi ^2$ distribution (Fig. 10). As a result, its ControlSHAP estimates were occasionally more variable than the original KernelSHAP values themselves.

Fig. 9.

Variance reductions of ControlSHAP across 40 samples with bootstrapped and least-squares estimates of KernelSHAP variance and covariance. Shapley estimates computed assuming independent features.

Fig. 10.

Correlation between model and Taylor approximation, grouped and least squares. Same input and feature across 50 iterations, with logistic regression on the bank dataset.