Explaining individual predictions when features are dependent: More accurate approximations to Shapley values

1. Introduction

• Interpretability matters because complex models are used in consequential applications, while limited transparency raises practical trust and legal concerns.Source: §1 · p. 1
• The paper focuses on model-agnostic local methods that learn simple, interpretable explanations for a specific prediction.Source: §1 · p. 1
• Shapley values allocate the difference between a prediction and the global average across features, unlike LIME’s local-average baseline; Lundberg & Lee (2017a) found stronger human agreement for Shapley explanations.Source: §1 · p. 2
• Kernel SHAP is computationally efficient but assumes feature independence, which can produce very wrong explanations even for linear models; Tree SHAP is tree-specific and may also be inaccurate.Tree SHAP accounts for some, but not all, dependence, according to the introduction.Source: §1 · p. 2
• The paper extends Kernel SHAP to handle dependent features and proposes clustering Shapley values for dependent features to improve individual-prediction explanations.The authors describe this as the first Shapley-based individual-prediction explanation method to account for feature dependence.Source: §1 · p. 2

2.1. The exact Shapley value and cooperative game theory.

• The Shapley value allocates a cooperative game’s total gains by averaging each player’s marginal contribution over coalitions, weighted across possible coalition formation orders.It is characterized as the unique distribution satisfying the listed desirable properties.Source: §2 · p. 3, §2.1 · pp. 3–4

2.2. Shapley values for prediction explanation.

• For prediction explanation, Shapley values decompose f(x*) into feature contributions relative to the global average, using expected model output conditional on observed features to define coalition worth.Without an explicit formula for general dependent-feature, nonlinear models, exact computation involves 2^M subsets; the paper identifies appropriate expectation estimation as its main methodological contribution.Source: §2.2 · p. 4, §2.1 · p. 4

2.3. Kernel SHAP.

• Kernel SHAP (Lundberg & Lee, 2017a) combines a computationally tractable Shapley-value approximation with a method for estimating coalition worth v(S).The paper restates the method in two parts: first computing approximate Shapley values when v(S) is known, then estimating v(S).Source: §2.3 · p. 5

2.3.1. Approximated weighted least squares.

• Kernel SHAP expresses Shapley values as the solution to a weighted least-squares problem and approximates it by sampling subsets according to the Shapley kernel weights.The empty and full feature subsets are excluded from sampling and appended separately with large weights.Source: §2.3.1 · p. 5

2.3.2. Estimating the contribution function under feature independence.

• Kernel SHAP estimates coalition worth under feature independence by replacing the conditional distribution of missing features with their marginal distribution and sampling from training data.The sampled training features are drawn independently of the observed feature values.Source: §2.3.2 · p. 6

3. Incorporating dependence into the Kernel SHAP method

• The authors propose Gaussian, Gaussian-copula, empirical conditional, and hybrid approaches to estimate p(x̄S|xS = x*), replacing independent sampling in Kernel SHAP.The independence assumption enters when approximating the integral; under strong dependence, it can produce a completely wrong answer.Source: §3 · p. 6

3.1. Multivariate Gaussian distribution.

• Under a multivariate Gaussian assumption, the conditional distribution is Gaussian, with full mean and covariance estimated from training data for conditional sampling.Those conditional samples are then used to approximate the integral in Kernel SHAP.Source: §3.1 · pp. 6–7

3.2. Gaussian copula.

• When marginals are non-Gaussian, the Gaussian copula approach models them empirically and represents dependence with a Gaussian latent distribution, enabling conditional sampling.It transforms empirical marginals to Gaussian features, samples conditionally, then maps samples back to the original marginals.Source: §3.2 · p. 7

3.3. Empirical conditional distribution.

