Model independent feature attributions: Shapley values that uncover non-linear dependencies

Recognising dependence: Example 1

For example, the empirical R² characteristic function ${\hat{R}}_{\mathbf{y}}$ is given by (4)${\displaystyle {\hat{R}}_{\mathbf{y}}\left(S\right)=1-\frac{\left|\rho \left(\mathbf{y},\mathbf{X}{|}_S\right)\right|}{\left|\rho \left(\mathbf{X}{|}_S\right)\right|},}$ where ρ is the empirical Pearson correlation matrix.

Regardless of whether we use a population measure or an estimate, the R² measures only the linear relationship between the response (i.e., labels) Y and features X. This implies the R² may perform poorly as a measure of dependence in the presence of non-linearity. The following example from a non-linear DGP demonstrates this point.

Suppose the features X_j, j ∈ [d] are independently uniformly distributed on [ − 1, 1]. Given a diagonal matrix A = diag(a₁, …, a_d), let the response variable Y be determined by the quadratic form (5)${\displaystyle Y=X^{T}AX=a_1{X}_1^2+\dots +a_d{X}_d^2.}$ Then, the covariance Cov(Y, X_j) = 0 for all j ∈ [d]. This is because ${\displaystyle Cov\left(X^{T}AX,X_j\right)=\sum _{j=1}^{d}Cov\left(X_j^2,X_j\right)=0,}$ since E[X_j] = 0 and $\text{E}\left[X_j^3\right]=0$. In Fig. 1, we display the X₄ cross section of 10, 000 observations generated from Eq. (5) with d = 5 and A = diag(0, 2, 4, 6, 8), along with the least squares line of best fit and associated R² value. We visualize the results for the corresponding Shapley decomposition in Fig. 2. As expected, we see that the R² is not able to capture the non-linear dependence structure of Eq. (5), and thus neither is its Shapley decomposition.

We note that improvements on the results in Figs. 1 and 2 can be obtained by choosing a suitable linearising transformation of the features or response prior to calculating R², but such a transformation is not known to be discernible from data in general, except in the simplest cases.

Distance correlation and affine invariant distance correlation

The distance correlation, and its affine invariant adaptation, were both introduced by Székely, Rizzo & Bakirov (2007). Unlike the Pearson correlation, the distance correlation between Y and X is zero if and only if Y and X are statistically independent. However, the distance correlation is equal to 1 only if the dimensions of the linear spaces spanned by Y and X are equal, almost surely, and Y is a linear function of X.

First, the population distance covariance between the response Y and feature vector X is defined as a weighted L₂ norm of the difference between the joint characteristic function¹ , f_YX and the product of marginal characteristic functions f_Yf_X. In essence, this is a measure of squared deviation from the assumption of independence, i.e., the hypothesis that f_YX = f_Yf_X.

The empirical distance covariance $V_n^2$ is based on Euclidean distances between sample elements, and can be computed from data matrices Y, X as (6)${\displaystyle {\hat{V}}^2\left(\mathbf{Y},\mathbf{X}\right)=\sum _{i,j=1}^{n}A{\left(\mathbf{Y}\right)}_{ij}A{\left(\mathbf{X}\right)}_{ij},}$ where the matrix function A(W) for W ∈ {Y, X} is given by ${\displaystyle A{\left(\mathbf{W}\right)}_{ij}=B{\left(\mathbf{W}\right)}_{ij}-\frac{1}{n}\sum _{i=1}^{n}B{\left(\mathbf{W}\right)}_{ij}-\frac{1}{n}\sum _{j=1}^{n}B{\left(\mathbf{W}\right)}_{ij}+\frac{1}{n^2}\sum _{i,j=1}^{n}B{\left(\mathbf{W}\right)}_{ij},}$ where ||⋅|| denotes the Euclidean norm, and B(W) is the n × n distance matrix with B(W)_ij = ||w_i − w_j||,  where w_i denotes the ith observation (row) of W. Here, Y is in general a matrix of observations, with potentially multiple features. Notice the difference between Y and y, where the latter is the (single column) label vector introduced in “Attributed dependence on labels”.

