remarkRemark hypothesisHypothesis claimClaim Derivative-based Shapley valueHui
Duan and Giray Ökten ex_supplement

# Derivative-based Shapley value for global sensitivity analysis and machine
learning explainability

Hui Duan Department of Mathematics, Florida State University, Tallahassee FL,
USA<EMAIL_ADDRESS>Giray Ökten Department of Mathematics, Florida
State University, Tallahassee FL, USA<EMAIL_ADDRESS>
###### Abstract

We introduce a new Shapley value approach for global sensitivity analysis and
machine learning explainability. The method is based on the first-order
partial derivatives of the underlying function. The computational complexity
of the method is linear in dimension (number of features), as opposed to the
exponential complexity of other Shapley value approaches in the literature.
Examples from global sensitivity analysis and machine learning are used to
compare the method numerically with activity scores, SHAP, and KernelSHAP.

###### keywords:

Shapley value; activity scores; derivative-based global sensitivity measures;
Sobol’ sensitivity indices; global sensitivity analysis; machine learning
explainability

##  1 Introduction

Global sensitivity analysis is the study of how uncertainty in the output of a
model can be allocated to uncertainties in the model input. The applications
of sensitivity analysis cover a wide range of disciplines, from natural
sciences to engineering, and social sciences to mathematical sciences. Several
methods for global sensitivity analysis have been introduced in the
literature, including Sobol’ sensitivity indices, derivative-based global
sensitivity measures (DGSM), activity scores, and Shapley values. In this
paper, we introduce a derivative-based Shapley value (DerSHAP) approach for
global sensitivity analysis and machine learning explainability.

The paper is organized as follows. In Section 2 we discuss Sobol’ sensitivity
indices, derivative-based global sensitivity measures, and activity scores. In
Section 3 we review the Shapley value, and in Section 4 we introduce our
method, the derivative-based Shapley value (DerSHAP). In Section 5 we use
DerSHAP to perform global sensitivity analysis of an Ebola model and an HIV
model, and compare it numerically with the activity score approach. We also
apply DerSHAP to two datasets from machine learning, Amazon stock data and
Boston housing data, and compare it with SHAP and KernelSHAP. Conclusions
follow in the last section.

##  2 Global Sensitivity Analysis

In this section, we review some of the background material from global
sensitivity analysis. We will give a review of Sobol’ sensitivity indices,
derivative based global sensitivity measures (DGSM), and activity scores.

###  2.1 Sobol’ Sensitivity Indices

Consider a $d$ dimensional input vector $\textbf{x}=(x_{1},\ldots,x_{d})$ with
the index set $D=\\{1,2,\ldots,d\\}$. Let $f(\textbf{x})$ be a square-
integrable function defined on $(0,1)^{d}$. The ANOVA decomposition of
$f(\textbf{x})$ is expressed as

$f(\textbf{x})=\sum_{u\subseteq D}f_{u}(\textbf{x}^{u}),$

where $f_{u}(\textbf{x}^{u})$ is the component function that only depends on
$\textbf{x}^{u}$. For the empty set, we have $f_{\emptyset}=\int
f(\textbf{x})d\textbf{x}$.

If we assume x has uniform distribution on $(0,1)^{d}$, we can write

$\mathbb{E}[f(\textbf{x})]=\int_{(0,1)^{d}}f(\textbf{x})d\textbf{x},$

and

$\text{Var}(f(\textbf{x}))=\sigma^{2}=\int_{(0,1)^{d}}f^{2}(\textbf{x})d\textbf{x}-\mathbb{E}[f(\textbf{x})]^{2}.$

Due to the orthogonality of the ANOVA decomposition, the variance of $f$ can
be written as

$\sigma^{2}=\sum_{u\subseteq D}\sigma_{u}^{2},$

where $\sigma_{u}^{2}$ is the variance for the component function $f_{u}$:

$\sigma_{u}^{2}=\int_{(0,1)^{d}}f_{u}^{2}(\textbf{x})d\textbf{x}-\left(\int_{(0,1)^{d}}f_{u}(\textbf{x})d\textbf{x}\right)^{2}=\int_{(0,1)^{d}}f_{u}^{2}(\textbf{x})d\textbf{x}.$

The Sobol’ sensitivity indices for the subset $u$ are defined from the
variances of component functions

$\underline{S}_{u}=\frac{1}{\sigma^{2}}\sum_{v\subseteq
u}\sigma_{v}^{2}=\frac{\underline{\tau}_{u}}{\sigma^{2}}\text{ and
}\overline{S}_{u}=\frac{1}{\sigma^{2}}\sum_{v\bigcap
u\neq\emptyset}\sigma_{v}^{2}=\frac{\overline{\tau}_{u}}{\sigma^{2}},$

where $\underline{S}_{u}$ is called the lower Sobol’ sensitivity index (or,
the main effect) and $\overline{S}_{u}$ is called the upper Sobol’ sensitivity
index (or, the total effect). If $\underline{S}_{u}$ is close to 1, then the
parameters $\textbf{x}^{u}$ are viewed as very important to the model. In
practical applications, we often consider the case when $u$ is a singleton
$\\{i\\}$. If $\overline{S}_{\\{i\\}}$ is close to 0, we consider $x_{i}$ to
be an unimportant variable, and freeze it to its mean value to simplify the
model $f$.

Sobol’ sensitivity indices can be generalized to independent random vectors
$\textbf{x}=(x_{1},\ldots,x_{d})$ with any distribution functions
$F_{1}(x_{1}),\ldots,F_{d}(x_{d})$. See Kucherenko et. al (2012) [6] for
details.