• The empirical conditional approach weights training observations by their distance from the explained instance on features in S, then uses the weights to approximate the integral.The method uses a scaled Mahalanobis distance and Gaussian-kernel-like weights; incorporating weights directly in the integral uses each training observation only once.Source: §3.3 · pp. 7–8
• Bandwidth σ trades bias against variance; AICc can select it, but selection for every subset and instance is computationally intensive.The approximate AICc method shares σ across subsets of equal size, while σ → ∞ recovers the original Kernel SHAP method.Source: §3.3 · pp. 8–9
• The empirical approach works well when dim(xS) ≤ D*, but Gaussian or copula methods outperform it at higher conditioning dimensions, motivating a hybrid strategy.The hybrid uses empirical conditionals for subsets with dim(xS) ≤ D* and a parametric method for the rest; the threshold may be chosen by cross-validation.Source: §3.4 · p. 9

4. Experiments

• Because exact Shapley values require O(2^M) integrals, experiments use low and moderate dimensions; in three dimensions, exact values are tractable and the exact Kernel SHAP variant is used.The comparison includes Gaussian, copula, empirical, and combined approximations, with empirical-AICc-exact limited to three dimensions and combined approaches to ten dimensions.Source: §4 · pp. 9–10

4.1. Evaluation measures.

• The experiments compare MAE-based skill scores against original Kernel SHAP, with all methods limited to 1,000 training samples per feature combination and test observation.The skill score is 1 for a perfect approximation, 0 for original Kernel SHAP, and negative when a method performs worse than that reference.Source: §4.1 · p. 10

4.2.4. Linear model; 4.2.5. Piecewise constant model.

• The three-feature experiments combine Gaussian, Generalized Hyperbolic, and Gaussian-mixture features with linear or piecewise-constant sampling models.The piecewise-constant model is fitted with XGBoost using default hyperparameters, histogram tree learning, and 50 boosting iterations.Source: §4.2 · p. 10, §4.2.4 · p. 11, §4.2.5 · p. 12

4.2.1. Gaussian distributed features.

• The three-feature study varies Gaussian-feature correlation from 0 to 0.98 to represent increasing positive dependence.The features follow a multivariate Gaussian distribution with covariance matrix Σ(ρ).Source: §4.2.1 · p. 11

4.2.2. Skewed and heavy-tailed distributed features; 4.2.3. Multi-modal distributed features.

• The study also tests skewed, heavy-tailed Generalized Hyperbolic features and a two-component Gaussian mixture with varying distance between modes.The GH skewness coefficient κ ranges from 1 to 10; the mixture distance parameter γ ranges from 0.5 to 10.Source: §4.2.2 · p. 11, §4.2.3 · p. 11

4.2.6. Results.

• For Gaussian features, the Gaussian method generally performed best; for bimodal features empirical methods led, while Gaussian or copula methods were preferable at larger skewness.In the linear Gaussian-feature experiment, all methods beat original Kernel SHAP once correlation exceeded 0.05, although original Kernel SHAP worked well for independent features.Source: §4.2.6 · p. 12
• TreeSHAP was only slightly better than original Kernel SHAP in piecewise-constant experiments D and E, and worse in F.This was surprising because TreeSHAP is supposed to handle dependence better than original Kernel SHAP.Source: §4.2.6 · pp. 12–13
• Across the three-dimensional experiments, dependence-aware methods outperformed original Kernel SHAP, the best method varied by feature distribution, and approximate-AICc empirical methods were recommended.The approximate-AICc results were fairly similar to, and sometimes better than, those from the slower exact-AICc method.Source: §4.2.6 · p. 13

Dimension 10