The empirical distance correlation $\widehat{ℛ}$ is given by (7)${\displaystyle {\widehat{ℛ}}^2\left(\mathbf{Y},\mathbf{X}\right)=\frac{{\hat{V}}^2\left(\mathbf{Y},\mathbf{X}\right)}{\sqrt{{\hat{V}}^2\left(\mathbf{Y},\mathbf{Y}\right){\hat{V}}^2\left(\mathbf{X},\mathbf{X}\right)}},}$ for ${\hat{V}}^2\left(\mathbf{Y},\mathbf{Y}\right){\hat{V}}^2\left(\mathbf{X},\mathbf{X}\right)\ne 0$, and $\widehat{ℛ}\left(\mathbf{Y},\mathbf{X}\right)=0$ otherwise. For our purposes, we define the distance correlation characteristic function estimator (8)${\displaystyle {\hat{D}}_{\mathbf{y}}\left(S\right)={\widehat{ℛ}}^2\left(\mathbf{y},\mathbf{X}{|}_S\right).}$ A transformation of the form x↦Ax + b for a matrix A and vector b is called affine. Affine invariance of the distance correlation is desirable, particularly in the context of hypothesis testing, since statistical independence is preserved under the group of affine transformations. When Y and X are first scaled as ${\mathbf{Y}}^{\prime }=\mathbf{Y}{S}_{\mathbf{Y}}^{-1∕2}$ and ${\mathbf{X}}^{\prime }=\mathbf{X}{S}_{\mathbf{X}}^{-1∕2}$, the distance correlation $\hat{V}\left({\mathbf{Y}}^{\prime },{\mathbf{X}}^{\prime }\right)$, becomes invariant under any affine transformation of Y and X (Székely, Rizzo & Bakirov, 2007, Section 3.2). Thus, the empirical affine invariant distance correlation is defined by (9)${\displaystyle {\widehat{ℛ}}^{\prime }\left(\mathbf{Y},\mathbf{X}\right)=\widehat{ℛ}\left(\mathbf{Y}{S}_{\mathbf{Y}}^{-1∕2},\mathbf{X}{S}_{\mathbf{X}}^{-1∕2}\right),}$ and we define the associated characteristic function estimator ${\hat{D}}_{\mathbf{y}}^{\prime }$ in the same manner as Eq. (8). Monte Carlo studies regarding the properties of these measures are given by Székely, Rizzo & Bakirov (2007).

Hilbert–Schmidt independence criterion

The Hilbert Schmidt Independence Criterion (HSIC) is a kernel-based independence criterion, first introduced by Gretton et al. (2005a). Kernel-based independence detection methods have been adopted in a wide range of areas, such as independent component analysis (Gretton et al., 2007). The link between energy distance-based measures, such as the distance correlation, and kernel-based measures, such as the HSIC, was established by Sejdinovic et al. (2013). There, it is shown that the HSIC is a certain formal extension of the distance correlation.

The HSIC makes use of the cross-covariance operator, C_YX, between random vectors Y and X, which generalises the notion of a covariance. The response Y and feature vector X are each mapped to functions in a Reproducing Kernel Hilbert Spaces (RKHS), and the HSIC is defined as the Hilbert–Schmidt (HS) norm $\left|\right|C_{YX}|{|}_{HS}^2$ of the cross-covariance operator between these two spaces (Gretton et al., 2005b; Gretton et al., 2007; Gretton et al., 2005a). Given two kernels ℓ, k, associated to the RKHS of Y and X, respectively, and their empirical evaluation matrices L, K with row i and column j elements ℓ_ij = ℓ(y_i, y_j) and k_ij = k(x_i, x_j), where y_i, x_i denote the ith observation (row) in data matrices X and Y, respectively, the empirical HSIC can be calculated as (10)${\displaystyle \widehat{\mathrm{HSIC}}\left(\mathbf{Y},\mathbf{X}\right)=\frac{1}{n^2}\sum _{i,j}^n{k}_{ij}{\ell }_{ij}+\frac{1}{n^4}\sum _{i,j,q,r}^n{k}_{ij}{\ell }_{qr}-\frac{2}{n^3}\sum _{i,j,q}^n{k}_{ij}{\ell }_{iq}.}$ As in “Distance correlation and affine invariant distance correlation”, notice the difference between Y and y, where the latter is the (single column) label vector introduced in “Attributed dependence on labels”. Intuitively, this approach endows the cross-covariance operator with the ability to detect non-linear dependence, and the HS norm measures the combined magnitude of the resulting dependence. For a thorough discussion of positive definite Kernels, with a machine learning emphasis, see the work of Hein & Bousquet (2004).