###  2.2 Derivative Based Global Sensitivity Measures

If the partial derivatives of $f$ exist, computationally more efficient global
sensitivity measures can be obtained. For example, Campolongo et al. (2007)
[1] introduced a global sensitivity measure based on
$\int_{(0,1)^{d}}\left|\frac{\partial f(\textbf{x})}{\partial
x_{i}}\right|d\textbf{x}$, and Sobol’ & Kucherenko (2009) [17] introduced a
global sensitivity measure based on $\int_{(0,1)^{d}}(\frac{\partial
f(\textbf{x})}{\partial x_{i}})^{2}d\textbf{x}$. Here we will discuss the
latter measure further. The derivative based global sensitivity measure (DGSM)
of Sobol’ & Kucherenko (2009) [17] is defined as

(1) $v_{i}=\int_{(0,1)^{d}}(\frac{\partial f(\textbf{x})}{\partial
x_{i}})^{2}d\textbf{x}=\mathbb{E}[(\frac{\partial f(\textbf{x})}{\partial
x_{i}})^{2}].$

DGSM can be estimated by Monte Carlo as

$v_{i}\approx\hat{v}_{i}=\frac{1}{N}\sum_{j=1}^{N}(\frac{\partial
f(\textbf{x}^{(j)})}{\partial x_{i}})^{2},$

where $\textbf{x}^{(j)},j=1,\ldots,N$ is a random sample from the uniform
distribution on $(0,1)^{d}$.

Sobol’ (2011) [16] showed that if $f(\operatorname{\textbf{x}})$ is a linear
function on each of its components $x_{i}$, then

(2) $\overline{S}_{i}=\frac{1}{12}\frac{v_{i}}{\sigma^{2}}.$

In general, the following inequality was obtained in Sobol’ & Kucherenko
(2009) [17]

(3) $\overline{S}_{i}\leq\frac{1}{\pi^{2}}\frac{v_{i}}{\sigma^{2}}.$

DGSM can be generalized to nonuniform distributions, where the random vector
$\operatorname{\textbf{x}}=(x_{1},\ldots x_{d})$ has independent components
with marginal distributions $F_{1},\ldots,F_{d}$, and each $F_{i}$ has a
corresponding density function $f_{i}$. Sobol’ & Kucherenko (2009) [17] and
Kucherenko and Iooss (2014) [5] obtained several results in this generalized
setting. In particular, the following generalization of inequality (3) was
obtained in [5]:

(4) $\overline{S}_{i}\leq D(F_{i})\frac{v_{i}}{\sigma^{2}},$

where

$D(F_{i})=4\left[\sup_{x\in\mathbb{R}}\frac{\min(F_{i}(x),1-F_{i}(x))}{f_{i}(x)}\right]^{2}.$

###  2.3 Activity Scores

The activity score was introduced by Constantine and Diaz (2017) [2], and it
is based on the active subspace method (Constantine et. al (2014) [3]). The
active subspace method finds the important directions of a function using the
gradient information, and uses these important directions to reduce the domain
of the function to a subspace. The activity score is a global sensitivity
measure obtained from this information. We follow [2] to describe the method.

Consider a square-integrable function $f(\operatorname{\textbf{x}})$ defined
on $(-1,1)^{d}$, with finite partial derivatives. Assume the partial
derivatives of $f$ are square-integrable. Let $p(\operatorname{\textbf{x}})$
be the uniform probability density function on $(-1,1)^{d}$. The gradient of
$f$, $\nabla f(\operatorname{\textbf{x}})\in\mathbb{R}^{d}$, is given by

$\nabla f(\operatorname{\textbf{x}})=\left[\frac{\partial
f(\textbf{x})}{\partial x_{1}},\frac{\partial f(\textbf{x})}{\partial
x_{2}},\ldots,\frac{\partial f(\textbf{x})}{\partial x_{d}}\right]^{T}.$

Define the matrix $C$ as

(5) $C=\mathbb{E}[\nabla f(\textbf{x})\nabla f(\textbf{x})^{T}],$

where the expectation is computed using the density function $p$. Consider the
eigenvalue decomposition of $C$

$C=\textbf{W}\Lambda\textbf{W}^{T},$

where $\textbf{W}=[\textbf{w}_{1},\ldots,\textbf{w}_{d}]$ is the $d\times d$
orthogonal matrix of eigenvectors, and
$\Lambda=\text{diag}(\lambda_{1},\ldots,\lambda_{d})$ with
$\lambda_{1}\geq\ldots\geq\lambda_{d}\geq 0$ is the diagonal matrix of
eigenvalues in descending order. Matrix $C$ can be approximated using the
Monte Carlo method:

$C\approx\hat{C}=\frac{1}{N}\sum_{j=1}^{N}(\nabla_{\operatorname{\textbf{x}}}f(\operatorname{\textbf{x}}^{({j})}))(\nabla_{\operatorname{\textbf{x}}}(f(\operatorname{\textbf{x}}^{(j)}))^{T},$

where the
$\operatorname{\textbf{x}}^{(1)},\ldots,\operatorname{\textbf{x}}^{(N)}$ is a
random sample of size $N$ from the density $p$.

If the eigenvalues $\lambda_{m+1},\ldots,\lambda_{d}$ are sufficiently small,
then the active subspace method approximates $f(\textbf{x})$ with a lower
dimensional function $g(\textbf{W}_{1})$ where $\textbf{W}_{1}$ is the
$d\times m$ matrix containing the first $m$ eigenvectors. The dimension of $g$
is $m$, whereas the dimension of $f$ is $d$.

The activity score for the $i$th parameter is defined as

(6) $\alpha_{i}(m)=\sum_{j=1}^{m}\lambda_{j}w_{ij}^{2},$

where $\textbf{w}_{j}=[w_{1j},\ldots,w_{dj}]^{T}$ is the $j$th eigenvector,
$m\leq d$, and $i=1,\ldots,d$.

Constantine and Diaz (2017) [2] showed that the activity scores are bounded by
DGSM

$\alpha_{i}(m)\leq v_{i},\text{ }i=1,\ldots,d,$

and the inequality becomes an equality if $m=d$.

The following inequality between the Sobol’ sensitivity index and activity
scores was also established in [2]:

$\overline{S}_{i}\leq\frac{1}{4\pi^{2}}\frac{\alpha_{i}(m)+\lambda_{m+1}}{\sigma^{2}}.$

This result can be generalized to non-uniform distributions over $R^{d}$ using
Eqn. (4).

######  Lemma 2.1.

For any random vector $\operatorname{\textbf{x}}$ of independent components
with marginal distributions $F_{i}$ and density $f_{i}$, $i=1,\ldots,d$,
distributions $F_{1},\ldots,F_{d}$, we have

$\overline{S}_{i}\leq
D(F_{i})\frac{(\alpha_{i}(m)+\lambda_{m+1})}{\sigma^{2}}.$

##  3 Shapley Values

Shapley (2016) [15] introduced the Shapley value to solve the question of
distributing a reward among multiple team members fairly in cooperative game
theory. Owen (2014) [12] used the Shapley value approach to quantify the
importance of input variables of a function to the function output. This
approach is particularly valuable when the input variables are dependent (Owen
& Prieur (2017) [13]).

Consider a function $f$ of $d$ variables, $f(x_{1},\ldots,x_{d})$, and let
$\mathcal{D}=\\{1,2,\ldots,d\\}$ be the index set. Let imp be a function from
the powerset of $\mathcal{D}$ to positive real numbers, with the
interpretation that for any subset $u$ of $\mathcal{D}$, $\text{imp}(u)$ is a
measure of importance of input variables with indices in $u$, and
$\text{imp}(\emptyset)=0$. The question we want to answer is how much of the
importance of the whole set of input variables can be attributed to each input
variable. To this end, let $\phi_{i}$ denote the contribution of input
variable $x_{i}$ to $\text{imp}(\mathcal{D})$ (in the terminology of
cooperative game theory, this function is called the Shapley value). Shapley
showed that for any given function imp, there is a corresponding function
$\phi_{i}$ satisfying the following conditions:

  1. 1.

Efficiency: $\sum_{i=1}^{d}\phi_{i}=\text{imp}(\mathcal{D})$;

  2. 2.

Dummy: If $\text{imp}(u\cup\\{i\\})=\text{imp}(u)$ for all
$u\subseteq\mathcal{D}$, then $\phi_{i}=0$;

  3. 3.

Symmetry: If $\text{imp}(u\cup\\{i\\})=\text{imp}(u\cup\\{j\\})$ for all
$u\subseteq\mathcal{D}$ such that $u\cap\\{i,j\\}=\emptyset$, then
$\phi_{i}=\phi_{j}$;

  4. 4.

Additivity: If $\text{imp},\text{imp}^{\prime}$ are two measures of importance
with corresponding functions $\phi$, $\phi^{\prime}$, then the importance
measure $\text{imp}+\text{imp}^{\prime}$ has the corresponding function
$\phi_{i}+\phi_{i}^{\prime}$.

This unique solution is given by

(7)
$\phi_{i}=\frac{1}{d}\sum_{u\subseteq-\\{i\\}}{d-1\choose|u|}^{-1}(\text{imp}(u\cup\\{i\\})-\text{imp}(u)),$

where $|u|$ denotes the cardinality of $u$, and $-\\{i\\}$ is the complement
of $\\{i\\}$.

For black-box models, Lundberg and Lee (2017) [10] proposed the Shapley
additive explanation (SHAP). This method uses conditional expectations to
represent the imp function and has gained popularity in machine learning
(Broeck et al. (2022) [19]), despite some authors having raised questions
about its explainability (Kumar et al. (2020) [7]). Lundberg et al. (2018) [9]
introduced TreeSHAP to compute the Shapley values for tree-based machine
learning models with reduced complexity. Mase et al. (2019) [11] proposed the
cohort Shapley to measure the similarity of a specific subject to the target
subject. Seiler et al. (2021) [14] defined uniqueness Shapley based on the
size of the cohort. They aimed to interpret the importance of a variable in
identifying a specific subject from a set of data in the machine learning
models.

##  4 Shapley Value Based on Derivatives

In this section we introduce an importance function, $\text{imp}(u)$, for
subsets $u$ of $\mathcal{D}$ based on the partial derivatives of $f$, similar
in spirit, to the DGSM and activity score measures of global sensitivity
analysis. We assume that the first-order partial derivatives of $f$ are
finite, and square-integrable. To simplify the notation, we denote the partial
derivative of $f(\textbf{x})$ with respect to $x_{i}$ as

$\frac{\partial f(\textbf{x})}{\partial x_{i}}=\partial f_{i}.$

Define the derivative-based importance of subset $u$ as follows:

(8) $\text{imp}(u)=\sum_{i\in u}\sum_{j\in u,j\geq i}\left|\mathbb{E}[\partial
f_{i}\partial f_{j}]\right|=\sum_{i\in u}\mathbb{E}[\partial
f_{i}^{2}]+\sum_{i,j\in u,j>i}\left|\mathbb{E}[\partial f_{i}\partial
f_{j}]\right|.$

Note that the first summation of the right-hand side can be written in terms
of the activity scores

$\sum_{i\in u}\mathbb{E}[\partial f_{i}^{2}]=\sum_{i\in u}\alpha_{i}(d),$

and thus we can write

$\text{imp}(u)=\sum_{i\in u}\alpha_{i}(d)+\sum_{i,j\in
u,j>i}\left|\mathbb{E}[\partial f_{i}\partial f_{j}]\right|.$

For the empty set, we set $\text{imp}(\emptyset)=0$, and for the full index
set, $\text{imp}(\mathcal{D})$ is the sum of entries of the upper triangular
matrix of $|C|$ (see Eqn. (5)):

$\text{imp}(\mathcal{D})=\sum_{i=1}^{d}\sum_{j\geq
i}^{d}\left|\mathbb{E}[\partial f_{i}\partial f_{j}]\right|.$

######  Theorem 4.1.

The Shapley value for the derivative-based importance function is

$\phi_{i}=\mathbb{E}[\partial f_{i}^{2}]+\frac{1}{2}\sum_{j=1,j\neq
i}^{d}\left|\mathbb{E}[\partial f_{i}\partial f_{j}]\right|.$

######  Proof 4.2.

The efficiency condition follows from

$\displaystyle\sum_{i=1}^{d}\phi_{i}=\sum_{i=1}^{d}\mathbb{E}[\partial
f_{i}^{2}]+\sum_{i=1}^{d-1}\sum_{j=i+1}^{d}\left|\mathbb{E}[\partial
f_{i}\partial f_{j}]\right|=\text{imp}(\mathcal{D}).$

To prove the dummy condition, assume $\text{imp}(u\cup\\{i\\})=\text{imp}(u)$
for all subsets $u$. Also assume we can write $\mathcal{D}=v\cup w\cup\\{i\\}$
where $v,w$ are nonempty disjoint sets not containing $i$. Observe that
$\text{imp}(v\cup\\{i\\})=\text{imp}(v)$ implies

$\mathbb{E}[\partial f_{i}^{2}]+\sum_{k<i,k\in v}\left|\mathbb{E}[\partial
f_{k}\partial f_{i}]\right|+\sum_{k>i,k\in v}\left|\mathbb{E}[\partial
f_{i}\partial f_{k}]\right|=\mathbb{E}[\partial f_{i}^{2}]+\sum_{k\in
v}\left|\mathbb{E}[\partial f_{k}\partial f_{i}]\right|=0.$

Likewise, $\text{imp}(w\cup\\{i\\})=\text{imp}(w)$ implies

$\mathbb{E}[\partial f_{i}^{2}]+\sum_{k\in w}\left|\mathbb{E}[\partial
f_{k}\partial f_{i}]\right|=0.$

Adding the last two equations we obtain

$2\phi_{i}=2\mathbb{E}[\partial f_{i}^{2}]+\sum_{k\neq
i,k\in\mathcal{D}}\left|\mathbb{E}[\partial f_{k}\partial f_{i}]\right|=0.$

To prove symmetry assume
$\text{imp}(u\cup\\{k_{1}\\})=\text{imp}(u\cup\\{k_{2}\\})$ for all subsets
$u$ that do not contain $k_{1}$ and $k_{2}$. Assume we can write
$\mathcal{D}=v\cup w\cup\\{k_{1},k_{2}\\}$ where $v,w$ are nonempty disjoint
sets not containing $k_{1}$ and $k_{2}$. (If $v$ and $w$ are the empty set and
a singleton, the condition holds trivially.) Note that
$\text{imp}(v\cup\\{k_{1}\\})=\text{imp}(v\cup\\{k_{2}\\})$ implies

$\mathbb{E}[\partial f_{k_{1}}^{2}]+\sum_{i\in v}\left|\mathbb{E}[\partial
f_{k_{1}}\partial f_{i}]\right|=\mathbb{E}[\partial f_{k_{2}}^{2}]+\sum_{i\in
v}\left|\mathbb{E}[\partial f_{k_{2}}\partial f_{i}]\right|$

and $\text{imp}(w\cup\\{k_{1}\\})=\text{imp}(w\cup\\{k_{2}\\})$ implies

$\mathbb{E}[\partial f_{k_{1}}^{2}]+\sum_{i\in w}\left|\mathbb{E}[\partial
f_{k_{1}}\partial f_{i}]\right|=\mathbb{E}[\partial f_{k_{2}}^{2}]+\sum_{i\in
w}\left|\mathbb{E}[\partial f_{k_{2}}\partial f_{i}]\right|.$

Adding these equations we get $\phi_{k_{1}}=\phi_{k_{2}}$.

The derivative-based Shapley (DerSHAP) value involves the summation of $d$
terms where each summand is an expectation that can be approximated using
Monte Carlo. Therefore the complexity of the method is $O(dN)$ where $N$ is
the Monte Carlo sample size. In contrast, the original Shapley value
formulation requires a summation of $2^{d}$ terms (Eqn. (7)). The Shapley
values introduced by Owen & Prieur (2017)[13], Song et al. (2016) [18], and
Lundberg and Lee (2017) [10] all have exponential complexity. KernelSHAP
(Lundberg and Lee (2017) [10]) provides a more efficient algorithm than SHAP,
where it approximate SHAP values using linear regression. Nevertheless, it
still has exponential complexity. TreeSHAP(Lundberg et al. (2018) [9]) employs
the same important function as SHAP but has polynomial complexity $O(TLD^{2})$
for tree-based models, where $L$ is the number of leaves for an ensemble tree,
$T$ is the number of trees in the ensemble, and $D$ is the maximum tree depth.

If the eigenvalues $\lambda_{k+1},\ldots,\lambda_{d}$ of the matrix $|C|$
(Eqn. (5)) are small, then we can approximate the Shapley values by setting
the small eigenvalues to zero. Let the eigenvalue decomposition of $|C|$ be
given as $|C|=\textbf{W}\Lambda\textbf{W}^{T}$, and write

$|C|=\sum_{i=1}^{d}\lambda_{i}w_{i}w_{i}^{T},$

where $w_{i}$ are the orthonormal eigenvectors. Assume $\lambda_{i}<\epsilon$
for $i=k+1,\ldots,d$, and define

$\tilde{C}=\sum_{i=1}^{k}\lambda_{i}w_{i}w_{i}^{T},$

where $k\leq d$. The diagonal entries of $|C|-\tilde{C}$ are

$[\left|C\right|-\tilde{C}]_{j,j}=\sum_{i=k+1}^{d}\lambda_{i}w_{j,i}^{2}\leq\sum_{i=k+1}^{d}\lambda_{i}.$

The Frobenius norm of $|C|-\tilde{C}$ satisfies

$\|\left|C\right|-\tilde{C}\|_{F}^{2}=\|\sum_{i=k+1}^{d}\lambda_{i}w_{i}w_{i}^{T}\|_{F}^{2}=\sum_{i=k+1}^{d}\lambda_{i}^{2}\leq(d-k)\epsilon^{2}.$

Let $\Phi$ be the vector of derivative-based Shapley values, that is,
$\Phi=[\phi_{1},\ldots,\phi_{d}]^{T}$. Let $e$ be the column vector of 1’s. We
have

$\Phi=\frac{1}{2}|C|e+\frac{1}{2}\operatorname{diag}(|C|)e.$

Define $\tilde{\Phi}$ as the vector of “approximate” Shapley values, where the
matrix $|C|$ is replaced by $\tilde{C}$:

$\tilde{\Phi}=\frac{1}{2}\tilde{C}e+\frac{1}{2}\operatorname{diag}(\tilde{C})e.$

The error between $\Phi$ and $\tilde{\Phi}$ is

$\Phi-\tilde{\Phi}=\frac{1}{2}(|C|-\tilde{C})e+\frac{1}{2}\operatorname{diag}(|C|-\tilde{C})e,$

and can be bounded by

$\displaystyle\|\Phi-\tilde{\Phi}\|_{2}$
$\displaystyle=\|\frac{1}{2}(|C|-\tilde{C})e+\frac{1}{2}\operatorname{diag}(|C|-\tilde{C})e\|_{2}$
$\displaystyle\leq\frac{1}{2}\|(|C|-\tilde{C})e\|_{2}+\frac{1}{2}\|\operatorname{diag}(|C|-\tilde{C})e\|_{2}$
$\displaystyle\leq\frac{1}{2}\sqrt{d-k}\epsilon\sqrt{d}+\frac{1}{2}(d-k)\epsilon\sqrt{d}$
$\displaystyle\leq(d-k)\epsilon\sqrt{d},$

where $k\leq d$, and in the second inequality we used the fact
$\|(|C|-\tilde{C})e\|_{2}\leq\|\left|C\right|-\tilde{C}\|_{F}\|e\|_{2}$.

##  5 Numerical Results

In this section we apply the derivative-based Shapley value (DerSHAP) to
perform global sensitivity analysis of an Ebola model and an HIV model. We
will compare the method with the activity scores approach in these
applications. In both applications the model parameters have independent
random inputs, and thus the activity scores approach is applicable. We will
then apply DerSHAP to two datasets from machine learning, Amazon stock data
and Boston housing data, and compare DerSHAP with SHAP and KernelSHAP. In
these problems the input parameters are not independent. In all the numerical
results, we report the normalized activity scores and normalized Shapley
values. Normalization simply means dividing each value (activity score or
Shapley value) by the sum of all values. The sum of the normalized activity
scores and the sum of the normalized Shapley values are equal to one.

SHAP assumes feature independence, and KernelSHAP assumes feature independence
and model linearity. In the numerical examples, we use SHAP Python package
(2021)111https://github.com/slundberg/shap to implement both SHAP and
KernelSHAP.

###  5.1 The Ebola Model

This is a modified SEIR model for the spread of Ebola in Western Africa. A
detailed description of the model is given in Diaz et al. (2018) [4] where the
activity score method was used to perform the global sensitivity analysis of
the model.

The SEIR model is a system of seven differential equations which describe the
dynamics of disease within a population. We are interested in estimating the
basic reproduction number, which is a metric measuring how many new cases of
disease each current case causes. Using the steady state solution and the next
generation matrix, Diaz et al. (2018) [4] shows that the basic reproduction
number $R_{0}$ is

$R_{0}=\left(\beta_{1}+\frac{\beta_{2}\rho_{1}\gamma_{1}}{\omega}+\frac{\beta_{3}}{\gamma_{2}}\psi\right)/\left(\gamma_{1}+\psi\right).$

Table 1, which is from [4], displays the model parameters and their
distributions, obtained using data from Liberia and Sierra Leone.

Parameter | Liberia | Sierra Leone
---|---|---
$\beta_{1}$ | $U(.1,.4)$ | $U(.1,.4)$
$\beta_{2}$ | $U(.1,.4)$ | $U(.1,.4)$
$\beta_{3}$ | $U(.05,.2)$ | $U(.05,.2)$
$\rho_{1}$ | $U(.41,1)$ | $U(.41,1)$
$\gamma_{1}$ | $U(.0276,.1702)$ | $U(.0275,.1569)$
$\gamma_{2}$ | $U(.081,.21)$ | $U(.1236,.384)$
$\omega$ | $U(.25,.5)$ | $U(.25,.5)$
$\psi$ | $U(.0833,.7)$ | $U(.0833,.7)$
Table 1: Parameters of Ebola model

We use the activity scores and DerSHAP to compute the importance of each
parameter. The matrix $C$, which is used by both methods, is approximated with
Gauss-Legendre quadrature with eight points in each of the eight dimensions of
the parameter space as done in [4]. Fig. 1 plots the activity scores and
DerSHAP values for the model parameters for the Liberia and Sierra Leone
datasets.

(a) Activity scores - Liberia

(b) Activity scores - Sierra Leone

(c) DerSHAP values - Liberia

(d) DerSHAP values - Sierra Leone

Figure 1: Normalized activity scores and DerSHAP values for the Ebola model

Although the methods agree on the importance of some parameters, there are
some significant differences. DerSHAP values for $\beta_{2},\rho_{1},\omega$
are much larger than their activity scores for both data sets. Parameter
$\beta_{3}$ is less important than $\gamma_{2}$ according to the activity
scores for Sierra Leone data, but they are nearly equally important according
to DerSHAP. These differences are due to the fact that DerSHAP takes into
account the 2nd order interactions $\partial f_{i}\partial f_{j}$ for $i\neq
j$, whereas the activity score only considers $\partial f_{i}^{2}$. In terms
of computational complexity, the activity score method computes $d$ terms for
$d$ inputs, where each term is an expectation, approximated by numerical
quadrature (as in this problem) or Monte Carlo, whereas DerSHAP requires the
computation of $d$ expectations for one single input and $d+d(d-1)/2$
expectations for all $d$ inputs.

###  5.2 HIV Model

Here we consider the Human Immunodeficiency Virus (HIV) model in Loudon and
Pankavich (2016) [8]. The HIV model is described by the following differential
equations:

$\displaystyle\frac{dT}{dt}$
$\displaystyle=s_{1}+\frac{p_{1}}{C_{1}+V}TV-\delta_{1}T-(K_{1}V+K_{2}M_{I})T,$
$\displaystyle\frac{dT_{I}}{dt}$
$\displaystyle=\psi(K_{1}V+K_{2}M_{I})T+\alpha_{1}T_{L}-\delta_{2}T_{I}-K_{3}T_{I}CTL,$
$\displaystyle\frac{dT_{L}}{dt}$
$\displaystyle=(1-\psi)(K_{1}V+K_{2}M_{I})T-\alpha_{1}T_{L}-\delta_{3}T_{L},$
$\displaystyle\frac{dM}{dt}$ $\displaystyle=s_{2}+K_{4}MV-
K_{5}MV-\delta_{4}M,$ $\displaystyle\frac{dM_{I}}{dt}$
$\displaystyle=K_{5}MV-\delta_{5}M_{I}-K_{6}M_{I}CTL,$
$\displaystyle\frac{dCTL}{dt}$
$\displaystyle=s_{3}+(K_{7}T_{I}+K_{8}M_{I})CTL-\delta_{6}CTL,$
$\displaystyle\frac{dV}{dt}$
$\displaystyle=K_{9}T_{I}+K_{10}M_{I}-K_{11}TV-(K_{12}+K_{13})MV-\delta_{7}V,$

where $T(t)$ is the CD4+ T-cell population, $T_{I}$ is the actively infected
T-cell population, $T_{L}$ represents latently-infected T-cells, $M$ is
macrophages, $M_{I}$ is infected macrophages, $CTL$ is cytotoxic lymphocytes,
and $V$ is virions. This model has 27 input parameters

${\small\textbf{x}=(s_{1},s_{2},s_{3},p_{1},c_{1},s_{1},K_{1},K_{2},\cdots,K_{13},\delta_{1},\delta_{2},\cdots,\delta_{7},\alpha_{1},\psi).}$

The output of the model is the T-cell population, $T(t)$. The distributions of
the parameters are given in Loudon and Pankavich (2016) [8]. The authors
considered T-cell count at many time points, and used the active subspace
method to obtain a reduced model at each time. Here we will fix one particular
time, $t=24$ days, and apply the activity scores and DerSHAP to identify the
important input parameters.

Fig. 2 plots the activity scores and DerSHAP values for the HIV model. Both
methods agree on the most important inputs: $(K_{1},K_{9},\delta_{7},\psi)$.
However, DerSHAP values for many other input variables are nonzero. In
particular, the parameter $\delta_{2}$ has an activity score of less than
$5\%$, but a DerSHAP value of more than $5\%$. The level $5\%$ is important
since some researchers suggest freezing input parameters that have importance
less than $5\%$. At this threshold, the Shapley value approach would not
freeze $\delta_{2}$, but the activity score approach would.

(a) Activity scores

(b) DerSHAP values

Figure 2: Normalized activity scores and DerSHAP values for the HIV model

###  5.3  Amazon Stock Data

The data, which is from Kaggle222https://www.kaggle.com/code/gcmadhan/amazon-
stock-prediction-nn-fbprophet, is the Amazon stock price data from Jan. 4,
2010 to Dec 30, 2016, with 1233 observations in total. The training data
contains the opening price, closing price, low price, high price and volume of
the stock with 863 observations. The problem is to estimate the closing price
(the target), from the remaining data (the features), using a machine learning
algorithm. In the notation of the previous sections, the output of the model,
$f(\textbf{x})$, is the closing price, and x is the input features. In this
example we use a linear regression model to obtain $f(\textbf{x})$, and use
the Amazon stock data to train the model.

Table 2 displays the correlation between each input feature and closing price
for the whole data. The opening price, low price, and high price have a large
correlation with the stock’s closing price. The volume of the stock has a
smaller negative correlation.

Features | Correlation with target
---|---
opening (Feature 0) | 0.9996
low (Feature 1) | 0.9998
high (Feature 2) | 0.9998
volume (Feature 3) | -0.2386
Table 2: Correlation between input features and target (closing price)

In the numerical results, we assume the input features are dependent and
follow the normal distribution. We use the Python package linear
SHAP333https://shap.readthedocs.io/en/latest/generated/shap.explainers.Linear.html
by Lundberg, which claims to compute SHAP values accounting for correlations.
We compute the expectations $\mathbb{E}[\partial f_{i}^{2}]$ and
$\mathbb{E}[\partial f_{i}\partial f_{j}]$ in DerSHAP using Monte Carlo, based
on 1,000 random normal variates where the covariance matrix and the mean
vector are obtained from the data. We use forward differences with an
increment of $10^{-6}$ to estimate the partial derivatives. Fig. 3 plots the
DerSHAP and SHAP values for the linear regression model. Both methods agree on
the importance of inputs: the most important ones are high price, low price,
open price, and the least important one is volume.

(a) DerSHAP

(b) SHAP

Figure 3: DerSHAP and SHAP values for dependent Amazon stock data

Table 3 presents the computational time444Both methods are implemented using
Python, on a MacBook Pro M2 Chip. We use the code given by [10] for SHAP. to
obtain the results in Fig. 3. The computational time for both methods are
similar.

Method | Time(seconds)
---|---
DerSHAP | 0.18
SHAP | 0.22
Table 3: Computation time comparison - Amazon stock data

###  5.4 Boston Housing Data

The data, which is from Kaggle555https://www.kaggle.com/c/boston-housing,
consists of 489 data points where there are eight features: LSTAT (% lower
status of the population), INDUS (proportion of non-retail business acres per
town), NOX (nitric oxides concentration (parts per 10 million)), PTRATIO
(pupil-teacher ratio by town), RM (average number of rooms per dwelling, TAX
(full-value property-tax rate per $10,000), DIS (weighted distances to five
Boston employment centres), AGE (proportion of owner-occupied units built
prior to 1940). The target, $f(\textbf{x})$, is MEDV (median value of owner-
occupied homes in $1000s). Using the data we train two different models,
support vector regression (SVR) and gradient boosting regression (GBR). The
mean and standard deviation of the mean square error for the SVR model is
$-0.04(+/-0.03)$ and the GBR model is $-0.03(+/-0.02)$.

For DerSHAP, we first assume the input features are independent and follow
normal distribution, with parameters estimated from data. We use 10,000 random
samples to estimate the expectations, and use forward differences with an
increment of $10^{-6}$ to estimate the partial derivatives. For KernelSHAP, we
first set the training sample for the explainer to 50. Fig. 4 plots the
DerSHAP and KernelSHAP values for the SVR model. Both methods agree on the
most important inputs: RM, LSTAT, and the three least important inputs (less
than 5% in DerSHAP): NOX, AGE, PTRATIO. The third and fourth important inputs
for DerSHAP are TAX and INDUS, whereas for KernelSHAP, they are INDUS and DIS.
If we increase the training sample for KernelSHAP from 50 to 490, the third
and fourth important variables for KernelSHAP matches the results of DerSHAP.
The computing time for KernelSHAP for 490 samples is, however, 3.7 hours,
whereas DerSHAP gets the same information in 6.81 seconds.

(a) DerSHAP

(b) KernelSHAP

Figure 4: DerSHAP and SHAP values for Boston housing data - SVR Model

Fig. 5 plots the DerSHAP and KernelSHAP values for the GBR model. DerSHAP
identifies the most important inputs as: RM, PTRATIO and LSTAT. However,
KernelSHAP identifies the most important inputs as: LSTAT, RM, PTRATIO. DIS
has a DerSHAP value of nearly zero, but its KernelSHAP value is not
negligible. The conclusions of DerSHAP and KernelSHAP are clearly different,
and the differences remain even if we raise the size of the training sample
for KenelSHAP to 490 (unlike the SVR case, where KernelSHAP converged to the
outcomes of DerSHAP.)

(a) DerSHAP

(b) KernelSHAP

Figure 5: DerSHAP and SHAP values for Boston housing data - GBR Model

Table 4 presents the computational time of DerSHAP and KernelSHAP with the SVR
and GBR models. The computational time for KernelSHAP is about factors of 122
and 539 higher than DerSHAP. Recall that KernelSHAP uses a training sample of
size 50, and DerSHAP uses a Monte Carlo sample size of 10,000, for each model,
in these comparisons.

Model | Method | Time(seconds)
---|---|---
SVR | DerSHAP | 6.81
| KernelSHAP | 3676.19
GBR | DerSHAP | 9.05
| KernelSHAP | 1106.08
Table 4: Computation time comparison - Boston housing problem

An advantage of DerSHAP is it can be used with dependent data, unlike
KernelSHAP. Next we model the input features as correlated normal random
variables, which is the appropriate choice, using the data to estimate the
correlation matrix and the mean vector. Fig. 6 displays the DerSHAP values for
the dependent inputs, for the SVR and GBR models. The importance ranking of
features under DerSHAP is similar between the independent and dependent input
assumptions, with one notable exception. In the case of independence
assumption with the SVR model, the DerSHAP value for DIS is 6.9%, however,
when we model features as dependent inputs, the value for DIS becomes 4.3% - a
negligible value for input freezing.

(a) DerSHAP (dependent inputs - SVR)

(b) DerSHAP (dependent inputs - GBR)

Figure 6: DerSHAP values for dependent inputs for Boston housing data for SVR
and GBR models.

##  6 Conclusions

DerSHAP is computationally efficient: it has linear complexity as opposed to
exponential complexity of other general purpose Shapley algorithms. It can be
used with independent and dependent inputs. Unlike DGSM and activity score,
DerSHAP incorporates 2nd order interactions via mixed partial derivatives. The
efficiency, flexibility, and the ability to include higher order interactions
make DerSHAP a competitive method for global sensitivity analysis and machine
learning explainability.

## References

  * [1] F. Campolongo, J. Cariboni, and A. Saltelli, An effective screening design for sensitivity analysis of large models, Environmental Modelling & Software, 22 (2007), pp. 1509–1518.
  * [2] P. G. Constantine and P. Diaz, Global sensitivity metrics from active subspaces, Reliability Engineering & System Safety, 162 (2017), pp. 1–13.
  * [3] P. G. Constantine, E. Dow, and Q. Wang, Active subspace methods in theory and practice: applications to kriging surfaces, SIAM Journal on Scientific Computing, 36 (2014), pp. A1500–A1524.
  * [4] P. Diaz, P. Constantine, K. Kalmbach, E. Jones, and S. Pankavich, A modified SEIR model for the spread of Ebola in Western Africa and metrics for resource allocation, Applied Mathematics and Computation, 324 (2018), pp. 141–155.
  * [5] S. Kucherenko and B. Iooss, Derivative based global sensitivity measures, arXiv preprint arXiv:1412.2619, (2014).
  * [6] S. Kucherenko, S. Tarantola, and P. Annoni, Estimation of global sensitivity indices for models with dependent variables, Computer Physics Communications, 183 (2012), pp. 937–946.
  * [7] I. E. Kumar, S. Venkatasubramanian, C. Scheidegger, and S. Friedler, Problems with Shapley-value-based explanations as feature importance measures, in International Conference on Machine Learning, PMLR, 2020, pp. 5491–5500.
  * [8] T. Loudon and S. Pankavich, Mathematical analysis and dynamic active subspaces for a long term model of HIV, Mathematical Biosciences & Engineering, 14 (2016), pp. 709–733.
  * [9] S. M. Lundberg, G. Erion, H. Chen, A. DeGrave, J. M. Prutkin, B. Nair, R. Katz, J. Himmelfarb, N. Bansal, and S.-I. Lee, Consistent individualized feature attribution for tree ensembles, arXiv preprint, (2018).
  * [10] S. M. Lundberg and S.-I. Lee, A unified approach to interpreting model predictions, Advances in neural information processing systems, 30 (2017).
  * [11] M. Mase, A. B. Owen, and B. Seiler, Explaining black box decisions by Shapley cohort refinement, arXiv preprint arXiv:1911.00467, (2019).
  * [12] A. B. Owen, Sobol’ indices and Shapley value, SIAM/ASA Journal on Uncertainty Quantification, 2 (2014).
  * [13] A. B. Owen and C. Prieur, On Shapley value for measuring importance of dependent inputs, SIAM/ASA Journal on Uncertainty Quantification, 5 (2017).
  * [14] B. B. Seiler, M. Mase, and A. B. Owen, What makes you unique?, Electronic Journal of Statistics, 17 (2023), pp. 1–18.
  * [15] L. S. Shapley, A value for n-person games, in Contribution to the Theory of Games II (Annals of Mathematics Studies 28), H. W. Kuhn and A. W. Tucker, eds., vol. 69, Princeton University Press, Princeton, NJ, 1953, pp. 307–317.
  * [16] I. M. Sobol’, On derivative-based global sensitivity criteria, Mathematical Models and Computer Simulations, 3 (2011), pp. 419–423.
  * [17] I. M. Sobol’ and S. Kucherenko, Derivative based global sensitivity measures and their link with global sensitivity indices, Mathematics and Computers in Simulation, 79 (2009), pp. 3009–3017.
  * [18] E. Song, B. L. Nelson, and J. Staum, Shapley effects for global sensitivity analysis: Theory and computation, SIAM/ASA Journal on Uncertainty Quantification, 4 (2016), pp. 1060–1083.
  * [19] G. Van den Broeck, A. Lykov, M. Schleich, and D. Suciu, On the tractability of SHAP explanations, Journal of Artificial Intelligence Research, 74 (2022), pp. 851–886.