• The linear and piecewise-constant 10D experiments extended earlier 3D setups, excluding the effect of the 10th feature on the outcome.The linear model used ordinary linear regression, while the piecewise-constant model used XGBoost with the same settings as the 3D case.Source: §4.3 · p. 16
• The final simulation used a 10D Generalized Hyperbolic feature distribution chosen to resemble the real-data example, with the piecewise-constant model as its sampling model.Source: §4.3 · p. 16
• The 10D experiments omitted empirical-AICc-exact because it was computationally intensive and similar to empirical-AICc-approx in 3D; the final experiment used empirical-0.1.Combined approaches used empirical estimates for subsets of dimension ≤3, with bandwidth σ set by approximate AICc in the linear case or fixed at 0.1.Source: §4.3 · p. 16
• For Gaussian features, results largely matched the 3D counterparts: Gaussian generally performed best, combined empirical/Gaussian or copula approaches worked well, and TreeSHAP beat Kernel SHAP only at small or medium dependence.For piecewise-constant models, TreeSHAP was worse than Kernel SHAP under high dependence.Source: §4.3 · pp. 16–17, Figure 8 · p. 17, Figure 9 · p. 17
• For skewed, heavy-tailed features, empirical-0.1+Gaussian had slightly better MAE values and skill scores than empirical-0.1+copula, and all suggested approaches outperformed Kernel SHAP and TreeSHAP.The 10D exact Shapley values were computed using Monte Carlo integration rather than numerical integration.Source: §4.3 · pp. 16–17

4.4. Real data example.

• The transaction data yielded 28 standardized features with skewed, heavy-tailed distributions and groups of highly correlated features.Features summarized daily balances, transaction counts, and transfers over 365 days; most pairwise rank correlations were near zero, but some groups had high mutual correlations.Source: §4.4 · p. 18
• An XGBoost model trained on 12,696 observations achieved test AUC 0.88 on 1,921 observations.The model used 50 trees with default parameter settings.Source: §4.4 · p. 18
• For inn std mean and sk std mean, the two methods’ Shapley estimates corresponded closely, but they differed greatly for tr std and br max, which were more dependent.The authors suggest that the original Kernel SHAP approach performs well for the nearly independent features and poorly for the more dependent ones.Source: §4.4 · p. 18

5. Clustering of Shapley values corresponding to feature dependence

• Grouping dependent features can surface their combined Shapley value among the top ten even when none ranks individually, helping present explanations more correctly.The authors suggest ranking absolute Shapley values, but note that grouping or clustering may ease interpretation when features are dependent.Source: §5 · p. 18
• The authors recommend clustering once on training data or sufficiently large test data, since individual-specific clusters could make explanations across people opaque.They describe their preferred dependence measures and clustering methods as noncanonical, allowing other choices to suit an application.Source: §5 · p. 19

5.1. Dependence measure.

• They use Kendall’s τ rank correlation and define dissimilarity as D_j,k = 1 − |τ_j,k|, grouping positive or negative associations; perfect correlation gives 0 and independence 1.The absolute value makes the dissimilarity depend on association strength regardless of its sign.Source: §5.1 · p. 19, §5.2 · p. 19

5.2. Clustering method and cluster size.

• They apply hierarchical agglomerative clustering with complete linkage and choose cluster count by a user-set dissimilarity level or the minimum Kelley-Gardner-Sutcliffe penalty.In their example, α = 0.1 yielded twelve groups.Source: §5.2 · pp. 19–20, §5.3 · p. 21
• For two example predictions, empirical-0.1+Gaussian produced group rankings regarded as more correct than original Kernel SHAP, whose group-12 attribution even reversed direction in example 2.The authors regard the dependence-aware method as more trustworthy based on their simulation experiments.Source: §5.3 · p. 21

6. Summary and discussion