Calculating the HSIC requires selecting a kernel. The Gaussian kernel is a popular choice that has been subjected to extensive testing in comparison to other kernel methods [see, e.g., Gretton et al., 2005a]. For our purposes, we define the empirical HSIC characteristic function by (11)${\displaystyle {\hat{H}}_y\left(S\right)=\widehat{\mathrm{HSIC}}\left(\mathbf{y},\mathbf{X}{|}_S\right),}$ and use a Gaussian kernel. Figure 2 shows the Shapley decomposition of $\hat{H}$ amongst the features generated from Eq. (5), again with d = 5 and A = diag(0, 2, 4, 6, 8). The decomposition has been normalised for comparability with the other measures of dependence presented in the figure. The HSIC can also be generalised to provide a measure of mutual dependence between any finite number of random vectors (Pfister et al., 2016).

Demonstration with concept drift

We illustrate the ADR and ADL techniques together with a simple and intuitive synthetic demonstration involving concept drift, where the DGP changes over time, impacting the mean squared prediction error (MSE) of a deployed XGBoost model. The model is originally trained with the assumption that the DGP is static, and the performance of the model is monitored over time with the intention of detecting violations of this assumption, as well as attributing any such violation to one or more features. A subset of the deployed features can be selected for scrutiny, by considering removal only of those selected features from the model. To highlight this, our simulated DGP has 50 features, and we perform diagnostics on 4 out of those 50 features.

For comparison, we compute SAGE values of the model mean squared error (Covert, Lundberg & Lee, 2020) and we compute the mean SHAP values of the logarithm of the model loss function (Lundberg et al., 2020). We will refer to the latter as SHAPloss. For SAGE and SHAPloss values, we employ a DGP similar to Eq. (14), with sample size 1000, but with X_i ≡ 0 for all i > 4. These features were nullified for tractability of the SAGE computation, since, unlike for Sunnies, the authors are not aware of any established method for selectively computing SAGE values of a subset of the full feature set. SAGE and SHAPloss were chosen for their popularity and ability to provide global feature importance scores.

At the initial time t = 0, we define the DGP as a function of temporal increments t ∈ ℕ∪{0}, (14)${\displaystyle Y=X_1+X_2+\left(1+\frac{t}{10}\right)X_3+\left(1-\frac{t}{10}\right)X_4+\sum _{i=5}^{50}{X}_i,}$ where X_i ∼ N(0, 4), for i = 1, 2, 3, 4, and X_i ∼ N(0, 0.05), for 5 ≤ i ≤ 50. Features 1 through 4 are the most effectual to begin with, and we can imagine that these were flagged as important during model development, justifying the additional diagnostic attention they enjoy after deployment. We see from Eq. (14) that, after deployment, i.e., during periods 1 ≤ t ≤ 10, the effect of X₄ decreases linearly to 0, while the effect of X₃ increases proportionately over time. In what follows, these changes are clearly captured by the residual and response dependence attributions of those features, using the DC characteristic function.

The results, with a sample size of n = 1000, from the DGP in Eq. (14), are presented in Fig. 3. According to the ADL (top), X₄ shows early signs of significantly reduced importance ${\varphi }_4\left({\hat{C}}_{\mathbf{y}}\right)$, as X₃ shows an increase in importance ${\varphi }_3\left({\hat{C}}_{\mathbf{y}}\right)$, which is roughly symmetrical to the decrease in ${\varphi }_4\left({\hat{C}}_{\mathbf{y}}\right)$. The ADR (bottom) show early significant signs that X₃ is disproportionately affecting the residuals, with high ${\varphi }_3\left({\hat{C}}_{\mathbf{e}}\right)$. The increase in residual attribution ${\varphi }_4\left({\hat{C}}_{\mathbf{e}}\right)$ is also evident, though the observation ${\varphi }_4\left({\hat{C}}_{\mathbf{e}}\right)<{\varphi }_3\left({\hat{C}}_{\mathbf{e}}\right)$ suggests that the drift impact from X₃ is the larger of the two.

