Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
84,092 characters · 24 sections · 129 citation commands
Variable importance without impossible data
Black-box models commonly attain state of the art prediction accuracy as measured by predictions of held-out response values. In many settings, it is not enough to know that they are accurate. We want to interpret them too. We might derive a useful scientific insight from understanding how a prediction works, or we might spot a flaw in the model that convinces us that it will not generalize well beyond our data set. We might also find that the algorithm works in a way that is unfair to some people.
A first step in interpreting a black-box model is to discern what are the important variables in that model. This is already a hard problem, and as we note below there are several large branches of the literature devoted to this problem. Most of the methods in those areas work by making changes to some, but not necessarily all, of the input variables in a model and taking note of how the predictions change in response. A severe problem with those approaches is that they will use unlikely or even impossible combinations of inputs. For instance, if one combines the youngest age of a criminal defendant with the greatest observed number of prior arrests, the result is an extremely unlikely synthetic data point. Combining features from different observations can lead to synthetic points describing people who graduated from high school before they were born, or who live in both Idaho and Los Angeles County. In medical settings, one can get impossible patients with systolic blood pressures below their diastolic blood pressure, or whose minimum O$_2$ saturation level is above their average.
Figure (ref) shows what this can look like when levels of one variable are permuted with respect to levels of another variable. The left panel uses the well known Boston housing data of harr:rubi:1978. The right panel uses the bodyfat data from the Data Analysis and Story Library (\url{https://dasl.datadescription.com}) after converting the waist measurements there from inches to centimeters. We see that the synthetic variable combinations can be quite unlike the real data. Those combinations may simply have nothing to do with any parts of Boston or with any likely human physique. This poses a problem: a black-box prediction function evaluated at such points is then being used far from where it was trained. Such extrapolations are at a minimum unreliable, and potentially meaningless. Furthermore, there are adversarial attacks on variable importance measures that use such inputs. slac:etal:2020 report how an algorithm that makes unfair predictions at observed data values can be rated as fair due to predictions made at synthetic data points formed by mixing and matching inputs. This applies to two of the leading methods in machine learning: the Local Interpretable Model agnostic Interpretation (LIME) of ribe:etal:2016 and the SHapley Additive exPlanations (SHAP) of lund:lee:2017. In contrast, the TreeSHAP method lund:2020 uses only observed values; we describe below how it differs from our proposal. For now, we note that it is only defined for tree structured predictions.
For the reasons above, we need a variable importance measure that does not use impossible data. While there are algorithmic methods to judge when a data point is `out of distribution,' there is not a reliable way to say exactly where the support of a distribution ends. In this article, we present the cohort Shapley measure from mase:owen:seil:2019 as one solution. Like many of the methods we discuss, this algorithm is grounded in axioms of economic game theory from shap:1953. However, the way the underlying game value is defined allows us to only use observed combinations of predictor variables.
This article is organized as follows: Section (ref) compares properties of cohort Shapley to some widely adopted variable importance methods and introduces a running example for illustration. Section (ref) provides a brief discussion of the problems of variable importance and points to literature surveys in three areas: statistics, machine learning and global sensitivity analysis. At first sight, variable importance looks to be easy because all the counterfactuals we might want to consider can usually be directly computed. At a second look, variable importance is actually extremely difficult because it seeks to identify causes of effects instead of effects of causes. We explain that distinction making reference to dawi:musi:2021. Section (ref) presents functional ANOVA, Sobol' indices from global sensitivity analysis, and Shapley value from game theory. Section (ref) describes how Shapley value has been used to quantify variable importance and presents the cohort Shapley method. Section (ref) illustrates cohort Shapley on an algorithmic fairness problem using data from angw:2016 on the COMPAS algorithm for predicting recidivism. The key issue there is whether Black defendants were treated unfairly. For every defendant, cohort Shapley can estimate the importance to their prediction of each of their feature variables. We are able to do this even though the algorithm was trained without using race as a variable. Furthermore, we are able to do this just using the algorithm's predictions, which is an advantage because the algorithm is proprietary and not available to us. We can aggregate importance over subsets of observations, consistent with an additivity axiom. We also show how to use the Bayesian bootstrap of rubi:1981 to get confidence intervals for variable importance of individuals and of aggregates. Section (ref) has some final remarks. Python code for the cohort Shapley method is available at \url{https://github.com/cohortshapley/cohortshapley}.
Here we compare cohort Shapley to some other methods that we will describe in later sections. Variable importance methods differ primarily in how they define importance. We base our comparisons on which useful properties they have or do not have, as shown in Table (ref). We have phrased the properties in a way that makes positive values desirable. Attributing importance to variables that are not in the model is controversial. We think that the ability to do that is desirable; cohort Shapley allows that but does not require it. TreeSHAP is fully automatic. Baseline Shapley and most other methods require defining a baseline point to compare the target point to. Cohort Shapley also requires a user-provided notion of variable similarity. We have listed LIME as requiring a second black box because it has a regularization process in determining how local an observation must be to be similar to the target. Kernel SHAP is listed as not computing the Shapley value because it computes a kernelized approximation to it. Similarly, on-manifold SHAP and conditional kernel SHAP do not compute the Shapley value; cohort and baseline Shapley do compute it. We include Monte Carlo versions of them because they are consistent for the Shapley value as computation increases. LIME requires the choice a surrogate model and a kernel, so we do not consider it to be automatic.
Of the methods in Table (ref), only cohort Shapley can be used on the COMPAS data that consist only of the predictions from a proprietary algorithm with no information on the model itself. It is similarly the only method that works with residuals to compare actual versus predicted levels of recidivism. The residuals are important because they let us study how an algorithm unfairly over-estimates or under-estimates a risk. Also if a variable is important for the residual it suggests that the model did not use that variable well. Cohort Shapley is one of only two that can attribute importance to an unused variable and one of only two that never uses unobserved data.
All the methods that we compare in Table (ref) can be used without retraining the model on different feature sets. They are aimed at attributing a specific point's prediction to one of their features, which is the task we consider in our numerical examples. Retraining with and without certain features present has been used, for example by strumbelwithretrainingj2009explaining and by hooker2019benchmark for their RemOve And Retrain (ROAR) method.
Table (ref) shows a small example that we will use to illustrate cohort Shapley. It has three binary predictors, two data distributions (one dependent and one independent) and two response functions (one additive and one with interaction). The dependent data distribution features an impossible combination. We can think of a setting where one binary variable specifies whether somebody has any children and a second whether their first child is male. Then the $(0,1)$ combination is impossible. Or we can think of a setting for the Boston housing data describing a region with many large lots and also high crime which while not logically impossible is quite improbable. In the given example, we have made a three variable combination impossible. We have also left variable three out of $f_1$. Variable 3 could be something like the race of a person which is not used by the algorithm at all, but might have some importance from association with variables that are used.
We begin with some notation. We will consider a data set with $n$ observations indexed by $i\in1{:}n\equiv\{1,\dots,n\}$. For observation $i$ there is a $d$-tuple $\boldsymbol{x}_i=(x_{i1},\dots,x_{id})$ of predictor variables. The component variables $x_{ij}$ do not always have to be real-valued so $\boldsymbol{x}_i$ is not necessarily a vector in $\mathbb{R}^d$. The data set also includes a response $y_i$ for observation $i$. The models we consider predict $y_i$ from $\boldsymbol{x}_i$. We focus on the case $y_i\in\mathbb{R}$. Sometimes we omit the observation's index and speak of predicting $y$ from $\boldsymbol{x}$. In some contexts we care about explaining the predictions for one specific target who then has index $t$.
In our formal discussions, we will need some special notation to study the changes that methods make to $\boldsymbol{x}$. For $u\subseteq1{:}d$ we let $\boldsymbol{x}_u$ be the tuple of variables $x_j$ for $j\in u$. We will use $|u|$ for the cardinality of $u$ and $-u=1{:}d\setminus u$ for the complement of $u$. For tuples $\boldsymbol{x}$ and $\tilde\boldsymbol{x}$ we let $\boldsymbol{x}_u{:}\tilde\boldsymbol{x}_{-u}$ be the hybrid point $\boldsymbol{x}'$ with $\boldsymbol{x}'_j=\boldsymbol{x}_j$ for $j\in u$ and $\boldsymbol{x}'_j=\tilde\boldsymbol{x}_j$ for $j\not\in u$. We abbreviate $\boldsymbol{x}_{\{j\}}$ and $\boldsymbol{x}_{-\{j\}}$ to $\boldsymbol{x}_j$ and $\boldsymbol{x}_{-j}$, respectively.
Our goal is to quantify the importance of the individual variables in a model that predicts $y$ from $\boldsymbol{x}$. Importance can be interpreted several ways. For instance, changing $x_{ij}$ could have a causal impact on $y_i$, removing predictor $x_j$ from a model to predict $y$ could degrade the accuracy of those predictions, or changing $x_{ij}$ could bring a large change to the prediction $\hat y_i$ jian:owen:2003. These are distinct issues about how $x_j$ relates to $y$. We are most interested in understanding specific decisions made by a given algorithm. For instance, which variables best explain why an algorithm may have recommended against lending money to a target applicant $t$, or may have recommended in favor of sending patient $t$ to an intensive care unit? For problems of this type, it is the third measure of importance that we want to study. In particular, an algorithm might have placed significant importance on a predictor known to have no causal effect, and we would want to detect that by showing how important that variable was to the predictions. It is for this goal that we especially want a method that avoids using impossible input combinations. In our conclusions we consider the consequences of including impossible combinations when judging whether variable $j$ is important for accurate prediction.
The models we consider predict $y$ by $\hat y = f(\boldsymbol{x})\in\mathbb{R}$. Importance of a variable can be `local', meaning that we want to explain $\hat y_t=f(\boldsymbol{x}_t)$ for a target point $t\in1{:}n$. It can also be global, meaning that we want to quantify the importance of variables in aggregate for a set of points.
To measure the importance of $x_j$ to $f(\boldsymbol{x})$ it is usual to change $x_j$, so $\boldsymbol{x}$ becomes $\tilde\boldsymbol{x}$, and to record the corresponding change $f(\tilde\boldsymbol{x})-f(\boldsymbol{x})$. There are an enormous number of ways to do this. We have to choose the value or values of $x_j$ to start from, and the value or values to change it to. Then we have to decide whether or not to also change $x_{j'}$ for $j'\ne j$ while we are changing $x_j$, and if we choose to change $x_{j'}$ we must specify how. When multiple changes are under consideration, we must find a way to aggregate them and then compare or rank the aggregate measures for variables $j\in1{:}d$.
Here are a few concrete examples from among many choices. Classical sensitivity analysis based on $|\partial f(\boldsymbol{x}_0)/\partial x_j|$ starts with a default point $\boldsymbol{x}_0\in\mathbb{R}^d$ and considers infinitesimally small changes there, assuming that $f$ is differentiable and the $x_j$ have been properly scaled relative to each other. The method of morr:1991 makes small (local) changes of size $\Delta_j$ to $x_j$ from a fine grid of starting points for $x_j$ with $x_{j'}$ for $j'\ne j$ sampled independently and uniformly over their values. It looks at the mean, standard deviation and even the whole cumulative distribution function of the resulting $f(\tilde\boldsymbol{x})-f(\boldsymbol{x})$ values. camp:cari:salt:2007 propose to aggregate $|f(\tilde\boldsymbol{x})-f(\boldsymbol{x})|$ instead. Global sensitivity analysis raza:etal:2021, about which we say more below, considers not just local changes but all possible changes from random $\boldsymbol{x}$ to random $\tilde\boldsymbol{x}$ with some but not all $\tilde x_j=x_j$, aggregating the changes by second moments. The variable importance measure in random forests of brei:2001 makes a random permutation of $x_{1j},\dots,x_{nj}$ with respect to the other components of $\boldsymbol{x}_1,\dots,\boldsymbol{x}_n$ and records how prediction accuracy changes in response to permuting the $j$th variable.
There is not room to survey all variable importance methods and we must therefore leave many of them out and focus on the ones most closely related to our proposed method. wei:lu:song:2015 provide a comprehensive survey of variable importance measures in statistics, organized around a taxonomy. They cite 197 references of which 24 are themselves surveys. The global sensitivity analysis literature is focused on climate models, computer aided design models and similar tools that often have a strong physical sciences emphasis. In addition to the recent survey of raza:etal:2021, with over 350 references, there are also textbooks of salt:etal:2008 and dave:gamb:ioos:prie:2021. moln:2018 surveys the explainable AI (XAI) literature. Two prominent methods discussed there are LIME ribe:etal:2016 and SHAP lund:lee:2017. A survey of methods based on Shapley value appears in sund:najm:2019, which appeared as we were completing mase:owen:seil:2019.
Our objective is to study $f$ to explain it, and this is easier than studying a ground truth quantity like $\mathbb{E}(Y\!\mid\! \boldsymbol{x})$ because $f$ is in our computer while $\mathbb{E}(Y\!\mid\! \boldsymbol{x})$ must ordinarily be estimated by gathering real-world data. A question about what changing $x_j$ to $\tilde x_j$ would do to $f(\cdot)$ can be directly answered by computing $f(\boldsymbol{x}_{-j}{:}\tilde\boldsymbol{x}_j)-f(\boldsymbol{x})$ for any setting $\boldsymbol{x}_{-j}$ of the other $d-1$ variables. Problems about causal inference when directed towards studying which variables cause $f$ to change are amenable to direct computational solutions. On the other hand, explaining why $f(\boldsymbol{x})=1$ instead of $f(\boldsymbol{x})=0$ is harder. There could be several $j\in1{:}d$ that would have given rise to $f(\boldsymbol{x}_{-j}{:}\tilde \boldsymbol{x}_j)=0$ had we changed $x_j$ to some other value $\tilde x_j$. In some problems every $j$ could have changed $f$. We might also have to consider the effects of changing more than one component of $\boldsymbol{x}$ at a time. We then have to consider which changes to one or more components $\boldsymbol{x}$ are most worthy of consideration.
The variable importance problem is one of studying `causes of effects' and not `effects of causes'. This distinction has been made repeatedly by Paul Holland; see, for example, holl:1986,holl:1988 citing mill:1851 and other philosophers. Holland's view is that statistical analysis is better suited to studying effects of causes. There is an excellent discussion of the difference between causes of effects and effects of causes in dawi:musi:2021 that also considers legal uses, such as a case about attributing a person's development of type II diabetes to past consumption of Lipitor. They consider at length the difficulties of using potential outcomes rubi:1974 or structural causal models pear:2009 in settings where there is only one variable being considered as the possible cause of some effect. In Section 16, they mention briefly that the setting with multiple putative causes brings further difficulties. That is precisely the problem one faces in variable importance settings. The following quote is from their conclusions:
Identifying causes of effects is a problem of attribution. We seek to measure the extent to which different variables $x_j$ have contributed to the value of $f(\boldsymbol{x})$. The big differences between attribution methods come down to their different definitions, which have different strengths and weaknesses. We cannot expect a `bakeoff' where algorithms compete to estimate some known true attributions to settle a discussion about which definition is most appropriate. This means that the common task framework dono:2019 that has been so effective in improving prediction methods has limited utility for attribution. A second issue that Holland raises about `causes of effects' is that the cause we identify might itself be an effect of some other cause. In causal inference, we can avoid adjusting for variables whose values became known after an intervention. For attribution, it is not so clear how to separate proximate causes from some prior causes. We might also have to think of even earlier ultimate causes of those prior causes.
We are familiar with the problem of attribution in real life. If a candidate lost an election by $10{,}000$ votes and some issue caused an unfavorable vote swing of $11{,}000$ votes, then we could at first reasonably attribute her loss to that issue. However, there may be ten such issues that individually or in combination brought a swing of over $10{,}000$ votes. It is then less clear that the first issue is `the reason'. The attribution issue has been intensely studied in advertising. In media mix modeling (MMM), users are told what percentage of their sales to attribute to each advertising channel; see chan:perr:2017 for a discussion of challenges in MMM. In online advertising, Shapley values are now being used instead of simply attributing a sale or other customer conversion to the last click prior to that event berm:2018. Shapley values are also being used in financial profit-and-loss attribution moeh:boyd:ang:2021.
One maxim in causal inference is that there is `no causation without manipulation'; see, for instance, holl:1986. Under this maxim, if we cannot in principal intervene to change a variable then we cannot make claims about its causal effect, at least not via the potential outcomes framework where both a treatment and an alternative must be possible. Judea Pearl disagrees with this maxim. For instance bollen2013eight remark that this principal would lead to a conclusion that the moon does not cause tides. In problems of attribution, we cannot avoid discussing the importance of variables, such as a specific person's birth year or ancestry, that an investigator cannot intervene to change.
Now we turn to quantifying the effect of $x_j$ on $f(\boldsymbol{x})$. The setting is simplest for an additive model of independent random inputs. Suppose that
where $x_j$ has support $\mathcal{X}_j$ and $x_1,\dots,x_d$ are independent random quantities with $f_j:\mathcal{X}_j\to\mathbb{R}$. Global importance measures can then be defined in terms of semi-norms of the $f_j$. We could choose $\mathrm{var}(f_j(x_j))$ or $\mathrm{var}(f_j(x_j))^{1/2}$ or $\sup_{x_j\in\mathcal{X}_j}f_j(x_j)-\inf_{x_j\in\mathcal{X}_j}f_j(x_j)$. If $\mathcal{X}_j\subseteq\mathbb{R}$ and $f_j$ has derivative $f_j'$, then we could choose $\sup_{x_j\in\mathcal{X}_j}|f'(x_j)|$ or the total variation $\int_{\mathcal{X}_j}|f'(x_j)|\,\mathrm{d} x_j$ or $\mathbb{E}(|f'(x_j)|)$. Semi-norms make sense because replacing $f_j(x_j)$ by $f_j(x_j)+c$ and $f_{j'}(x_{j'})$ by $f_{j'}(x_{j'})-c$ does not reasonably change the importance of $x_j$ and $x_{j'}$. Implicit in each of these measures is a choice of how to combine the effects of changes in $f_j$ and hence in $f$ over pairs $x_j,\tilde x_j\in\mathcal{X}_j$. These various semi-norms quantify importance in different ways, but the setting of equation (ref) is the one where it is most straightforward to choose a definition for a given problem for global and also local variable importance. The effect of changing $x_j$ is unrelated to the value of $x_{j'}$ and given a value for $x_j$ we have no reason to direct our attention about $x_{j'}$ to any specific part of $\mathcal{X}_{j'}$.
Variable importance becomes much more difficult when the additivity condition (ref) does not hold. In that case the variables have interactions, and the importance of those interactions has to be shared somehow among the variables that contribute. A second complication is that the distribution of $\boldsymbol{x}$ strongly affects which variable is important. To illustrate this effect, note that according to the CDC \url{https://www.cdc.gov/measles/about/faqs.html} about 3% of people immunized against measles will get measles if they are exposed to it. Let's assume that anybody not exposed does not get measles and also, since it is very contagious, that 90% of people not immunized will get it if exposed. Table (ref) illustrates how an interaction can complicate the relative importance of variables. If almost everybody is in the upper left cell, then exposure is the more important variable. If almost everybody is in the lower right cell then immunity is the more important variable. Note that this second complication from interactions is present already in the setting of independent input variables. For a thorough account of the phenomenon of interaction, see cox:1984 and dego:cox:2007.
While interactions can be conveniently handled by Sobol' indices based on the ANOVA decomposition as we will see in Section (ref), dependence among the components of $\boldsymbol{x}$ is far harder to deal with. In that case, the reasonable changes to $x_j$ can depend on all components in $\boldsymbol{x}_{-j}$ and their combinations. The Boston housing data illustration in Figure (ref) shows a setting with two variables where at least one must be near a default value. An extreme version of this problem is given in owen:prie:2017 where $\boldsymbol{x}\in\{0,1\}^2$ but the only possible values of $\boldsymbol{x}$ are $\{0,0\}$, $\{0,1\}$ and $\{1,0\}$. If $\boldsymbol{x}=\{0,1\}$ then it is not meaningful to change $x_1$ while leaving $x_2$ unchanged, so such a change cannot contribute to the importance of variable $1$. Higher dimensional patterns will be much harder to detect.
The challenges brought by dependence are not easy to handle via the ANOVA. See chas:gamb:prie:2012 and owen:prie:2017. When there are dependencies among variables, we can turn from ANOVA to the Shapley value that we describe in Section (ref). Shapley value opens up some conceptually attractive solutions that nonetheless have computational challenges such as estimating an exponentially large number of conditional expectations.
Our approach of studying changes to $f$ does not capture every variable importance measure in widespread use. For instance the importance of predictors in spike-and-slab regression models is often measured by the posterior probability that $\beta_j\ne0$ which is not directly obtainable from $f(\boldsymbol{x})=\widehat\mathbb{E}(y\!\mid\! \boldsymbol{x})$ scot:vari:2014. Such measures provide a useful aggregate for the role of $x_j$ but do not easily explain local decisions. Similarly, counts of the number of times a variable is used to define a split in a tree or forest model do not have a direct interpretation in terms of the value of $f$.
In this section we describe some of the tools we need in order to present the cohort Shapley measure. These are general analysis of variance (ANOVA) decomposition, Sobol' indices, and Shapley value.
When the components of $\boldsymbol{x}$ are independent and $\mathbb{E}(f(\boldsymbol{x})^2)<\infty$, there is a convenient analysis of variance (ANOVA) decomposition of $f$. This ANOVA can be used to define variance components and subsequently Sobol' indices. Sobol' indices provide a very powerful way to quantify global importance of variables. They do not provide local explanations, and, as we will see below, they have difficulty with dependent data settings.
The familiar ANOVA used in experimental design applies to tabular data defined in terms of categorical $x_j$. It goes back to fish:mack:1923. A generalization to $\boldsymbol{x}\sim\mathsf{U}[0,1]^d$ was used by hoef:1948, sobo:1969 and efro:stei:1981 for U-statistics, integration and the jackknife, respectively. The ANOVA decomposes $\sigma^2=\mathrm{var}(f(\boldsymbol{x}))$ into components for each $u\subseteq1{:}d$. For $d=\infty$ see lss.
To understand the generalized ANOVA, we recall how the orginal one works. We are given $f(\boldsymbol{x})$ where now $x_j$ takes $k_j$ levels and we have all $N=\prod_{j=1}^dk_j$ possible evaluations of $f$. The $x_j$ are then independent random variables under equal weighting of those $N$ values. We first compute a grand mean $\mu$ averaging all the values $f(\boldsymbol{x}_1),\dots,f(\boldsymbol{x}_N)$. Next, to define the main effect of variable $j$, we subtract $\mu$ from each $f(\boldsymbol{x}_i)$ and at each level of $x_j$ we average those differences over $N/k_j$ settings of the other $d-1$ variables. The idea is that only differences merit an explanation. For an interaction between variables $j$ and $j'$, we subtract the grand mean and both main effects from $f(\boldsymbol{x}_i)$ and for each of the $k_jk_{j'}$ combinations of $x_j$ and $x_{j'}$ we average the results over all $N/(k_jk_{j'})$ settings of the other $d-2$ variables. The general rule is like this: we don't attribute to $\boldsymbol{x}_u$ what can be explained by $\boldsymbol{x}_v$ for $v\subsetneq u$, then to get a function of $\boldsymbol{x}_u$ we average over $\boldsymbol{x}_{-u}$.
The general ANOVA replaces averages by expectations. It begins with a grand mean function $f_{\varnothing}(\boldsymbol{x})$ which is everywhere equal to the constant $\mu\equiv\mathbb{E}( f(\boldsymbol{x}))$. The main effect for variable $j$ is the function $$f_{\{j\}}(\boldsymbol{x}) = \mathbb{E}( f(\boldsymbol{x})-\mu\!\mid\! \boldsymbol{x}_j)=\mathbb{E}(f(\boldsymbol{x})\!\mid\! \boldsymbol{x}_j)-\mu.$$ For a subset of variables $u\subseteq1{:}d$, we define
The first expression in (ref) shows that we do not attribute to $\boldsymbol{x}_u$ what can be explained by $\boldsymbol{x}_v$ for any proper subset $v$ of $u$. The second expression follows because $f_v$ is not random given $\boldsymbol{x}_u$. While $f_u$ is defined on the whole domain of $\boldsymbol{x}$ its value only depends on $\boldsymbol{x}_u$.
When $|u|>1$, the effect $f_u$ is a $|u|$-fold interaction. The effect $f_u$ has variance component $\sigma^2_u =\mathrm{var}(f_u(\boldsymbol{x}))$. We easily find that $$ f(\boldsymbol{x}) = \sum_{u\subseteq1{:}d}f_u(\boldsymbol{x})\quad\text{and}\quad \sigma^2=\sum_{u\subseteq1{:}d}\sigma^2_u. $$
Sobol' indices sobo:1993 measure the importance of the subvector $\boldsymbol{x}_u$ not just individual variables $x_j$. There are two versions, the closed index and the total index, respectively, $$ \underline\tau^2_u = \sum_{v\subseteq u}\sigma^2_v\quad\text{and}\quad\overline\tau^2_u = \sum_{v:v\cap u\ne{\varnothing}}\sigma^2_v. $$ These are usually divided by $\sigma^2$ to give a proportion of variance explained. We easily see that $0\leqslant \underline\tau^2_u \leqslant \overline\tau^2_u = \sigma^2-\underline\tau^2_{-u}$. The total index $\overline\tau^2_u$ includes interactions between $\boldsymbol{x}_u$ and $\boldsymbol{x}_{-u}$ that the closed index $\underline\tau^2_u$ excludes. If $\underline\tau^2_u$ is relatively large then we interpret $\boldsymbol{x}_u$ as important while if $\overline\tau^2_u$ is small then $\boldsymbol{x}_u$ can be interpreted as unimportant. One can show that $\underline\tau^2_u =\mathrm{var}( \mathbb{E}(f(\boldsymbol{x})\!\mid\! \boldsymbol{x}_u))$ and $\overline\tau^2_u = \mathbb{E}(\mathrm{var}(f(\boldsymbol{x})\!\mid\! \boldsymbol{x}_{-u}))$.
Sobol' indices have been intensely studied in the global sensitivity analysis literature. They are very conveniently estimated by sampling methods based on identities like $$ \underline\tau^2_u = \mathbb{E}\bigl( f(\boldsymbol{x})f(\boldsymbol{x}_u{:}\boldsymbol{z}_{-u})\bigr)-\mu^2 \quad\text{and}\quad \overline\tau^2_u = \frac12\mathbb{E}\bigl((f(\boldsymbol{x})-f(\boldsymbol{x}_{-u}{:}\boldsymbol{z}_{u}))^2\bigr) $$ due to sobo:1993 and jans:1999, respectively.
Sobol' indices provide a very good global measure of how important a variable or even a set of variables is. They are not aimed at variable importance at a specific point $\boldsymbol{x}$, so they are not suited to local importance.
We mentioned earlier that the ANOVA does not extend well to dependent data settings. Since Sobol' indices are based on the ANOVA, this causes difficulties for them. The ANOVA for dependent data has been considered by chas:gamb:prie:2012 building on work of hook:2007 who builds on that of ston:1994. One of the challenges is that some ways to define ANOVA lead to negative analogues of variance components for dependent data. A key condition in some of the definitions is that the density $p(\boldsymbol{x})$ must be bounded below by $c\times p_u(\boldsymbol{x}_u)\times p_{-u}(\boldsymbol{x}_{-u})$ for some $c>0$ and all $\boldsymbol{x}$. This rules out distributions with `holes' where $p(\boldsymbol{x})=0$ on a rectangle within the support of $p$. It also rules out multivariate Gaussian distributions with nonzero correlations.
The two aforementioned problems with Sobol' indices do not apply to Shapley value which we consider next.
The Shapley value shap:1953 is used in cooperative game theory to define a fair allocation of rewards to team members who have jointly produced some value. It has seen many uses in defining variable importance measures. See sund:najm:2019 and moln:2018 for surveys of uses in XAI and song:nels:stau:2016 in global sensitivity analysis. The analogy is to view the variables $x_j$ in a model $f(\boldsymbol{x})$ as team members cooperating on a goal such as predicting $y$.
In Shapley's formulation, there is a set $1{:}d$ of entities. Together they produce a value that we denote by $\nu(1{:}d)$. We suppose that we can find the value $\nu(u)$ that would have been produced by any $u\subseteq1{:}d$. It is customary to assume that $\nu({\varnothing})=0$. Shapley proposes four quite reasonable axioms on the basis of which there is then a unique value $\phi_j=\phi_j[\nu]$ that should be attributed to player $j$ on the team. A key quantity in Shapley's formulation is the incremental value that would be brought to a team $u$ if player $j\not\in u$ were to join it. This incremental value of player $j$ given that players in $u$ are already on the team is $$ \nu(j\!\mid\! u) = \nu( u\cup\{j\})-\nu(u). $$
Shapley's four axioms are:
The second and third axioms are compelling. moln:2018 suggests that the first axiom might even be a requirement in settings where one is required by law to provide an explanation. See also doshi2017accountability. The fourth axiom is more debatable but is usually included in variable importance applications and we will use it. Under these four axioms there is a unique solution
An intuitive explanation of (ref) is as follows: Suppose that we build a team from ${\varnothing}$ to $1{:}d$ in one of $d!$ possible orders and track the incremental value that arises at the moment when player $j$ joins this team. Then $\phi_j$ is the average of those $d!$ incremental values. Actually averaging all permutations would have a cost that is $O(d!)$. The direct calculation of (ref) depends on only $d2^{d-1}$ incremental values: each of $2^d$ sets $u$ has $d$ neighboring sets that differ by inclusion or deletion of exactly one of the $d$ team members. Every pair $(u,j)$ with $u\subset1{:}d$ and $j\not\in u$ is counted twice this way.
For illustration, suppose that $d=3$ as in our running example. Then there are $6$ permutations of the three variables. Variable $1$ enters first in two of them, following ${\varnothing}$. It enters last in two of them, following $\{2,3\}$. It enters once after just $\{2\}$ is present and once after just $\{3\}$ is present. Therefore
and the other $\phi_j$ are similar.
The Shapley values $\phi_j$ are unchanged if we add a constant $c$ to all of the values $\nu(u)$. This follows because $c$ cancels out of the incremental values. As remarked above, it is common to assume that $\nu({\varnothing})=0$. Sometimes it is simpler to define a value function in a way that does not force $\nu({\varnothing})=0$. In such cases the Shapley values attribute the incremental value $\nu(1{:}d)-\nu({\varnothing})$ of the whole team to the members of that team.
The cost to compute Shapley values exactly grows exponentially in $d$. The representation via permutations opens up a natural Monte Carlo strategy of using an average over a sample of randomly generated permutations. Monte Carlo is routinely used in settings such as integration where an exact computation would have infinite cost. Its effectiveness depends on the variance of the incremental values and not on the number $d2^{d-1}$ of incremental values needed for exact calculation. michalak2013efficient and mitchell2022sampling describe some value functions where they find that Monte Carlo is inefficient. For cohort Shapley, Monte Carlo would be unfavorable if only a vanishing fraction of the incremental values $\nu(j\!\mid\! u)$ were non-negligible. We consider this unlikely and have not seen it in our experience.
The integrated gradients method of sundararajan2017axiomatic is based on the Aumann-Shapley value of auma:shap:1974. Where the Shapley value averages difference over paths on the edges of $\{0,1\}^d$, integrated gradients average a gradient over the diagonal of $[0,1]^d$ which can be much faster to compute.
In this section, we show how Shapley values have been used to study variable importance. Then we present the cohort Shapley algorithm that only uses actually observed data for variable importance. We also include a discussion of some tradeoffs and choices that come up when using Shapley values for variable importance. There are many different ways to formulate a game with which to apply the Shapley formulations. The axioms have considerable logical force, conditional on the value function being a reasonable one. As a result, the important differences among methods based on Shapley values stem from the choice of value function.
lind:mere:gold:1980 propose a measure of the importance of including a variable in a linear regression model. Their proposal is equivalent to using the Shapley value with $\nu(u)$ is proportion of the variance explained by a linear model (usually denoted $R^2$) using the variables $\boldsymbol{x}_u$. stru:kono:2010 used Shapley value to explain the value of $f(\boldsymbol{x}_t)$ in a local explanation for target observation $t$. In their equation (1), they use the value function $$ \nu(u) = \mathbb{E}( f(\boldsymbol{x}_{t,u}{:}\boldsymbol{x}_{-u}))-\mathbb{E}(f(\boldsymbol{x})) $$ where the first expectation is taken with respect to independent discrete random variables $\boldsymbol{x}_j$ for $j\not\in u$ and the second one has all $\boldsymbol{x}_j$ from independent discrete distributions. Those discrete distributions could be the empirical marginal distributions. They also consider random sampling approximation that allows continuous distributions.
For $\boldsymbol{x}$ with independent components, owen:2014 takes $\nu(u)$ to be the variance explained by $\boldsymbol{x}_u$, $\mathrm{var}( \mathbb{E}(f(\boldsymbol{x})\!\mid\! \boldsymbol{x}_u))=\underline\tau^2_u$. The result is that $$ \phi_j = \sum_{u\subseteq1{:}d} \frac{1\{j\in u\}\sigma^2_u}{|u|}. $$ Every variance component $\sigma^2_u$ is shared equally among the $|u|$ variables that contribute to it. It follows that the easily estimated Sobol' indices bracket the Shapley value: $$ \underline\tau^2_j\leqslant \phi_j\leqslant \overline\tau^2_j.$$ plis:rabi:borg:2021 note that the upper bound can be improved to $(\overline\tau^2_j+\underline\tau^2_j)/2$.
song:nels:stau:2016 propose the use of Shapley value with $\nu(u)=\mathrm{var}(\mathbb{E}(f(\boldsymbol{x})\!\mid\! \boldsymbol{x}_u))$ for a setting where components of $\boldsymbol{x}$ may be dependent. They present a computational algorithm based on permutations and they show that the same Shapley values arise using $\nu(u)=\mathbb{E}(\mathrm{var}(f(\boldsymbol{x})\!\mid\! \boldsymbol{x}_{-u}))$. owen:prie:2017 describe how Shapley values solve the conceptual problems that dependence causes for Sobol' indices, though computational challenges remain.
Baseline Shapley (see sund:najm:2019) explains $f(\boldsymbol{x}_t)-f(\boldsymbol{x}_b)$ where $t$ is a target point and $\boldsymbol{x}_b$ is a `baseline' vector of predictors. It could be made of default values or it could be $\bar\boldsymbol{x} =(1/n)\sum_{i=1}^n\boldsymbol{x}_i$, although some components of $\bar\boldsymbol{x}$ might not be possible feature values. In baseline Shapley, the value function is $\nu(u)=f(\boldsymbol{x}_{t,u}{:}\boldsymbol{x}_{b,-u})$ and unobserved hybrid values $\boldsymbol{x}_{t,u}{:}\boldsymbol{x}_{b,-u}$ are to be expected.
Many of the alternatives to baseline Shapley use a `conditional expectation Shapley' with $$\nu(u) =\mathbb{E}( f(\boldsymbol{x}_t)\!\mid\! \boldsymbol{x}_{t,u})$$ differing according to which distribution of $\boldsymbol{x}_{t,-u}$ given $\boldsymbol{x}_{t,u}$ they use to define the expectation. For instance the above mentioned choice of stru:kono:2010 is of this type assuming independent predictors. The quantitative input influence of datt:2016 takes $\nu(u)=\mathbb{E}(f(\boldsymbol{x}_u{:}\boldsymbol{z}_{-u})\!\mid\!\boldsymbol{x}_u)$ under a sampling model where $\boldsymbol{z}_{-u}$ has the empirical marginal distribution of $\boldsymbol{x}_{-u}$. That will generally produce unobserved hybrid points. Kernel SHAP lund:lee:2017, similarly uses an independence assumption on inputs to compute a conditional Shapley value.
The TreeSHAP method lund:2020 is designed to be especially computationally efficient on tree structured models but can only be used on such models. If the prediction is a weighted sum of tree predictions then, using additivity of the Shapley value one can take a weighted sum of Shapley values of single trees. For a single tree, the value is $\mathbb{E}(f(\boldsymbol{x})\!\mid\!\boldsymbol{x}_u)$ where the distribution of $f(\boldsymbol{x}_{-u})$ given $\boldsymbol{x}_u$ is uniform over all the values of $f(\boldsymbol{x}_i)$ in all the leaves of the tree that are compatible with $\boldsymbol{x}_u$. It only uses observed data and it cannot attribute importance to variables that are not used by the model.
frye:2021 share our concerns about using impossible data. They use a variational autoencoder (VAE) to represent the $d$ predictors as being scattered around a lower dimensional manifold. They compute $\mathbb{E}(f(\boldsymbol{x})\!\mid\! \boldsymbol{x}_u)$ using their model to define the distribution of $\boldsymbol{x}_{-u}$ given $\boldsymbol{x}_u$, but not every point $\boldsymbol{x}$ they use will have been observed, so whether those are reasonable depends on performance of the chosen VAE. Similarly, aas:2019 propose a conditional kernel SHAP that uses weights from an ellipsoidally shaped kernel based on a Mahalanobis distance between target and reference points using the sample covariance. They need a scaling parameter for each subset of input variables. Their kernel puts positive weight on unobserved variable combinations.
Finally, we note that ghor:zou:2019 use Shapley value to measure the value of individual data points in a model. They compare this to leave-one-out methods and leverage scores and find that it helps to identify outliers and corrupted data. For instance, data points with a negative Shapley value for model accuracy are worthy of inspection; perhaps the response value was mislabeled or there is some other error in that data.
Some more uses of Shapley value in XAI are surveyed in sund:najm:2019, moln:2018 and kuma:2020.
The cohort Shapley method quantifies variable importance using only actually observed data values. Our description here is based on mase:owen:seil:2019.
For target $t$ and variable $j$ we define a function $s_{t,j}$ where $s_{t,j}(\boldsymbol{x}_{ij})=1$ if $\boldsymbol{x}_{ij}$ is similar to $\boldsymbol{x}_{tj}$ and is zero otherwise. It must always be true that $s_{t,j}(\boldsymbol{x}_{tj})=1$. For binary variables, similarity is just whether $\boldsymbol{x}_{ij}=\boldsymbol{x}_{tj}$. For real valued features $j$, similarity could be that $|\boldsymbol{x}_{ij}-\boldsymbol{x}_{tj}|\leqslant\delta$ or $|\boldsymbol{x}_{ij}-\boldsymbol{x}_{tj}|\leqslant \delta |\boldsymbol{x}_{tj}|$ for some $\delta>0$, or we could discretize $\boldsymbol{x}_{ij}$ values into a set of ranges and take observations to be similar if they are in the same range. This last choice makes similarity transitive, which the other choices do not impose.
For observation $t$ and $u\subseteq1{:}d$, define the cohort
with $C_t({\varnothing})=1{:}n$ by convention. The value function in cohort Shapley is $$ \nu(u) =\frac1{|C_t(u)|}\sum_{i\in C_t(u)}f(\boldsymbol{x}_i), $$ for local explanation for observation $t$. These values are all well defined because $t\in C_t(u)$ always holds. In the random permutation framework, we begin with the cohort $C_t({\varnothing})$ of all observations whose mean is then the global mean $(1/n)\sum_{i=1}^nf(\boldsymbol{x}_i)$. We then refine the cohort, imposing similarity with respect to all $d$ variables in order. An important variable is one that brings a relatively large change to the cohort mean when we refine on it. The cohort Shapley values $\phi_j$ explain the difference between the mean of $f(\boldsymbol{x}_i)$ over the fully refined cohort $C_t(1{:}d)$ and the global mean. If all input variables have been binned, then cohort Shapley becomes a conditional expectation Shapley defined in terms of the binned variables.
In applying cohort Shapley, one is forced to define what it means for $\boldsymbol{x}_{ij}$ to be similar to $\boldsymbol{x}_{tj}$. This is difficult for continuous variables. Modest numbers of levels are usually recommended in similar settings that require `binning' a continuous variable. gelm:park:2009 use means over the largest and smallest of just three groups to approximate a regression slope with small loss of efficiency. coch:1968 finds diminishing returns to using more than five strata in sampling problems. Cochran's results are based on linear regression and nonlinear relationships would reasonably require a few more strata. For a modestly nonlinear setting like propensity scores for logistic regression, neuh:thie:ruxt:2018 find diminishing returns at around ten strata.
While binning a variable can seem arbitrary we think that it is more transparent than employing yet another black box, such as a VAE, to approximate empirical conditional distributions.
We can understand TreeSHAP in terms of cohorts. For a single tree it uses cohorts defined by a union of leaves of that tree. Such a union of leaves need not define a spatially connected set when $\mathcal{X} =\mathbb{R}^d$. For aggregates of multiple trees the notion of similarity will generally differ between the trees in that aggregate. Our concern with TreeSHAP is that it uses a notion of variable similarity defined in part by the response values it is fitting. This makes it harder to interpret or explain the underlying similarity concept. In a high stakes setting, such as whether a given algorithm was unfair to one or more people, we would be interested to know who those people were compared to. It is also not able to attribute importance to a variable that the model does not use.
Categorical variables with many levels may need to be coarsened similarly to continuous variables. If all of the variables are categorical or have been made categorical by stratification then cohort Shapley matches conditional expectation Shapley based on the empirical distribution of the data.
Table (ref) shows the cohorts of target observation $1$ for our running example. Potential observation 7 is unobserved in one of the data distributions and so it is removed from the two cohorts that it would have appeared in had all $8$ $\boldsymbol{x}$ values been observed. That table also records the cohort means for our two distributions and two example functions, $f_1(\boldsymbol{x}) = x_1 - 2 x_2$ and $f_2(\boldsymbol{x}) = x_1 x_2 (1 - x_3)$.
Table (ref) shows cohort Shapley values for our running example. For independent data and $f_1(\boldsymbol{x})=x_1-2x_2$ the average value of $f(\boldsymbol{x}_i)$ over all data (cohort $C_1({\varnothing})$) is $\nu({\varnothing})=-1/2$. For target observation $t=1$ and cohort $C_1(\{1,2,3\})=\{1\}$ we get $\nu({\{1,2,3\}})=0$. Therefore we have to explain why $f(\boldsymbol{x}_t)$ is $1/2$ unit above average. There are six equally weighted paths from ${\varnothing}$ to $\{1,2,3\}$. As described at equation (ref), we get
Similarly, $\phi_2=1$ and $\phi_3=0$. So the target observation $t=1$ has $f$ that is $1/2$ unit above average and the cohort Shapley attributions are that this arises from a $-1/2$ effect of having $x_1=0$ offset by a $+1$ effect of having $x_2=0$ while $x_3$ does nothing. We see that $\phi_j$ takes account of not just the sign of the coefficient of $x_j$ but also whether $x_{tj}$ is at the low or high value of $x_j$.
If $\boldsymbol{x}=(1,1,0)$ is impossible, as in distribution $P_*$, then that introduces some amount of dependence between $x_3$ and the other variables. Now $\phi_3\ne0$ for the target observation and the linear function $f_1=x_1-2x_2$ that does not involve $x_3$. For $f_2=x_1x_2(1-x_3)$ under the independence distribution $P$, the target observation has $f_2(\boldsymbol{x}_t)=0$ while the average of $f_2$ is $1/8$. The Shapley attributions are $-1/12$ for both $x_1$ and $x_2$ and $+1/24$ for $x_3$ owing to the factor $1-x_3$. Under $P_*$, $f_2$ is only nonzero at the impossible combination and then all $\phi_j=0$.
There are several difficult choices to consider in rating variable importance generally and some difficult choices specific to methods based on Shapley value. These stem from `importance' meaning different things in different contexts. kuma:2020 present several qualms about the use of Shapley value in any of its forms for variable importance.
One consequential decision is whether to attribute any importance to a variable $x_j$ that is not used at all in $f(\boldsymbol{x})$. In a case like this, changes to $x_{tj}$ cannot possibly change the prediction $\hat y_t = f(\boldsymbol{x}_t)$ so it might seem unreasonable to attribute importance to $x_{tj}$. Suggesting that person $t$ change $x_{tj}$ to qualify for a loan would be bad advice. However there are settings in algorithmic fairness where a protected variable such as race or gender that is not used by $f(\cdot)$ might still play a role through its association with other variables that were used. In those applications, it is essential to consider the possibility that $x_j$ could be important. Fairness by unawareness is not considered reliable dwor:2012. Cohort Shapley can be used to identify importance for variables not used by the function. To do that, one includes that variable in the set of predictors under study. We illustrate this in Section (ref), in an example involving race.
A potential problem with Shapley values arises when two variables $x_j$ and $x_{j'}$ are very strongly related. This can cause them to `share their importance' which might make both of them appear unimportant. Or, as kuma:2020 note, if one variable denotes a protected category and the other does not, including both in the Shapley calculations could make the protected variable look less important than it really is.
Algorithmic fairness makes use of variable importance. For instance, it is unfair when a protected variable such as race or gender is improperly important in a decision that affects a person. In this section, we conduct a detailed investigation of the data from angw:2016 regarding the fairness of the COMPAS tools and their use in recidivism prediction. The focal issue is whether COMPAS is unfair to Black defendants, and if so, by how much and to which defendants? The results in this section were presented earlier in fairnessshapley, which this article supercedes.
There are several aspects of cohort Shapley that make it very suitable to analyzing fairness of the COMPAS data:
Cohort Shapley is the unique method from Table (ref) to satisfy condition 2 and one of only two that satisfy condition 1. The uncertainty in condition 6 is distinct from uncertainty due to random sampling of permutations. In this analysis we are able to compute cohort Shapley exactly and the uncertainty arises from viewing the defendants as a sample from a larger population. This setting is simpler to bootstrap because the model was trained on a separate data set.
Just as there are many ways to define variable importance, there are multiple ways to define what fairness means. Some of those definitions are mutually incompatible and some of them differ from legal definitions. Here we present just a few of the issues as a prelude to studying the COMPAS data. See corb:goel:2018, chou:roth:2018, berk:2018, and frie:etal:2019 for surveys. We do not make assertions about which definitions are preferable.
For $y,\hat y\in\{0,1\}$, let $n_{y\hat y}$ be the number of observations with $y_i=y$ and $\hat y_i=\hat y$. These four counts and their derived properties can be computed for any subset of the observations. The false positive rate (FPR) is $n_{01}/n_{0\text{\tiny$\bullet$}}$, where a bullet indicates that we are summing over the levels of that index. We ignore uninteresting corner cases such as $n_{0\text{\tiny$\bullet$}}=0$; when there are no observations with $y=0$ then we have no interest in the proportion of them with $\hat y=1$. The false negative rate (FNR) is $n_{10}/n_{1\text{\tiny$\bullet$}}$. The prevalence of the trait under study is $p=n_{1\text{\tiny$\bullet$}}/n_{\text{\tiny$\bullet$}\text{\tiny$\bullet$}}$. The positive predictive value (PPV) is $n_{11}/n_{\text{\tiny$\bullet$}1}$. As chou:2017 notes, these values satisfy
See also klei:mull:ragh:2016.
Equation (ref) shows how some natural definitions of fairness conflict. FPR and FNR describe $\hat y\!\mid\! y$, while PPV describes $y\!\mid\!\hat y$. If two subsets of observations have the same PPV but different prevalences $p$, then they cannot also match up on FPR and FNR. Fairness in $y\!\mid\! \hat y$ terms and fairness in $\hat y\!\mid\! y$ terms can only coincide in trivial settings such as when $\hat y=y$ always or empirically unusual settings with equal prevalence between observations having different values of a protected variable.
There is some debate about when or whether using protected variables can lead to improved fairness. See xian:2020 who gives a summary of legal issues surrounding fairness and corb:etal:2017 who study whether imposing calibration or other criteria might adversely affect the groups they are meant to help.
datt:2016 attributes demographic parity of the prediction $\hat y$ to input variables through aggregated Shapley value. We quantify individual level bias on protected variables using cohort Shapley and then aggregate to a measure for any group of interest. Cohort Shapley can measure quantities of various fairness definitions that are conventionally quantified by group level statistical measures.
The COMPAS Core Risk and Needs Assessment tool from Northpointe Inc.\ includes the General Recidivism Risk Scale (GRRS) and Violent Recidivism Risk Scale (VRRS) which were investigated by angw:2016. A more complete description of these scales can be found in bren:etal:2009 and nort:2019. We will focus on the GRRS data collected and published by angw:2016 in the analysis that follows when referring to the “COMPAS data”. Each defendant is rated into one of ten deciles with higher deciles considered higher risk of reoffending. angw:2016 investigated whether that algorithm was biased against Black people. They obtained data for defendants in Broward County Florida, including the COMPAS decile, the defendants' race, age, gender, number of prior arrests and whether the crime for which they were charged is a felony or not. angw:2016 describe how they processed their data including how they found followup data on offences committed and how they matched those to defendants for whom they had prior COMPAS scores. They also note that race was not one of the variables used in the COMPAS predictions, and the original model for those predictions was trained on more features than those available and included in this Broward County data. angw:2016 look at pre-trial defendants and measure recidivism as being arrested for another crime within two years, and we will use the same metrics and data for our analysis so that it can be compared to the original work. This example is controversial. For a more complete discussion of the appropriateness of using this data to evaluate COMPAS including the important discrepancies between pre-trial vs probation or parole use and the time frame for reoffending, see rudi:etal:2020, flor:etal:2016, diet:etal:2016, and jack:mend:2020. While this data is therefore likely inadequate for the purposes of evaluating the true merits of COMPAS, we include it as an example here given its place in the literature as a known benchmark for meaningful comparison.
angw:2016 find that COMPAS is biased because it gave a higher rate of false positives for Black defendants and a higher rate of false negatives for White defendants. flor:etal:2016 and diet:etal:2016 disagree, raising the issue of $\hat y\!\mid\! y$ fairness versus $y\!\mid\! \hat y$ fairness. The prevalence of reoffences differed between Black and White defendants forcing $y\!\mid\!\hat y$ and $\hat y \!\mid\! y$ notions to be incompatible. tan:etal:2018 find some evidence of racial bias when comparing an interpretable surrogate model trained to predict $\hat{y}$ to one trained on $y$, but stop short of attributing that bias to COMPAS itself instead pointing to the difference between features used to train COMPAS and those present in the data set. agar:etal:2021 also find evidence of racial bias in an interpretable surrogate model trained to predict COMPAS scores. rudi:etal:2020 find no evidence that the COMPAS scores depend on race except through age and criminal history. Further they find that the same racial FPR/FNR discrepancies detected by angw:2016 when looking at the COMPAS scores would exist when using a model that relied solely on age.
Our cohort Shapley analysis below finds statistically significant racial effects. White defendants tend to have higher than predicted rates of recidivism while Black defendants tend to have lower than predicted rates of recidivism. The magnitudes of those effects are in the range of 3.5 to 5.5 percent when measured by the average of $y-\hat y$. Consistent with some of the earlier findings we see larger effects of race on the false positive and false negative rates.
Beyond these aggregate findings, our analysis provides local (i.e., individual level) racial effects. This would be useful in looking at a specific case of interest or, as we show below, for finding that the impact of race varies by gender. While some previous work like angw:2016 and rudi:etal:2020 have examined specific individuals or pairs of individuals as examples, this work is the first to present a method that quantifies this local fairness systematically.
Following chou:2017 we focus on just Black and White defendants. That provides a sample of 5278 defendants from among the original 6172 defendants. As in that paper we record the number of prior arrests as a categorical variable with five levels: 0, 1--3, 4--6, 7--10 and $>$10. Following angw:2016, we record the defendants' ages as a categorical variable with three levels: $<$25, 25--45 and $>$45. Also, following chou:2017, we consider the prediction $\hat y_i$ to be $1$ if defendant $i$ is in deciles 5--10 and $\hat y_i=0$ for defendant $i$ in deciles 1--4.
Figure (ref) shows some conventional group fairness metrics for the COMPAS data set. Horizontal bars show group specific means and the vertical dashed lines show population means. We see that Black defendants had a higher average value of $\hat y$ than White defendants. Black defendants also had a higher average of $y$ but a lower average residual $y-\hat y$. Using $B$ and $W$ to denote the two racial groups, $\hat\mathbb{E}(y-\hat y\!\mid\! B)\doteq -0.035$ and $\hat\mathbb{E}(y-\hat y\!\mid\! W)\doteq 0.054$. The FPR was higher for Black defendants and the FNR was higher for White defendants.
In this section we use cohort Shapley to study some fairness issues in the COMPAS data, especially individual level metrics. We have selected what we found to be the strongest and most interesting findings related to race and gender from fairnessshapley. Histograms there also show effects due to age and the number of prior arrests.
We have computed Shapley impacts for these responses: $y_i$, $\hat y_i$, $y_i-\hat y_i$, $\mathrm{FP}_i = 1\{ y_i=0\,\&\, \hat y_i=1\}$ and $\mathrm{FN}_i = 1\{ y_1=1\,\&\, \hat y_i=1\}$. If $\mathrm{FP}_i=1$ then defendant $i$ received a false positive prediction. Note that the sample average value of $\mathrm{FP}_i$ $$\hat\mathbb{E}(\mathrm{FP}_i)=\frac{n_{01}}{n_{\text{\tiny$\bullet$}\text{\tiny$\bullet$}}} =\mathrm{FPR}\times \frac{n_{0\text{\tiny$\bullet$}}}{n_{\text{\tiny$\bullet$}\text{\tiny$\bullet$}}} =\mathrm{FPR}\times (1-p), \text{ and } \hat\mathbb{E}(\mathrm{FN}_i)=\mathrm{FNR}\times p$$
Figure (ref) shows histograms of Shapley impacts of race for the defendants in the COMPAS data. The first panel there shows a positive impact for every Black defendant and a negative one for every White defendant for the prediction $\hat y$. For the actual response $y$, the histograms overlap slightly. By additivity of Shapley value the impacts for $y-\hat y$ can be found by subtracting the impact for $\hat y$ from that for $y$ for each defendant $i=1,\dots,n$. The histograms of $y_i-\hat y_i$ show that the impact of race on the residual is typically positive for White defendants and negative for Black defendants.
The histograms of Shapley impacts for race overlap for the two metrics, FPR and FNR. There, a small number of adversely affected White defendants and beneficially affected Black defendants are observed. We can inspect the corresponding defendants to see which conditions are the exceptional cases. For example, the 469'th defendant (Black) and the 486'th defendant (White) have the same values for all other features. They both had 4-6 prior arrests, age $<$25, crime severity of felony, and both were female. In Section (ref) we give an example with a confidence interval for the Shapley values of an individual defendant. The same could be done for a difference between two defendant's Shapley values.
Figure (ref) shows histograms of Shapley impacts for race color-coded by the defendants' gender. We see that the impact on $\hat y$ is bimodal by race for both male and female defendants, but the effect is larger in absolute value for male defendants. The impacts for the response, the residual and false positive values do not appear to be bimodal by race for female defendants, but they do for male defendants. The impacts for false negatives do not appear bimodal for either race. It is clear from these figures that the race differences we see are much stronger among male defendants.
We take a particular interest in the residual, or error, $y_t-\hat y_t$. It equals $1$ for false negatives and $-1$ for false positives and $0$ when the prediction was correct. Table (ref) shows race and gender Shapley values for the residual $y_i-\hat y_i$ aggregated over all four race-gender subsets of defendants.
What we see there is that revealing that a defendant is Black tells us that, on average, that defendant's residual $y_t-\hat y_t$ is decreased by $3.6$%. By this measure, the Black defendants offend less than predicted. Revealing that the defendant is White increases the residual by $5.4$%. Revealing race makes very little difference to the residual averaged over male or over female defendants (of both races). Revealing gender makes a relatively large difference of $-8.2$% for female defendants and $+2.0\%$ for male defendants.
To judge the uncertainty in the values in Table (ref) we applied the Bayesian bootstrap of rubi:1981. That algorithm randomly reweights each data point by a unit mean exponential random variable. There is an asymptotic justification in newton1994approximate. Figure (ref) shows violin plots of 1000 bootstrapped cohort Shapley values. We see that the mean Shapley value of race on the residual is significantly negative for Black defendants and significantly positive for White defendants. These differences persist when disaggregated by gender. The average effect of race is near zero for male defendants and for female defendants due to race effects nearly cancelling.
Because the Bayesian bootstrap never fully deletes any of the data we can use it to quantify the statistical uncertainty in the cohort Shapley values for an individual defendant. Figure (ref) shows a Bayesian bootstrap violin plot of cohort Shapley values on the residual for the 2999'th defendant. This is the Black defendant whose Shapley value for race on the residual $y-\hat y$ was most negative, in the region overlapping Shapley values of White defendants. As one would expect there is greater uncertainty on the impacts for an individual than for an aggregate. Our analysis has not taken account of the fact that this defendant was selected based on their data.
The risk of being falsely predicted to reoffend involves two factors: having $y_i=0$ and having $\hat y_i=1$ given that $y_i=0$. FPR is commonly computed only over defendants with $y_i=0$. Accordingly, in this section we study it by subsetting the defendants to $\{ i\mid y_i=0\}$ and finding cohort Shapley value of $\hat y_i=1$ for the features. That provides a different way to quantify the importance of race on the incidence of false positives.
Figure (ref) shows the cohort Shapley impact of race conditioned for the FPR and FNR after working with subsets of defendants as described above. The distribution of impacts conditioned on race are clearly separated in the figure. This conditional analysis shows a stronger disadvantage for Black defendants than in Figure (ref).
There are many incompatible definitions of importance. Local importance measures address a causes of effects question, which means that we have to carefully consider which counterfactuals to use and how to use them. When attributing the value of $f(\boldsymbol{x}_i)$ to features $x_{ij}$, we have avoided using any impossible data by only using actually observed data. We think it is important to have such a choice as otherwise the attribution could be based on extremely unlikely or even impossible combinations. hook:ment:2019:tr also point to problems with mixing and matching variable inputs.
There could be good uses for possible but hitherto unobserved input combinations. Using only observed values means that we rely on data cleaning to have eliminated any impossible combinations. An alternative that uses input values that pass a `within distribution' test could be used but in our view it is exceedingly difficult to draw the line between possible and impossible values based on a sample. kuma:2020 also raise this point.
There may be a use for impossible data when explaining whether a predictor affects accuracy. For instance, when judging whether a predictor helps predictions, brei:2001 studies what happens when $x_{1j},\dots,x_{nj}$ are subjected to a random permutation $\pi(\cdot)$ with respect to the other variables, producing hypothetical values $\tilde\boldsymbol{x}_i=\boldsymbol{x}_{i,-j}{:}\boldsymbol{x}_{\pi(i),j}$. He then looks at the accuracy of predicting $y_i$ by $f(\tilde\boldsymbol{x}_i)$. If accuracy does not change much, then it is reasonable to suppose that $x_j$ has only a minimal role in predictions. If accuracy does change a lot however, it could be because the fitted model extrapolates poorly into the low probability or impossible regions of the space. Breiman's method appears then to be a `one way diagnostic', persuasive when it indicates low importance but otherwise not.
Despite these possibilities, we adopted a constraint of not using any unobserved values in order to define a measure of importance. As mentioned in the introduction this counters adversarial approaches that can make changes to predictions at unobserved points in order to change the fairness measures of an algorithm.
In defining cohort Shapley we have required a human in the loop to define a similarity measure for the variables. We acknowledge that this is a burden. However, the alternatives we have seen use additional black boxes or use the response itself to define similarity and we think those choices make the final result harder to understand. We think that people would not accept an automatic algorithmic declaration that some black box is or is not fair without having the ability to see how it defines the similar cases that should be treated similarly. Also, the human or humans making and discussing those choices can decide on similarity for one variable at a time. Those choices then define the conditional expectations needed by the algorithm.
In addition to transparency and avoiding impossible data, cohort Shapley brings some other advantages: the user can opt to measure importance of variables that were not used in the model, the method does not require access to the black box beyond the predictions it makes, and having defined cohorts we can study means or medians or other quantities of interest. It can also be applied to data on which the algorithm was not trained, as for instance with the COMPAS data. This allows one to quantify external validity properties directly. All conditional expectation Shapley value methods use estimates of conditional means and these become challenging as the conditioning dimension increases and the sample size decreases. For a setting where the black box is used on fresh data that it was not trained on, we have illustrated a convenient bootstrap method to quantify the uncertainty that arises from sampling observations.
With the COMPAS data we have illustrated how cohort Shapley values can be used to explore data surrounding a complex fairness issue. There are many more analyses that could be done, and we expect that domain experts would not be unanimous on which is best. We do not have the data necessary to compare fairness of COMPAS to fairness of human judges.
This work was supported by Hitachi, Ltd. and by the National Science Foundation grant IIS-1837931. We thank Sharad Goel, Jessica Hwang, Alexandra Chouldechova and Alice Xiang for helpful discussions. The opinions in this article are our own and not necessarily shared by those we acknowledge.