• The modified Kernel SHAP estimates conditional feature distributions using Gaussian, Gaussian-copula, empirical, or combined approaches.The nonparametric approach did better when conditioning on few features, while Gaussian and copula methods did better with more conditioned features; the authors regard the combined approach as most promising.Source: §6 · pp. 22–23
• In simulations, dependent-feature methods more accurately approximated true Shapley values than original Kernel SHAP and outperformed Tree SHAP on nonlinear models.The combined approach was applied to a 28-variable mortgage-default case; although true values were unknown, the authors report indications that it gives more sensible approximations than Tree SHAP and original Kernel SHAP.Source: §6 · p. 23
• Hierarchical clustering groups dependent features, whose individual Kernel SHAP values are summed and presented as group contributions.Figure 13 illustrates contributions from twelve feature groups for two mortgage-default predictions using waterfall charts and group-value plots.Source: Figure 13 · p. 22, §6 · p. 23
• The method’s main computational cost is nonparametric bandwidth selection; with fixed bandwidth, the combined approach matches original Kernel SHAP’s computational time.Graph-based community decomposition has been explored by Chen et al. (2018) and Li et al. (2018), but tabular data may require many conditional-independence tests.Source: §6 · p. 23
• For probability predictions, decomposing probabilities can yield Shapley-value sums outside [0,1], while log-odds contributions may be harder for people to interpret.The paper treats the choice of what to decompose as a practical rather than mathematical question.Source: §6 · p. 23
• Categorical features complicate the proposed approaches: nonparametric methods can use numerical encodings, while parametric methods remain more challenging.Options discussed include one-hot encoding, categorical or mixed-type distances (Huang, 1997, 1998; de Leon & Carriere (2005)), and entity embeddings (Guo & Berkhahn, 2016).Source: §6 · pp. 23–24

Appendix A. Shapley properties in the prediction explanation setting

• The baseline φ0 = E[f(x)] measures the prediction portion not due to features; features are “not important” for a prediction when |φ0| is large relative to their φj values.The baseline is the global average prediction.Source: Appendix A · p. 24
• Efficiency requires feature Shapley values to sum to f(x*) − E[f(x)], allocating the prediction’s difference from the global mean to features and enabling comparison across predictions.Formally, ∑_j φj = f(x*) − E[f(x)].Source: Appendix A · p. 24
• Symmetry requires two features with equal contributions across every combination with other features to receive equal Shapley values, avoiding inconsistent explanations.Source: Appendix A · p. 24
• The dummy-player property assigns zero Shapley value to any feature that never changes the prediction, preventing explanations for features without influence.Source: Appendix A · p. 24
• Linearity makes a feature’s value for a sum or linear combination of prediction functions equal to the sum of its values for each function, supporting interpretation of ensembles.The passage names Random Forests and other structurally simple ensemble models as examples.Source: Appendix A · p. 24
• Shapley values are the only additive explanation method satisfying all these properties (Lundberg & Lee, 2017a).Source: Appendix A · p. 24

• The appendix proves the explicit Shapley-value formula for linear regression with independent features, then shows how to obtain v(S) when features may be dependent.Source: Appendix B · p. 24

B.1. Linear model and independent features.

• With a linear regression model and independent features, the Shapley contribution of a feature is the same for every coalition of other features.This makes the Shapley formula simplify by summing identical marginal contributions with weights that sum to 1.Source: §B.1 · pp. 24–25
• The derivation computes v(S) and v(S ∪ {j}); linearity and feature independence are used in successive steps to obtain the expression for v(S).The passage attributes the transition from step (18) to (19) to linearity and from (19) to (20) to independence.Source: §B.1 · p. 25
• The proof’s simplified Shapley expression follows after the coalition-independent difference is combined with Shapley weights summing to 1.Source: §B.1 · p. 25

B.2. Linear model and dependent features.

• The needed quantity is a proper estimate of E[x_{S}|x_S=x*_S].Source: §B.2 · p. 25
• An analytical estimate of this conditional expectation can avoid time-consuming simulations.Source: §B.2 · p. 25
• Obtaining a proper analytical estimate is not straightforward in general.Source: §B.2 · p. 25

Appendix C. Generalized Hyperbolic Distribution

• The appendix introduces the density of a d-dimensional Generalized Hyperbolic random vector X.Source: Appendix C · p. 26
• The mean is µ + E(W)β and the covariance is E(W)Σ + Var(W)ββT; δ is the squared Mahalanobis distance and Kλ is a modified Bessel function.Source: Appendix C · p. 26
• For a partitioned vector, the conditional distribution of X2 given X1 = x1 is GH distributed, using a slightly different parameterization for the conditional case.The conditional parameterization follows McNeil et al. (2006); the conditional-distribution result is attributed to Wei et al. (2019).Source: Appendix C · p. 26