The resulting SAGE and mean SHAPloss values are presented in Fig. 4. Interestingly, the behaviours of SHAPloss and SAGE are (up to scale and translation) analogous to the behaviour of ADL, rather than ADR, despite the model-independence of ADL. A reason for this, in this example, is that the feature with higher (resp. lower) dependence on Y contributes less (resp. more) to the residuals. When interpreting the SHAPloss and SAGE outputs, it is important to note that the model loss is increasing with t, since the true planar trend in (Y, X₃, X₄) rotates away from the learned trend. So, while the SAGE and SHAPloss results may appear to make the paradoxical suggestion that X₃ is utilised better by the model at t = 10 compared with t = 0, this is not the case: SAGE and SHAPloss are not accounting for the change in model loss over time. The model loss may decrease more when marginalising a feature under a misspecified model, than under a model with lower overall loss.

Demonstration with misspecified model

To illustrate the ADL, ADP and ADR techniques, we demonstrate a case where the model is misspecified on the training set, due to model bias. The inadequacy of this misspecified model is then detected on the validation set. Unlike the example given in “Model attributed dependence”, the DGP is unchanging between the two data sets. The key technique used in this demonstration is the comparison of differences between ADL (calculated in the absence of any model) and ADP (calculated using the output of a fitted model), in order to identify any differences in the attributions between dependence on labels and the dependence on the predictions produced by the misspecified model. Such a comparison, between model absence and model outputs, is not possible using purely model-dependent Shapley values.

To make this example intuitive, we avoid using a complex model such as XGBoost, in favour of a linear regression model. Since the simulated DGP is also linear, this example allows a simple comparison between the correct model and the misspecified model. The DGP in this example is (15)${\displaystyle Y=X_1+X_2+5X_3{X}_4{X}_5+\epsilon ,}$ where X₁, X₂, X₃ ∼ N(0, 1) are continuous, ε ∼ N(0, 0.1) is a small random error, and X₄, X₅ ∼ Bernoulli(1∕2) are binary. Hence, we can make the interpretation that the effect of X₃ is modulated by X₄ and X₅, such that X₃ is effective, only if X₄ = X₅ = 1. For this demonstration, we fit a misspecified linear model EY = β₀ + βX, where $X^T={\left(X_i\right)}_{i=1}^d$ is the vector of features, and ${\beta }_0,\beta ={\left({\beta }_i\right)}_{i=1}^d$ are real coefficients. This is a simple case where the true DGP is unknown to the analyst, who therefore seeks to summarise the 80 marginal contributions from 5 features into 5 Shapley values.

Figure 5 shows the outputs for attributed dependence on labels, residuals and predictions, via ordinary least squares estimation. From these results, we make the following observations:

1. For X₃ the ADP is significantly higher than the ADL.

2. For X₄ and X₅ the ADP is significantly lower than the ADL.

3. For X₁, X₂ there is no significant difference between ADP and ADL.

4. For X₁, X₂ ADR is negative, while X₃, X₄, X₅ have positive ADR.

Observations (i) and (ii) suggest that the model EY = β₀ + βX overestimates the importance of X₃ and underestimates the importance of X₄ and X₅. Observations (iii) and (iv) suggest that the model may adequately represent X₁, X₂, X₃, but that X₃, X₄ and X₅ are significantly more important for determining structure in the residuals than X₁ and X₂. A residuals versus fits plot may be useful for confirming that this structure is present and of large enough magnitude to be considered relevant.

Having observed the result in Fig. 5, for the misspecified linear model EY = β₀ + βX, we now fit the correct model: EY = β₀ + βX + X₃X₄X₅, which includes the three-way interaction effect X₃X₄X₅. The results, shown in Fig. 6, show no significant difference between the ADL and ADP for any of the features, and no significant difference in ADR between the features.