# Feature-Specific Coefficients of Determination in Tree Ensembles

Zhongli Jiang  Dabao Zhang  Min Zhang Zhongli Jiang is with Purdue University,
West Lafayette, IN 47907 USA (email: jiang548@purdue.edu).Dabao Zhang and Min
Zhang are with University of California, Irvine, CA 92697 USA (email:
<EMAIL_ADDRESS>minzhang@hs.uci.edu).

###### Abstract

Tree ensemble methods provide promising predictions with models difficult to
interpret. Recent introduction of Shapley values for individualized feature
contributions, accompanied with several fast computing algorithms for
predicted values, shows intriguing results. However, individualizing
coefficients of determination, aka $R^{2}$, for each feature is challenged by
the underlying quadratic losses, although these coefficients allow us to
comparatively assess single feature’s contribution to tree ensembles. Here we
propose an efficient algorithm, Q-SHAP, that reduces the computational
complexity to polynomial time when calculating Shapley values related to
quadratic losses. Our extensive simulation studies demonstrate that this
approach not only enhances computational efficiency but also improves
estimation accuracy of feature-specific coefficients of determination.

Our proposed Q-SHAP algorithm significantly advances the interpretability of
tree ensemble methods by reducing the computational complexity of calculating
Shapley values for quadratic losses, such as $R^{2}$, to polynomial time. By
enabling an efficient and accurate decomposition of $R^{2}$, Q-SHAP fosters
greater understanding and transparency in machine learning applications,
particularly in biomedical fields where model interpretability is paramount.
Furthermore, Q-SHAP’s framework has the potential to extend to general loss
functions via approximation, opening new avenues for research and broadening
its applicability across various machine learning models. This work
significantly contributes to making complex models more understandable and
trustworthy, enhancing their adoption and impact across diverse domains.

Explainable artificial intelligence (XAI), Shapley values, coeffcients of
determination, tree ensembles

##  1 Introduction

Models built with tree ensembles are powerful but often complicated, making it
challenging to understand the influence of inputs. Feature importance plays a
critical role in demystifying these models and enhancing their
interpretability by assigning each input feature a score. This is crucial in
domains like healthcare and biomedicine, where trust and interpretation of the
model are essential [1, 2]. Common feature importance measures like gain can
be inconsistent [3] while permutation importance lacks theoretical foundations
[4].

Shapley values, derived from cooperative game theory and introduced by Shapley
[5], offer a robust method for the fair distribution of payoffs generated by a
coalition of players. This can be analogously applied to assess the
contribution of each feature in a machine learning model. It ensures that each
feature’s contribution is assessed by considering all possible combinations of
features, thereby providing a comprehensive understanding of feature impacts.
Recent applications of Shapley values have focused on local interpretation [6,
3, 7], where they are employed to examine the influence of individual features
on specific predictions. Nonetheless, there are numerous scenarios where
global importance is preferred, such as analysis of the role of a feature
across the entire dataset [8, 9].

Among the works that compute Shapley values in a global context, a popular
approach is to use model variance decomposition. Lipovetsky and Conklin [10]
decomposed $R^{2}$ in linear regression, offering consistent interpretations
even in the presence of multicollinearity. Owen and Prieur [11] also conducted
a conceptual analysis of Shapley values for the variance. However, computation
remains a significant challenge, as the calculation of Shapley values grows
exponentially with the number of features. To address this issue, several
Monte Carlo-based methods have been proposed to effectively reduce the
computational burden [12, 9, 13].

Although Monte Carlo-based, model-agnostic methods are more efficient than
brute-force approaches, they are still computationally intensive, especially
when dealing with high-dimensional data that requires extensive feature
permutation sampling to ensure consistency [6, 14]. This challenge has
prompted the development of methods that leverage the specific structures of
tree-based models. However, much of the focus has been on explaining
individual predictions, as seen with TreeSHAP [3], FastTreeSHAP [15], and
LinearTreeSHAP [16]. Bénard et al. [17] considered population-level importance
using $R^{2}$, specifically tailored for random forests [18].

To the best of our knowledge, there is no available method to calculate
Shapley values of quadratic losses by leveraging structures of decision trees
for fast computation. In this paper, we propose Q-SHAP, which can decompose
quadratic terms of predicted values of a decision tree into each feature’s
attribute in polynomial time. It leads to fast computation of feature-specific
$R^{2}$ for a decision tree. We also extend our approach to Gradient Boosted
Decision Trees.

The rest of the paper is structured as follows. In Section 2, we provide a
brief overview of Shapley values of $R^{2}$. In Section 3, we present our
proposed algorithm Q-SHAP to calculate Shapley values of $R^{2}$ in polynomial
time for single trees, and then extend the approach for tree ensembles in
Section 4. We justify the efficacy and efficiency of the algorithm using
extensive simulations in Section 5 and real data analysis in high dimension in
Section 6. We conclude the paper with a discussion in Section 7.

##  2 Shapley Values of $R^{2}$ for Individual Features

###  2.1 Model specification

Here we investigate a specific label $Y$ and its explainability by a full set
of $p$ features $X=(X_{1},X_{2},\cdots,X_{p})$. For any subset
$F\subseteq\mathcal{P}=\\{1,2,\cdots,p\\}$, we define the corresponding set of
features as $X_{F}=(X_{j})_{j\in F}$.

Suppose that, for any set of feature $X_{F}$, an oracle model $m_{F}$ can be
built such that, for any specific value $x=(x_{1},x_{2},\cdots,x_{p})$,

$m_{F}(x)=E[Y|X_{F}=(x_{j})_{j\in F}].$

The Shapley value of $j$-th feature, in terms of its contribution to the total
variation, is defined as

$\displaystyle\phi_{\rho^{2},j}$ $\displaystyle=$
$\displaystyle\frac{1}{p~{}var(m_{\emptyset})}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}$ (1)
$\displaystyle\times\left(var(m_{F\cup\\{j\\}})-var(m_{F})\right),$

where $|F|$ is the number of features in $F$. The term $var(m_{F\cup\\{j\\}})$
is the variance explained by feature set $F\cup\\{j\\}$ and the term
$var(m_{F})$ is the variance explained solely by set $F$. This definition is
analogous to Covert et al. [9] and Williamson and Feng [13]. By averaging over
all possible feature combinations, the Shapley values are the only solution
that satisfies the desired properties of symmetry, efficiency, additivity, and
dummy [5].

Note that $\sqrt{\phi_{\rho^{2},j}}$ can be used to describe the general
association of feature $j$ with the label $Y$. In fact, when a linear model is
under consideration with independent features $X_{1},X_{2},\cdots,X_{p}$, we
have

$\sqrt{\phi_{\rho^{2},j}}=|corr(Y,X_{j})|,$

i.e., $\sqrt{\phi_{\rho^{2},j}}$ measures the scale of Pearson’s correlation
between $Y$ and $X_{j}$. Therefore, $\sqrt{\phi_{\rho^{2},j}}$ extends the
concept of correlation from a linear model to a general nonlinear model.

###  2.2 Empirical Estimation

Suppose we have a set of data with sample size $n$ observed for both label and
features as

$\displaystyle\mathbf{Y}=(y_{1},y_{2},\cdots,y_{n});$
$\displaystyle\mathbf{X}=\left(\mathbf{X}_{\cdot 1},\mathbf{X}_{\cdot
2},\cdots,\mathbf{X}_{\cdot
p}\right)=\left(\mathbf{x}_{1\cdot}^{T},\mathbf{x}_{2\cdot}^{T},\cdots,\mathbf{x}_{n\cdot}^{T}\right)^{T}.$

Accordingly, we denote the observed data of features in subset $F$ as

$\mathbf{X}_{\cdot F}=(\mathbf{X}_{\cdot j})_{j\in F}.$

Suppose that, for each subset $F$ of features, a single optimal model
$\hat{m}_{F}$ is built on data $(\mathbf{Y},\mathbf{X}_{\cdot F})$. Then the
$i$-th label can be predicted with

$\hat{y}_{i}(\mathbf{X}_{\cdot F})=\hat{m}_{F}(\mathbf{x}_{i\cdot}).$

###  2.3 From $R^{2}$ to a Quadratic Loss

We will establish the connection of $R^{2}$ to a quadratic loss through
equation (1). We define the quadratic loss on the optimal model $\hat{m}_{F}$
as

$\displaystyle
Q_{F}=\sum_{i=1}^{n}\left(y_{i}-\hat{m}_{F}(\mathbf{x}_{i\cdot})\right)^{2}$
(2)

for any set of features $F$. With
$m_{\emptyset}(\mathbf{x}_{i\cdot})=\bar{y}$, we have

$Q_{\emptyset}=\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}.$

Following the law of total variance, we can estimate $var(m_{F})$ by

$\widehat{var}(m_{F})=\left(Q_{\emptyset}-Q_{F}\right)\big{/}n.$

Thus, an empirical estimate of (1) is

$\displaystyle\phi_{R^{2},j}$ $\displaystyle=$
$\displaystyle\frac{1}{pQ_{\emptyset}}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}$
$\displaystyle\times\left((Q_{\emptyset}-Q_{F\cup\\{j\\}})-(Q_{\emptyset}-Q_{F})\right)$
$\displaystyle=$
$\displaystyle-\frac{1}{pQ_{\emptyset}}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}(Q_{F\cup\\{j\\}}-Q_{F}),$

which is proportional to a Shapley value for the sum of squared errors, i.e.,
the quadratic loss in (2).

###  2.4 From Quadratic Loss to Q-SHAP

We now further reduce Shapley values of the sum of squared errors to Shapley
values of linear and quadratic terms of predicted values. Expanding the loss
function in (2), we can rewrite,

$\displaystyle\phi_{R^{2},j}$ $\displaystyle=$
$\displaystyle-\frac{1}{pQ_{\emptyset}}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}\sum_{i=1}^{n}\left(\hat{m}_{F\cup
j}^{2}(\mathbf{x}_{i\cdot})\right.$
$\displaystyle\left.-\hat{m}_{F}^{2}(\mathbf{x}_{i\cdot})-2(\hat{m}_{F\cup
j}(\mathbf{x}_{i\cdot})-\hat{m}_{F}(\mathbf{x}_{i\cdot}))y_{i}\right).$

To calculate this, we define the Shapley value for each sample $i$ as,

$\displaystyle\phi_{R^{2},j}(\mathbf{x}_{i\cdot})$ $\displaystyle=$
$\displaystyle-\frac{1}{pQ_{\emptyset}}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}\left(\hat{m}_{F\cup j}^{2}(\mathbf{x}_{i\cdot})\right.$
$\displaystyle\left.-\hat{m}_{F}^{2}(\mathbf{x}_{i\cdot})-2(\hat{m}_{F\cup
j}(\mathbf{x}_{i\cdot})-\hat{m}_{F}(\mathbf{x}_{i\cdot}))y_{i}\right)$
$\displaystyle=$
$\displaystyle-\frac{1}{pQ_{\emptyset}}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}\left(\hat{m}_{F\cup
j}^{2}(\mathbf{x}_{i\cdot})-\hat{m}_{F}^{2}(\mathbf{x}_{i\cdot})\right)$
$\displaystyle+\frac{2y_{i}}{pQ_{\emptyset}}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}\left(\hat{m}_{F\cup
j}(\mathbf{x}_{i\cdot})-\hat{m}_{F}(\mathbf{x}_{i\cdot})\right),$

which is a linear combination of two sets of Shapley values, i.e., Shapley
values of predicted value $\hat{m}_{F}$, which are ready to be calculated [3,
15, 16], and Shapley values of the quadratic term of predicted value
$\hat{m}_{F}^{2}$, i.e.,

$\displaystyle\phi_{\hat{m}^{2},j}(\mathbf{x}_{i\cdot})$ $\displaystyle=$
$\displaystyle\frac{1}{p}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}$ (3) $\displaystyle\times\left(\hat{m}_{F\cup
j}^{2}(\mathbf{x}_{i\cdot})-\hat{m}_{F}^{2}(\mathbf{x}_{i\cdot})\right),$

for which we will develop the algorithm Q-SHAP to calculate. For the rest of
the paper, we will focus on computing the Shapley values in Equation (3) in
polynomial time for tree-based models and carrying it over to calculate
feature-specific $R^{2}$.

##  3 The Algorithm Q-SHAP for Single Trees

Unlike most regression problems that can yield infinitely many predictions
across the diverse input space $\mathbf{X}$, decision trees restrict
predictions to a finite set, specifically to the values at each leaf node.
This nature revitalizes hope in previously unattainable solutions to Shapley
values in decision tree-based models. The main idea of our algorithm lies on
the fact that Shapley values with various targets such as predictions and
various loss functions, are essentially weighted functions of the leaf
predictions. While Lundberg et al. [14] suggests that explaining the loss
function for a “path-dependent” algorithm is challenging, we provide an exact
solution to decomposing quadratic losses using our Q-SHAP. Q-SHAP works for
tree ensembles. For simplicity, we here illustrate it with a single decision
tree.

###  3.1 Notations

We assume the underlying decision tree has the maximum depth at $D$ and a
total of $L$ leaves, and use $l$ to denote a specific leaf. We further
introduce a dot product for polynomials for subsequent calculation. For two
polynomials

$\displaystyle A(z)=\sum_{i=0}^{n}a_{i}z^{i},\ \
B(z)=\sum_{i=0}^{n}b_{i}z^{i},$

we define their dot product as

$A(z)\cdot B(z)=\sum_{i=0}^{n}a_{i}b_{i}.$

###  3.2 Distributing the prediction to leaves

Decision trees match each data point to one leaf for prediction. However, for
our prediction defined on any subset $F$, a data point $\mathbf{x}_{i\cdot}$
can fall into multiple leaves due to the uncertainty by unspecified features
$\mathcal{P}\backslash F$. We can calculate
$\hat{m}_{F}(\mathbf{x}_{i\cdot})$, following TreeSHAP, as the empirical mean
by aggregating the weighted prediction on each leaf,

$\hat{m}_{F}(\mathbf{x}_{i\cdot})=\sum_{l}\hat{m}_{F}^{l}(\mathbf{x}_{i\cdot}),$
(4)

where $\hat{m}_{F}^{l}$ is the weighted prediction from leaf $l$ in a tree
built on feature set $F$.

Given an oracle tree built on all available features, we try to recover the
oracle tree for a subset of features without rebuilding, following Bifet et
al. [16] and Karczmarz et al. [19]. For example, let us take the tree built on
two features $X_{1}$ and $X_{2}$ as shown in Figure 2. Figure 2 can be viewed
as the oracle tree built solely on feature $X_{1}$, and hence $X_{2}$ no
longer exists in the tree. Therefore, we replace it with a pseudo internal
node to preserve the structure of the original full oracle tree and pave the
way for further formulation.

Figure 1: Decision tree built on $X_{1}$ and $X_{2}$

Figure 2: Hypothetic decision tree with $X_{1}$ only

With a data point $\mathbf{x}_{i\cdot}=(40,25)$, we illustrate the calculation
in Equation (4) by first calculating the predicted value for the tree in
Figure 2,

$\displaystyle\hat{m}_{\\{1\\}}(\mathbf{x}_{i\cdot})$ $\displaystyle=$
$\displaystyle\hat{m}_{\\{1\\}}^{l_{1}}(\mathbf{x}_{i\cdot})+\hat{m}_{\\{1\\}}^{l_{2}}(\mathbf{x}_{i\cdot})+\hat{m}_{\\{1\\}}^{l_{3}}(\mathbf{x}_{i\cdot})$
$\displaystyle=$ $\displaystyle 0\times 10+0.5\times 20+0.5\times 30.$

When the tree is built with an additional feature $j=2$ as shown in Figure 2,
we have the predicted value

$\displaystyle\hat{m}_{\\{1,2\\}}(\mathbf{x}_{i\cdot})$ $\displaystyle=$
$\displaystyle\hat{m}_{\\{1,2\\}}^{l_{1}}(\mathbf{x}_{i\cdot})+\hat{m}_{\\{1,2\\}}^{l_{2}}(\mathbf{x}_{i\cdot})+\hat{m}_{\\{1,2\\}}^{l_{3}}(\mathbf{x}_{i\cdot})$
$\displaystyle=$ $\displaystyle\textbf{1}\times 0\times
10+\mathbf{0.5^{-1}}\times 0.5\times 20+\textbf{0}\times 0.5\times 30,$

where the bold numbers reweight $\hat{m}^{l}_{\\{1\\}}(\mathbf{x}_{i\cdot})$
for $\hat{m}^{l}_{\\{1,2\\}}(\mathbf{x}_{i\cdot})$.

Let us take a closer look at these weights, each corresponding to one leaf.
For leaf $l_{1}$, the weight is 1 since the newly added feature $X_{2}$ is not
involved in its path and the reweighted prediction remains as zero. For leaf
$l_{2}$, the reweighted prediction is lifted up with the weight inversely
proportional to the previous probability because $\mathbf{x}_{i\cdot}$ follows
its path to the leaf with probability 1. On the other hand, although the path
to leaf $l_{3}$ includes the newly added feature but $\mathbf{x}_{i\cdot}$
doesn’t follow this path, resulting in a weight at 0. Next we will generalize
such a reweighting strategy to calculate (4) for trees with different sets of
features.

Denote $F^{l}$ the features involved in the path to leaf $l$ and
$F^{l}(\mathbf{x}_{i\cdot})$ the subset of $F^{l}$ whose decision criteria are
satisfied by $\mathbf{x}_{i\cdot}$. Note that each feature $j\in F^{l}$ may
appear multiple times in the path to leaf $l$ so we denote $n^{l}_{j,c}$ the
number of samples passing through the node which is attached to the $c$-th
appearance. We similarly define $n^{l}_{j,c}(\mathbf{x}_{i\cdot})$ for each
feature $j\in F^{l}(\mathbf{x}_{i\cdot})$.

For any feature $j\in\mathcal{P}$, we can define the weight function based on
a partition of $\mathcal{P}$ into three subsets $F^{l}(\mathbf{x}_{i\cdot})$,
$F^{l}\backslash F^{l}(\mathbf{x}_{i\cdot})$, and $\mathcal{P}\backslash
F^{l}$,

$w_{j}^{l}(\mathbf{x}_{i\cdot})\triangleq\begin{cases}\prod_{c}{\frac{n^{l}_{j,c}}{n^{l}_{j,c}(\mathbf{x}_{i\cdot})}},&\text{if
}j\in F^{l}(\mathbf{x}_{i\cdot});\\\ 0,&\text{if }j\in F^{l}\backslash
F^{l}(\mathbf{x}_{i\cdot});\\\ 1,&\text{if }j\in\mathcal{P}\backslash
F^{l}.\end{cases}$

Therefore, for $j\notin F$, we have

$\hat{m}_{F\cup
j}^{l}(\mathbf{x}_{i\cdot})=w_{j}^{l}(\mathbf{x}_{i\cdot})\hat{m}_{F}^{l}(\mathbf{x}_{i\cdot}).$

Recursive application of the above formula leads to

$\hat{m}_{F}^{l}(\mathbf{x}_{i\cdot})=\prod_{k\in
F}w^{l}_{k}(\mathbf{x}_{i\cdot})\hat{m}_{\emptyset}^{l},$

where $\hat{m}_{\emptyset}^{l}=\hat{m}^{l}\frac{n^{l}}{n}$ with $n^{l}$ the
sample size at leaf $l$, $n$ the total sample size, and $\hat{m}^{l}$ the
predicted value at leaf $l$ based on the model built on all features.

When $F=\emptyset$, the above result reduces to

$\hat{m}_{\emptyset}(\mathbf{x}_{i\cdot})=\sum_{l}\hat{m}^{l}\frac{n^{l}}{n},$

so the optimal prediction is just the mean for all data points, which is
consistent with $\hat{m}_{\emptyset}(\mathbf{x}_{i\cdot})=\bar{y}$.

We can rewrite (3) as

$\displaystyle\phi_{\hat{m}^{2},j}(\mathbf{x}_{i\cdot})$ $\displaystyle=$
$\displaystyle\frac{1}{p}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}$
$\displaystyle\times\left(\left[\sum_{l}w_{j}^{l}(\mathbf{x}_{i\cdot})\prod_{k\in
F}w_{k}^{l}(\mathbf{x}_{i\cdot})\hat{m}^{l}_{\emptyset}\right]^{2}\right.$
$\displaystyle\left.-\left[\sum_{l}\prod_{k\in
F}w_{k}^{l}(\mathbf{x}_{i\cdot})\hat{m}^{l}_{\emptyset}\right]^{2}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{p}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}$
$\displaystyle\times\left(\sum_{l}\left(w_{j}^{l2}(\mathbf{x}_{i\cdot})-1\right)\hat{m}^{l2}_{\emptyset}\prod_{k\in
F}w_{k}^{l2}(\mathbf{x}_{i\cdot})\right)$
$\displaystyle+\frac{2}{p}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}$ $\displaystyle\times\left(\sum_{l_{1}\neq
l_{2}}(w_{j}^{l_{1}}(\mathbf{x}_{i\cdot})w_{j}^{l_{2}}\left(\mathbf{x}_{i\cdot})-1\right)\hat{m}_{\emptyset}^{l_{1}}\hat{m}^{l_{2}}_{\emptyset}\right.$
$\displaystyle\times\left.\prod_{k\in
F}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot})\right)$
$\displaystyle\triangleq$ $\displaystyle
T_{1,j}(\mathbf{x}_{i\cdot})+2T_{2,j}(\mathbf{x}_{i\cdot}).$

We further define, for leaves $l_{1}$ and $l_{2}$,

$\displaystyle T_{j}^{l_{1}l_{2}}(\mathbf{x}_{i\cdot})$ (5) $\displaystyle=$
$\displaystyle\frac{1}{p}\sum_{F\subseteq\mathcal{P}\backslash\\{j\\}}\begin{pmatrix}p-1\\\
|F|\end{pmatrix}^{-1}\bigg{(}(w_{j}^{l_{1}}(\mathbf{x}_{i\cdot})w_{j}^{l_{2}}(\mathbf{x}_{i\cdot})-1)$
$\displaystyle\times\hat{m}_{\emptyset}^{l_{1}}\hat{m}^{l_{2}}_{\emptyset}\prod_{k\in
F}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot})\bigg{)},$

and we have

$T_{2,j}(\mathbf{x}_{i\cdot})=\sum_{l_{1}\neq
l_{2}}T_{j}^{l_{1}l_{2}}(\mathbf{x}_{i\cdot}),\ \ \
T_{1,j}(\mathbf{x}_{i\cdot})=\sum_{l}T_{j}^{ll}(\mathbf{x}_{i\cdot}).$

Therefore, we will focus on the calculation of
$T_{j}^{l_{1}l_{2}}(\mathbf{x}_{i\cdot})$ in (5) the rest of this section.

We can reduce the calculation of $\prod_{k\in
F}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot})$ in (5)
by only calculating $\prod_{k\in
F_{-}}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot})$
with

$F_{-}=F\cap(F^{l_{1}}\cup F^{l_{2}}),$

because

$\prod_{k\in F\backslash
F_{-}}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot})=1.$

In combination with the proposition below, computation in (5) can be
dramatically reduced from the full feature set $\mathcal{P}$ to a set only
related to the corresponding leaves in a tree.

###### Proposition 1

For any well-defined $p,n,|F|$,

$\sum_{k=0}^{p-n}\frac{\binom{p-n}{k}}{p\binom{p-1}{|F|+k}}=\frac{1}{n\binom{n-1}{|F|}}.$

We leave the proof of Proposition 1 in Appendix .1. Further denote
$n_{12}=|F^{l_{1}}\cup F^{l_{2}}|$ and a polynomial of $z$,

$P^{l_{1}l_{2}}(z)=\prod_{k\in F^{l_{1}}\cup F^{l_{2}}\backslash
j}(z+w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot})).$

We then define a coefficient polynomial

$C_{n_{12}}(z)=\frac{1}{\binom{n_{12}-1}{0}}z^{0}+\frac{1}{\binom{n_{12}-1}{1}}z^{1}+\ldots+\frac{1}{\binom{n_{12}-1}{n_{12}-1}}z^{n_{12}-1}.$

###### Theorem 1

$\displaystyle T_{j}^{l_{1}l_{2}}(\mathbf{x}_{i\cdot})$ $\displaystyle=$
$\displaystyle\frac{1}{n_{12}}(w_{j}^{l_{1}}(\mathbf{x}_{i\cdot})w_{j}^{l_{2}}(\mathbf{x}_{i\cdot})-1)\hat{m}_{\emptyset}^{l_{1}}\hat{m}_{\emptyset}^{l_{2}}$
$\displaystyle\times[C_{n_{12}}(z)\cdot P^{l_{1}l_{2}}(z)].$

Proof. With Proposition 1, we can write (5) as

$\displaystyle T_{j}^{l_{1}l_{2}}(\mathbf{x}_{i\cdot})$ $\displaystyle=$
$\displaystyle\frac{1}{n_{12}}(w_{j}^{l_{1}}(\mathbf{x}_{i\cdot})w_{j}^{l_{2}}(\mathbf{x}_{i\cdot})-1)\hat{m}_{\emptyset}^{l_{1}}\hat{m}_{\emptyset}^{l_{2}}$
$\displaystyle\times\sum_{t=0}^{n_{12}-1}\frac{1}{\binom{n_{12}-1}{t}}\sum^{|F|=t}_{F\subseteq
F^{l_{1}}\cup F^{l_{2}}\backslash j}\prod_{k\in
F}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot}).$

We further notice that $\sum^{|F|=t}_{F\subseteq F^{l_{1}}\cup
F^{l_{2}}\backslash j}\prod_{k\in
F}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot})$ is the
coefficient of $z^{t}$ in polynomial $P^{l_{1}l_{2}}(z)$, hence the equation
holds with $C_{n_{12}}(z)$ adjusting the weight based on the size of set $F$.
$\blacksquare$

We only need to consider feature $j\in$$|F^{l_{1}}\cup F^{l_{2}}|$ as,
otherwise, we have $T_{j}^{l_{1}l_{2}}(\mathbf{x}_{i\cdot})=0$ following the
definition of $w_{j}^{l}(\mathbf{x}_{i\cdot})$. Note that, when there is a
feature in set $F$ that doesn’t belong to $F^{l_{1}}(\mathbf{x}_{i\cdot})\cap
F^{l_{2}}(\mathbf{x}_{i\cdot})\backslash j$, we have

$\prod_{k\in
F}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot})=0.$

Thus we can further simplify the term to

$\displaystyle T_{j}^{l_{1}l_{2}}(\mathbf{x}_{i\cdot})$ $\displaystyle=$
$\displaystyle\frac{1}{n_{12}}(w_{j}^{l_{1}}(\mathbf{x}_{i\cdot})w_{j}^{l_{2}}(\mathbf{x}_{i\cdot})-1)\hat{m}_{\emptyset}^{l_{1}}\hat{m}_{\emptyset}^{l_{2}}\sum_{t=0}^{n_{12}-1}\frac{1}{\binom{n_{12}-1}{t}}$
$\displaystyle\times\sum^{|F|=t}_{F\subseteq
F^{l_{1}}(\mathbf{x}_{i\cdot})\cap F^{l_{2}}(\mathbf{x}_{i\cdot})\backslash
j}\prod_{k\in
F}w_{k}^{l_{1}}(\mathbf{x}_{i\cdot})w_{k}^{l_{2}}(\mathbf{x}_{i\cdot}).$

Consequently, the evaluation of $P^{l_{1}l_{2}}(z)$ can be reduced to a much
smaller set.

###  3.3 The Algorithm

In this section, we will introduce our algorithms. Theorem 1 demonstrates that
we can construct a polynomial form of the NP-problem. Now we introduce a fast
and stable evaluation for the dot product of a coefficient polynomial $C(z)$
where we know the coefficients and a polynomial $P(z)$ with a known product
form, involved in Theorem 1.

###### Proposition 2

Let $\omega$ be a vector of the complex $n$-th roots of unity whose element is
$\exp(\frac{2k\pi i}{n})$ for $k=0,1.,\ldots,n-1$, c the coefficient vector of
$C(z)$, and IFFT the Inverse Fast Fourier Transformation. Then

$C(z)\cdot P(z)=P(\omega)^{T}\text{IFFT}(c).$

The proof of Proposition 2 is shown in Appendix .1. We facilitate the
computation via the complex roots of unity because of their numerical
stability and fast operations in matrix multiplications. Due to the potential
issue of ill condition, especially at large degrees, our calculation avoids
inversion of the Vandermonde matrices, although it has been proposed to
facilitate the computing by Bifet et al. [16]. In addition, for each sample
size $n$, we only need calculate IFFT$(c)$ once, up to order $D$ in
$O(n\log(n))$ operations, and the results can be saved for the rest of
calculation through Q-SHAP. Note that term $k$ and term $n-k$ in $P(w)$ are
complex conjugates, and, for a real vector $c$, IFFT$(c)$ also has the
conjugate property for paired term $k$ and term $n-k$. Consequently, the dot
product of $P(\omega)$ and IFFT$(c)$ inherits the conjugate property and its
imaginary parts are canceled upon addition. Therefore, we only need evaluate
the dot product at half of the $n$ complex roots.

We can aggregate the values of leaf combinations to derive the Shapley values
of squared predictions using Q-SHAP as in Algorithm 1 and then calculate the
Shapley values of $R^{2}$ using RSQ-SHAP as in Algorithm 2. The calculation of
feature-specific $R^{2}$ uses the iterative Algorithm 1 instead of a recursive
one. The time complexity of the algorithm is $O(L^{2}D^{2})$ for a single
tree, which is super fast when the maximum tree depth is not too large.

Algorithm 1 Q-SHAP

Q-SHAP($\mathbf{x}_{i\cdot}$)

Initialize $T[j]=0$ for $j=1,\cdots,p$

for $l_{1}$ $\in$ index set ${0,\ldots},L-1$ do

for $l_{2}$ $\in$ index set $l_{1},\ldots,L-1$ do

Let $n_{12}=|F^{l_{1}}\cup F^{l_{2}}|$

for $j\in F^{l_{1}}\cup F^{l_{2}}$ do

Let
$t[j]=\frac{1}{n_{12}}[w_{j}^{l_{1}}(\mathbf{x}_{i\cdot})w_{j}^{l_{2}}(\mathbf{x}_{i\cdot})-1]\hat{m}_{\emptyset}^{l_{1}}\hat{m}_{\emptyset}^{l_{2}}\times$

$\hskip 42.67912pt[C_{n_{12}}(z)\cdot P^{l_{1}l_{2}}(z)]$

if $l_{1}\neq l_{2}$ then

$T[j]=T[j]+2t[j]$

else

$T[j]=T[j]+t[j]$

end if

end for

end for

end for

return $T=(T[1],T[2],\cdots,T[p])$

Algorithm 2 RSQ-SHAP

RSQ-SHAP($j$) =
$-\frac{1}{Q_{\emptyset}}\Sigma_{i=1}^{n}\\{\textbf{Q-SHAP}(\mathbf{x}_{i\cdot})[j]$

$\hskip 110.96556pt-2y_{i}\textbf{SHAP}(\mathbf{x}_{i\cdot})[j]\\}$

##  4 The Algorithm Q-SHAP for Tree Ensembles from Boosting

Tree ensembles from Gradient Boosted Machines (GBM) [20] greatly improve
predictive performance by aggregating many weak learners [21, 22, 23]. Each
tree, say tree $k$, is constructed on the residuals from the previous tree,
i.e., tree $k-1$. We assume that there are a total of $K$ trees in the
ensemble and the quadratic loss by the first $k$ trees, with all features in
$\mathcal{P}$, is $Q_{\mathcal{P}}^{(k)}$. Denoting
$Q_{\mathcal{P}}^{(0)}=Q_{\emptyset}$, the $k$-th tree reduces the loss by

$\Delta Q_{\mathcal{P}}^{(k)}=Q_{\mathcal{P}}^{(k-1)}-Q_{\mathcal{P}}^{(k)},$

with the tree ensemble reducing the total loss by

$Q_{\emptyset}-Q_{\mathcal{P}}^{(K)}=\sum_{k=1}^{K}\Delta
Q_{\mathcal{P}}^{(k)}.$

Per our interest in feature-specific $R^{2}$, we resort to the quadratic loss
defined as the sum of squared errors in (2).

On the other hand, the $k$-th tree provides the prediction
$\hat{m}_{\mathcal{P}}^{(k)}(\mathbf{x}_{i\cdot})$. Therefore, the prediction
by the first $k$ trees can be recursively calculated as

$\hat{y_{i}}^{(k)}(\mathbf{X})=\hat{y}_{i}^{(k-1)}(\mathbf{X})+\alpha\hat{m}_{\mathcal{P}}^{(k)}(\mathbf{x}_{i\cdot}),$

where $\alpha$ is the learning rate and
$\hat{y}_{i}^{(0)}(\mathbf{X})\equiv\bar{y}$. Note that the residuals after
building $(k-1)$ tree are

$\\{r_{i}^{(k-1)}=y_{i}-\hat{y}^{(k-1)}_{i}(\mathbf{X}):i=1,2,\cdots,n\\},$

which are taken to build the $k$-th tree. Thus,

$\displaystyle\Delta Q_{\mathcal{P}}^{(k)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}(r_{i}^{(k-1)})^{2}-\sum_{i=1}^{n}(r_{i}^{(k-1)}-\alpha\hat{m}^{(k)}_{\mathcal{P}}(\mathbf{x}_{i\cdot}))^{2}$
$\displaystyle=$
$\displaystyle-\sum_{i=1}^{n}(\alpha^{2}\hat{m}^{(k)2}_{\mathcal{P}}(\mathbf{x}_{i\cdot})-2\alpha
r_{i}^{(k-1)}\hat{m}^{(k)}_{\mathcal{P}}(\mathbf{x}_{i\cdot})).$

Therefore, the decomposition of $\Delta Q_{\mathcal{P}}^{(k)}$ for feature-
specific Shapley values in the $k$-th tree is again subject to the
decomposition of two sets of values, i.e., both the predicted value
$\hat{m}_{\mathcal{P}}^{(k)}(\mathbf{x}_{i\cdot})$ and its quadartic term,
which can be similarly dealt with the algorithms proposed in the previous
section.

##  5 Simulation Study

One of the challenges in assessing methods that explain predictions is the
typical absence of a definitive ground truth. Therefore, to fairly demonstrate
the fidelity of our methodology, we must rely on synthetic data that allows
for the calculation of the theoretical Shapley values. Here we consider three
different models,

$\displaystyle a.$ $\displaystyle Y=4X_{1}-5X_{2}+6X_{3}+\epsilon;$
$\displaystyle b.$ $\displaystyle
Y=4X_{1}-5X_{2}+6X_{3}+3X_{1}X_{2}-X_{1}X_{3}+\epsilon;$ $\displaystyle c.$
$\displaystyle Y=4X_{1}-5X_{2}+6X_{3}+3X_{1}X_{2}-X_{1}X_{2}X_{3}+\epsilon.\ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $

All three features involved in the models are generated from Bernoulli
distributions,

$\displaystyle X_{1}\sim Bernoulli(0.6),$ $\displaystyle X_{2}\sim
Bernoulli(0.7),$ $\displaystyle X_{3}\sim Bernoulli(0.5).$

We also simulate additional nuisance features independently from
$Bernoulli(0.5)$ to make the total number of features $p=100$ and $p=500$,
respectively. The error term $\epsilon$ is generated from
$N(0,\sigma_{\epsilon}^{2})$ with $\sigma_{\epsilon}$ at three different
levels, i.e., $0.5$, $1$, and $1.5$. The theoretical values of total $R^{2}$
and feature-specific $R^{2}$ are shown in Table 1 of Appendix .2 as well as
indicated by dashed lines in Figure 3.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 3: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=1,000$, $p=100$, and $\sigma_{\epsilon}=1.5$. The dashed lines
show the theoretical $R^{2}$.

We evaluate the performance of three different methods, our proposed Q-SHAP,
SAGE by Covert et al. [9], and SPVIM by Williamson and Feng [13], in
calculating the feature-specific $R^{2}$ for the above three models with data
sets of different sample sizes at $n=500,1000,2000$, and $5000$. We use
package sage-importance for SAGE and package vimpy for SPVIM to calculate
feature-specific Shapley values of total explained variance, which are divided
by the total variance for corresponding feature-specific $R^{2}$ values.

For each setting, we generated 1,000 data sets. For each data set, we built a
tree ensemble using XGBoost [21] with tuning parameters optimized via 5-fold
cross-validation and grid search in a parameter space specified with the
learning rate in $\\{0.01,0.05,0.1\\}$ and number of estimators in
$\\{50,100,200,300,\cdots,1000\\}$. We fixed the maximum depth of Models a, b,
and c at 1, 2, and 3 respectively. Figure 3 show the calculated feature-
specific $R^{2}$ for the first three features as well as the sum of all
feature-specific $R^{2}$ for all three models with $n=1000$, $p=100$, and
$\sigma_{\epsilon}=1.5$. The results of the three models in other settings are
shown in Appendix .3. Overall, Q-SHAP provides a more stable and accurate
calculation of feature-specific $R^{2}$ than the other two methods.

We divide all features into two groups, signal features (the first three) and
nuisance features (the rest). For each group, we calculated the mean absolute
error (MAE) by comparing feature-specific $R^{2}$ values to the theoretical
ones in each of the 1,000 datasets and averaged MAE over the 1,000 datasets,
shown in Figure 4. Note that, by limiting memory to 2GB, SAGE can only report
$R^{2}$ for the data sets with sample size at 500 and 1,000.

Figure 4: The mean absolute error (MAE) averaged over 1,000 datasets with
$p=100$

For both signal and nuisance features, Q-SHAP and SAGE exhibit consistent
behavior across all models. In contrast, SPVIM tends to bias the calculation,
especially for small sample sizes, indicated by the rapid increase of MMAE
when the sample size goes down. Among signal features, Q-SHAP has better
accuracy than SAGE, followed by SPVIM in general. All methods tend to have
better accuracy when sample size increases.

For the nuisance features, only SPVIM is biased away from 0. On the other
hand, both Q-SHAP and SAGE have almost no bias for nuisance features across
different sample sizes. For all three methods, $R^{2}$ of signal features
tends to have a larger bias than nuisance features.

Figure 5: The running time in $\log_{10}$(seconds)

We compared the computational time of the three different methods by running
all algorithms in parallel on a full node consisting of two AMD<EMAIL_ADDRESS>with 128 cores and 256 GB memory. We unified the environment with the help of
a Singularity container [24] built under Python version 3.11.6. Due to the
large size of the simulation, we limit all methods to a maximum wall time of 4
hours per dataset on a single core, with memory limited to 2 GB. The running
times are shown in Figure 5. Both SAGE and SPVIM demanded a long time to
compute even with only 100 features. Q-SHAP is hundreds of times faster than
both SAGE and SPVIM in general and is the only method that can be completed
when the dimension is 500 in constrained computation time and memory.

##  6 Real Data Analysis

We illustrate the utility of Q-SHAP by applying it to predicting Gleason
score, a grading prognosis of men with prostate cancer, with gene expressions.
The dataset was gerenated by The Cancer Genome Atlas Program (TCGA) [25],
including 551 samples and 17,261 features. The Gleason score was retrieved
through TCGAbiolinks [26] and further adjusted for age and race as potential
confounding factors. The gene expression data was downloaded from UCSC Xena
[27] and preprocessed using SIGNET [28].

We first constructed the tree ensemble using XGBoost with tuning parameters
optimized via 5-fold cross-validation and random research in a parameter space
specified with the number of trees in
$\\{50,100,500,1000,1500,2000,2500,3000\\}$, maximum depth in
$\\{1,2,\cdots,6\\}$, and learning rate in $\\{0.01,0.05,0.1\\}$. As both SAGE
and SPVIM cannot manage this large number of features, we only applied Q-SHAP
to the tree ensemble built on 17,261 features, which took 11 minutes to
compute and reported the sum of all feature-specific $R^{2}$ at $99.98\%$
which is equivalent to the model $R^{2}$. The 15 highest feature-specific
$R^{2}$ values are reported in Figure 6.

To allow the application of both SAGE and SPVIM, we selected the top 100
features based the result of Q-SHAP, and rebuilt the tree ensemble with the
selected 100 features using XGBoost. The rebuilt tree ensemble reported a
total of $R^{2}$ at $0.93$. Calculating the 100 feature-specific $R^{2}$
values took about 10 seconds of Q-SHAP, 16 minutes of SAGE, and 67 hours of
SPVIM. As shown in Figures 7-9, SPVIM tends to overstate the feature-specific
$R^{2}$ although all of its feature-specific $R^{2}$ only sums up to 0.72,
much lower than 0.93 reported by Q-SHAP. On the other hand, SAGE tends to
underestimate the feature-specific $R^{2}$ and all its feature-specific
$R^{2}$ also sums up to only 0.73, but its top 15 features match well with
those by Q-SHAP. In summary, the real data analysis shows consistent results
with the simulation study, confirming that Q-SHAP is superior in both
computational time and the accuracy of feature-specific $R^{2}$.

Figure 6: Top 15 of 17,261 feature-specific $R^{2}$ by Q-SHAP Figure 7: Top 15
of 100 feature-specific $R^{2}$ by Q-SHAP Figure 8: Top 15 of 100 feature-
specific $R^{2}$ by SAGE Figure 9: Top 15 of 100 feature-specific $R^{2}$ by
SPVIM

##  7 Conclusion

The coefficient of determination, aka $R^{2}$, measures the proportion of the
total variation explained by available features. Its additive decomposition,
following Shapley [5], provides an ideal evaluation of each feature’s
attribute to explaining the total variation. However, the calculation of
corresponding Shapley values is an NP-hard problem, and further complicated by
complexities involved in building tree ensembles. Recently, several methods
[3, 15, 16] have been developed to leverage the structure of tree-based models
and provide computationally efficient algorithms to decompose the predicted
values. However, decomposing $R^{2}$ demands the decomposition of a quadratic
loss reduction by multiple trees. We have shown in Section 4 that we can
attribute the total loss reduction by the tree ensemble to each single tree,
and the tree-specific loss reductions are subject to further decomposition to
each feature. However, decomposing a quadratic loss of a single tree need work
with the squared terms of predicted values, invalidating previously developed
methods for predicted values. Thus we developed Q-SHAP algorithm to calculate
Shapley values of squared terms of predicted values in polynomial time. The
algorithm works not only for $R^{2}$ but also for general quadratic losses.
Ultimately, it may provide a framework for more general loss functions via
approximation.

## Acknowledgment

We gratefully acknowledge the financial support from NIH grant R01GM131491 and
NCI Grant R03 CA235363. The results shown in Section 6 are in part based upon
data generated by the TCGA Research Network: https://www.cancer.gov/tcga.

## References

  * Stiglic et al. [2020] G. Stiglic, P. Kocbek, N. Fijacko, M. Zitnik, K. Verbert, and L. Cilar, “Interpretability of machine learning-based prediction models in healthcare,” _Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery_ , vol. 10, no. 5, p. e1379, 2020.
  * Bussmann et al. [2021] N. Bussmann, P. Giudici, D. Marinelli, and J. Papenbrock, “Explainable machine learning in credit risk management,” _Computational Economics_ , vol. 57, no. 1, pp. 203–216, 2021.
  * Lundberg and Lee [2017a] S. M. Lundberg and S.-I. Lee, “Consistent feature attribution for tree ensembles,” _Workshop on Human Interpretability in Machine Learning_ , 2017\.
  * Ishwaran [2007] H. Ishwaran, “Variable importance in binary regression trees and forests,” _Elecctronic Journal of Statistics_ , vol. 1, pp. 519–537, 2007.
  * Shapley [1953] L. S. Shapley, “A value for N-person games,” _Contributions to the Theory of Games_ , vol. 2(28), p. 307–317, 1953.
  * Lundberg and Lee [2017b] S. M. Lundberg and S.-I. Lee, “A unified approach to interpreting model predictions,” _Advances in Neural Information Processing Systems_ , vol. 30, 2017.
  * Chau et al. [2022] S. L. Chau, R. Hu, J. Gonzalez, and D. Sejdinovic, “RKHS-SHAP: Shapley values for kernel methods,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 13 050–13 063, 2022.
  * Molnar [2020] C. Molnar, _Interpretable Machine Learning_. lulu.com, 2020.
  * Covert et al. [2020] I. Covert, S. M. Lundberg, and S.-I. Lee, “Understanding global feature contributions with additive importance measures,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 17 212–17 223, 2020.
  * Lipovetsky and Conklin [2001] S. Lipovetsky and M. Conklin, “Analysis of regression in game theory approach,” _Applied Stochastic Models in Business and Industry_ , vol. 17, no. 4, pp. 319–330, 2001.
  * Owen and Prieur [2017] A. B. Owen and C. Prieur, “On Shapley value for measuring importance of dependent inputs,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 5, no. 1, pp. 986–1002, 2017.
  * Song et al. [2016] E. Song, B. L. Nelson, and J. Staum, “Shapley effects for global sensitivity analysis: theory and computation,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 4, no. 1, pp. 1060–1083, 2016.
  * Williamson and Feng [2020] B. Williamson and J. Feng, “Efficient nonparametric statistical inference on population feature importance using Shapley values,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 10 282–10 291.
  * Lundberg et al. [2020] S. M. Lundberg, G. Erion, H. Chen, A. DeGrave, J. M. Prutkin, B. Nair, R. Katz, J. Himmelfarb, N. Bansal, and S.-I. Lee, “From local explanations to global understanding with explainable ai for trees,” _Nature Machine Intelligence_ , vol. 2, no. 1, pp. 56–67, 2020.
  * Yang [2021] J. Yang, “Fast Treeshap: Accelerating SHAP value computation for trees,” _Advances in Neural Information Processing Systems_ , vol. 34, 2021.
  * Bifet et al. [2022] A. Bifet, J. Read, C. Xu _et al._ , “Linear TreeShap,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 25 818–25 828, 2022\.
  * Bénard et al. [2022] C. Bénard, G. Biau, S. Da Veiga, and E. Scornet, “SHAFF: Fast and consistent SHApley eFfect estimates via random Forests,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2022, pp. 5563–5582.
  * Breiman [2001] L. Breiman, “Random Forests,” _Machine Learning_ , vol. 45, pp. 5–32, 2001\.
  * Karczmarz et al. [2022] A. Karczmarz, T. Michalak, A. Mukherjee, P. Sankowski, and P. Wygocki, “Improved Feature Importance Computation for Tree Models based on the Banzhaf Value,” in _Proceedings of the Thirty-Eight Conference on Uncertainty in Artificial Intelligence_ , 2022, pp. 969–979.
  * Friedman [2001] J. H. Friedman, “Greedy function approximation: A gradient boosting machine,” _Annals of Statistics_ , pp. 1189–1232, 2001.
  * Chen and Guestrin [2016] T. Chen and C. Guestrin, “XGBoost: A scalable tree boosting system,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2016, pp. 785–794.
  * Ke et al. [2017] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu, “Lightgbm: A highly efficient gradient boosting decision tree,” _Advances in Neural Information Processing Systems_ , vol. 30, 2017.
  * Prokhorenkova et al. [2018] L. Prokhorenkova, G. Gusev, A. Vorobev, A. V. Dorogush, and A. Gulin, “Catboost: unbiased boosting with categorical features,” _Advances in Neural Information Processing Systems_ , vol. 31, 2018.
  * Kurtzer et al. [2017] G. M. Kurtzer, V. Sochat, and M. W. Bauer, “Singularity: Scientific containers for mobility of compute,” _PLOS ONE_ , vol. 12, no. 5, p. e0177459, 2017\.
  * Weinstein et al. [2013] J. N. Weinstein, E. A. Collisson, G. B. Mills, K. R. Shaw, B. A. Ozenberger, K. Ellrott, I. Shmulevich, C. Sander, and J. M. Stuart, “The cancer genome atlas pan-cancer analysis project,” _Nature Genetics_ , vol. 45, no. 10, pp. 1113–1120, 2013.
  * Colaprico et al. [2016] A. Colaprico, T. C. Silva, C. Olsen, L. Garofano, C. Cava, D. Garolini, T. S. Sabedot, T. M. Malta, S. M. Pagnotta, I. Castiglioni _et al._ , “TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data,” _Nucleic Acids Research_ , vol. 44, no. 8, p. e71, 2016.
  * Goldman et al. [2020] M. J. Goldman, B. Craft, M. Hastie, K. Repečka, F. McDade, A. Kamath, A. Banerjee, Y. Luo, D. Rogers, A. N. Brooks _et al._ , “Visualizing and interpreting cancer genomics data via the Xena platform,” _Nature Biotechnology_ , vol. 38, no. 6, pp. 675–678, 2020.
  * Jiang et al. [2023] Z. Jiang, C. Chen, Z. Xu, X. Wang, M. Zhang, and D. Zhang, “SIGNET: Transcriptome-wide causal inference for gene regulatory networks,” _Scientific Reports_ , vol. 13, no. 1, p. 19371, 2023.
  * Gould [1972] H. W. Gould, _Combinatorial Identities_. Morgantown Printing and Binding Co, 1972.
  * Geddes et al. [1992] K. O. Geddes, S. R. Czapor, and G. Labahn, _Algorithms for Computer Algebra_. Springer Science & Business Media, 1992.

###  .1 Proofs

We first establish the following lemma.

###### Lemma 1

$\sum_{k=0}^{p-n}\binom{|F|+k}{k}\binom{p-1-|F|-k}{p-n-k}=\binom{p}{n}.$

Proof. Using Gould’s identity [29], we have

$\sum_{k=0}^{p-n}\binom{|F|+k}{k}\binom{p-1-|F|-k}{p-n-k}=\sum_{k=0}^{p-n}\binom{k}{k}\binom{p-k-1}{p-n-k}=\sum_{k=0}^{p-n}\binom{p-k-1}{n-1}=\binom{p}{n},$

where we used the Hockey-Stick Identity in the last step. $\hfill\blacksquare$

Proof of Proposition 1. Through expansion and Lemma 1, we have

$\displaystyle\sum_{k=0}^{p-n}\frac{\binom{p-n}{k}\binom{n-1}{|F|}}{p\binom{p-1}{|F|+k}}$
$\displaystyle=\frac{1}{p}\sum_{k=0}^{p-n}\frac{(|F|+k)!(p-n)!(p-1-|F|-k)!(n-1)!}{(p-1)!(p-n-k)!k!|F|!(n-1-|F|!)}$
$\displaystyle=\frac{1}{p}\frac{(p-n)!(n-1)!}{(p-1)!}\sum_{k=0}^{p-n}\binom{|F|+k}{k}\binom{p-1-|F|-k}{p-n-k}$
$\displaystyle=\frac{1}{p}\frac{(p-n)!(n-1)!}{(p-1)!}\binom{p-1}{n-1}\frac{p}{n}=\frac{1}{n}.$

$\hfill\blacksquare$

Proof of Proposition 2. We first rewrite the two polynomials

$C(z)=V(z)c,\ \ \ P(z)=V(z)a,$

where $V(z)$ is the Vandermonde matrix for vector $z$, and $c$ and $a$ are the
coefficients of polynomials $C(z)$ and $P(z)$, respectively. Then the inner
product

$C(z)\cdot P(z)=P(z)\cdot
C(z)=a^{T}c=(V(z)^{-1}P(z))^{T}c=P(z)^{T}(V(z)^{T})^{-1}c.$

Letting

$z=\omega,$

and noting that the Vandermonde matrix evaluated at $\omega$ is symmetric, we
have

$(V(\omega)^{T})^{-1}=V(\omega)^{-1}=\frac{1}{n}V(\omega^{-1}),$

whose multiplication with $c$ is just the Inverse Fast Fourier transformation
(IFFT) over $c$ [30]. Hence the proposition holds. $\blacksquare$

###  .2 Calculation of Theoretical Shapley Values of $R^{2}$

Here we will calculate the theoretical feature-specific $R^{2}$ values of the
following three models,

$\displaystyle a.$ $\displaystyle
Y=\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{3}+\epsilon;$ $\displaystyle b.$
$\displaystyle
Y=\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{3}+\beta_{4}X_{1}X_{2}+\beta_{5}X_{1}X_{3}+\epsilon;$
$\displaystyle c.$ $\displaystyle
Y=\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{3}+\beta_{4}X_{1}X_{2}+\beta_{5}X_{1}X_{2}X_{3}+\epsilon.$

With all three features generated from Bernoulli distribution, we have

$\displaystyle E(X_{i})=E(X_{i}^{2})=p_{i}\triangleq\mu_{i},$ $\displaystyle
var(X_{i})=p_{i}(1-p_{i})\triangleq\sigma_{i}^{2},$ $\displaystyle
cov(X_{1},X_{1}X_{2})=\sigma_{1}^{2}\mu_{2},$ $\displaystyle
cov(X_{1}X_{2},X_{1}X_{3})=\sigma_{1}^{2}\mu_{2}\mu_{3},$ $\displaystyle
cov(X_{1},X_{1}X_{2}X_{3})=\sigma_{1}^{2}\mu_{2}\mu_{3},$ $\displaystyle
cov(X_{1}X_{2},X_{1}X_{2}X_{3})=\sigma_{12}^{2}\mu_{3},$ $\displaystyle
var(X_{1}X_{2})=p_{1}p_{2}(1-p_{1}p_{2})\triangleq\sigma_{12}^{2},$
$\displaystyle
var(X_{1}X_{2}X_{3})=p_{1}p_{2}p_{3}(1-p_{1}p_{2}p_{3})\triangleq\sigma_{123}^{2}.$

For Model a, we have

$\displaystyle
var(Y)=\beta_{1}^{2}\sigma_{1}^{2}+\beta_{2}^{2}\sigma_{2}^{2}+\beta_{3}^{2}\sigma_{3}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{1}))=\beta_{2}^{2}\sigma_{2}^{2}+\beta_{3}^{2}\sigma_{3}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{2}))=\beta_{1}^{2}\sigma_{1}^{2}+\beta_{3}^{2}\sigma_{3}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{3}))=\beta_{1}^{2}\sigma_{1}^{2}+\beta_{2}^{2}\sigma_{2}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{1,2\\}}))=\beta_{3}^{2}\sigma_{3}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{1,3\\}}))=\beta_{2}^{2}\sigma_{2}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{2,3\\}}))=\beta_{1}^{2}\sigma_{1}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle E(var(Y|X_{\\{1,2,3\\}}))=\sigma_{\epsilon}^{2}.$

For Model b, we have

$\displaystyle
E(var(Y|X_{\phi})=\beta_{1}^{2}\sigma_{1}^{2}+\beta_{2}^{2}\sigma_{2}^{2}+\beta_{3}^{2}\sigma_{3}^{2}+\beta_{4}^{2}\sigma_{12}^{2}+\beta_{5}^{2}\sigma_{13}^{2}+2\beta_{1}\beta_{4}\sigma_{1}^{2}\mu_{2}+2\beta_{1}\beta_{5}\sigma_{1}^{2}\mu_{3}+2\beta_{2}\beta_{4}\sigma_{2}^{2}\mu_{1}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+2\beta_{3}\beta_{5}\sigma_{3}^{2}\mu_{1}+2\beta_{4}\beta_{5}\sigma_{1}^{2}\mu_{2}\mu_{3}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{1})=(\beta_{2}^{2}+2\beta_{2}\beta_{4}\mu_{1}+\beta_{4}^{2}\mu_{1})\sigma_{2}^{2}+(\beta_{3}^{2}+2\beta_{3}\beta_{5}\mu_{1}+\beta_{5}^{2}\mu_{1})\sigma_{3}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{2})=(\beta_{1}^{2}+2\beta_{1}\beta_{4}\mu_{2}+\beta_{4}^{2}\mu_{2})\sigma_{1}^{2}+\beta_{3}^{2}\sigma_{3}^{2}+\beta_{5}^{2}\sigma_{13}^{2}+2(\beta_{1}+\beta_{4}\mu_{2})\beta_{5}\sigma_{1}^{2}\mu_{3}+2\beta_{3}\beta_{5}\sigma_{3}^{2}\mu_{1}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{3})=(\beta_{1}^{2}+2\beta_{1}\beta_{5}\mu_{3}+\beta_{5}^{2}\mu_{3})\sigma_{1}^{2}+\beta_{2}^{2}\sigma_{2}^{2}+\beta_{4}^{2}\sigma_{12}^{2}+2(\beta_{1}+\beta_{5}\mu_{3})\beta_{4}\sigma_{1}^{2}\mu_{2}+2\beta_{2}\beta_{4}\sigma_{2}^{2}\mu_{1}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{1,2\\}})=(\beta_{3}^{2}+2\beta_{3}\beta_{5}\mu_{1}+\beta_{5}^{2}\mu_{1})\sigma_{3}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{1,3\\}})=(\beta_{2}^{2}+2\beta_{2}\beta_{4}\mu_{1}+\beta_{4}^{2}\mu_{1})\sigma_{2}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{2,3\\}})=(\beta_{1}^{2}+\beta_{4}^{2}\mu_{2}+\beta_{5}^{2}\mu_{3}^{2}+2\beta_{1}\beta_{4}\mu_{2}+2\beta_{1}\beta_{5}\mu_{3}+2\beta_{4}\beta_{5}\mu_{2}\mu_{3})\sigma_{1}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle E(var(Y|X_{\\{1,2,3\\}})=\sigma_{\epsilon}^{2}$

For Model c, we have

$\displaystyle
var(Y)=\beta_{1}^{2}\sigma_{1}^{2}+\beta_{2}^{2}\sigma_{2}^{2}+\beta_{3}^{2}\sigma_{3}^{2}+\beta_{4}^{2}\sigma_{12}^{2}+\beta_{5}^{2}\sigma_{123}^{2}+2\beta_{1}\beta_{4}\sigma_{1}^{2}\mu_{2}+2\beta_{1}\beta_{5}\sigma_{1}^{2}\mu_{2}\mu_{3}+$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
2\beta_{2}\beta_{4}\sigma_{2}^{2}\mu_{1}+2\beta_{2}\beta_{5}\sigma_{2}^{2}\mu_{1}\mu_{3}+2\beta_{3}\beta_{5}\sigma_{3}^{2}\mu_{1}\mu_{2}+2\beta_{4}\beta_{5}\sigma_{12}^{2}\mu_{3}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{1})=(\beta_{2}^{2}+2\beta_{2}\beta_{4}\mu_{1}+\beta_{4}^{2}\mu_{1})\sigma_{2}^{2}+\beta_{3}^{2}\sigma_{3}^{2}+\beta_{5}^{2}\mu_{1}\sigma_{23}^{2}+2(\beta_{2}+\beta_{4})\beta_{5}\mu_{1}\sigma_{2}^{2}\mu_{3}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+2\beta_{3}\beta_{5}\mu_{1}\sigma_{3}^{2}\mu_{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{2})=(\beta_{1}^{2}+2\beta_{1}\beta_{4}\mu_{2}+\beta_{4}^{2}\mu_{2})\sigma_{1}^{2}+\beta_{3}^{2}\sigma_{3}^{2}+\beta_{5}^{2}\mu_{2}\sigma_{13}^{2}+2(\beta_{1}+\beta_{4})\beta_{5}\mu_{2}\sigma_{1}^{2}\mu_{3}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+2\beta_{3}\beta_{5}\mu_{2}\sigma_{3}^{2}\mu_{1}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{3}))=\beta_{1}^{2}\sigma_{1}^{2}+\beta_{2}^{2}\sigma_{2}^{2}+(\beta_{4}^{2}+2\beta_{4}\beta_{5}\mu_{3}+\beta_{5}^{2}\mu_{3})\sigma_{12}^{2}+2\beta_{1}(\beta_{4}+\beta_{5}\mu_{3})\sigma_{1}^{2}\mu_{2}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+2\beta_{2}(\beta_{4}+\beta_{5}\mu_{3})\sigma_{2}^{2}\mu_{1}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{1,2\\}})=(\beta_{3}^{2}+2\beta_{3}\beta_{5}\mu_{1}\mu_{2}+\beta_{5}^{2}\mu_{1}\mu_{2})\sigma_{3}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{1,3\\}})=(\beta_{2}^{2}+\beta_{4}^{2}\mu_{1}+\beta_{5}^{2}\mu_{1}\mu_{3}+2\beta_{2}\beta_{4}\mu_{1}+2\beta_{2}\beta_{5}\mu_{1}\mu_{3}+2\beta_{4}\beta_{5}\mu_{1}\mu_{3})\sigma_{2}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle
E(var(Y|X_{\\{2,3\\}})=(\beta_{1}^{2}+\beta_{4}^{2}\mu_{2}+\beta_{5}^{2}\mu_{2}\mu_{3}+2\beta_{1}\beta_{4}\mu_{2}+2\beta_{1}\beta_{5}\mu_{2}\mu_{3}+2\beta_{4}\beta_{5}\mu_{2}\mu_{3})\sigma_{1}^{2}+\sigma_{\epsilon}^{2},$
$\displaystyle E(var(Y|X_{\\{1,2,3\\}})=\sigma_{\epsilon}^{2}.$

For all of the cases the Shapley values can be calculated as

$\displaystyle\phi_{1}$ $\displaystyle=$
$\displaystyle-\frac{1}{var(Y)}\left\\{\frac{1}{3}[E(var(Y|X_{1}))-var(Y)]+\frac{1}{6}[E(var(Y|X_{\\{1,2\\}})-E(var(Y|X_{2}))]\right.$
$\displaystyle\left.+\frac{1}{6}[E(var(Y|X_{\\{1,3\\}})-E(var(Y|X_{3}))]+\frac{1}{3}[E(var(Y|X_{\\{1,2,3\\}})-E(V(var|X_{\\{2,3\\}}))]\right\\},$
$\displaystyle\phi_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{var(Y)}\left\\{\frac{1}{3}[E(var(Y|X_{2}))-var(Y)]+\frac{1}{6}[E(var(Y|X_{\\{1,2\\}})-E(var(Y|X_{1}))]\right.$
$\displaystyle\left.+\frac{1}{6}[E(var(Y|X_{\\{2,3\\}})-E(var(Y|X_{3}))]+\frac{1}{3}[E(var(Y|X_{\\{1,2,3\\}})-E(var(Y|X_{\\{1,3\\}}))]\right\\},$
$\displaystyle\phi_{3}$ $\displaystyle=$
$\displaystyle-\frac{1}{var(Y)}\left\\{\frac{1}{3}[E(var(Y|X_{3}))-var(Y)]+\frac{1}{6}[E(var(Y|X_{\\{1,3\\}})-E(var(Y|X_{1}))]\right.$
$\displaystyle\left.+\frac{1}{6}[E(var(Y|X_{\\{2,3\\}})-E(var(Y|X_{2}))]+\frac{1}{3}[E(var(Y|X_{\\{1,2,3\\}})-E(var(Y|X_{\\{1,2\\}}))]\right\\}.$

Therefore, the theoretical feature-specific $R^{2}$ in the three models can be
evaluated and are shown in Table 1.

Table 1: Theoretical Feature-Specific $R^{2}$ in Simulated Models |  | $R^{2}$
---|---|---
Model | $\sigma_{\epsilon}$ | Total | $X_{1}$ | $X_{2}$ | $X_{3}$
| 0.50 | 0.9864 | 0.2094 | 0.2863 | 0.4907
a | 1.00 | 0.9477 | 0.2012 | 0.2750 | 0.4715
| 1.50 | 0.8894 | 0.1888 | 0.2581 | 0.4425
| 0.50 | 0.9860 | 0.4390 | 0.1341 | 0.4129
b | 1.00 | 0.9459 | 0.4212 | 0.1286 | 0.3961
| 1.50 | 0.8860 | 0.3945 | 0.1205 | 0.3710
| 0.50 | 0.9868 | 0.4288 | 0.1450 | 0.4130
c | 1.00 | 0.9491 | 0.4124 | 0.1395 | 0.3972
| 1.50 | 0.8925 | 0.3878 | 0.1312 | 0.3735

###  .3 Boxplots of the First Three Feature-Specific and Total $R^{2}$ Values

We have compare the preformance of three different methods, i.e., our proposed
Q-SHAP, SAGE by Covert et al. [9], and SPVIM by Williamson and Feng [13], in
calculating the feature-specific $R^{2}$ as well as the sum of all feature-
specific $R^{2}$ for the three models specified in Section 5, with different
settings, i.e., $n\in\\{500,1000,2000,5000\\}$, $p\in\\{100,500\\}$, and
$\sigma_{\epsilon}\in\\{0.5,1,1.5\\}$. The results are shown in Fig. 10-17.
Note that the results of SAGE are unavailable in Fig. 12-17 because it cannot
report those $R^{2}$ with our limited computational resources, and the results
of SPVIM are unavailable in Fig. 14-17 because it demands too much time to
complete the computation when $p=500$.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 10: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=500$, $p=100$. The dashed lines show the theoretical $R^{2}$.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 11: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=1000$, $p=100$. The dashed lines show the theoretical $R^{2}$.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 12: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=2000$, $p=100$. The dashed lines show the theoretical $R^{2}$.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 13: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=5000$, $p=100$. The dashed lines show the theoretical $R^{2}$.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 14: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=500$, $p=500$. The dashed lines show the theoretical $R^{2}$.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 15: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=1000$, $p=500$. The dashed lines show the theoretical $R^{2}$.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 16: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=2000$, $p=500$. The dashed lines show the theoretical $R^{2}$.

(a) $X_{1}$-specific $R^{2}$ (b) $X_{2}$-specific $R^{2}$ (c) $X_{3}$-specific
$R^{2}$ (d) Sum of all $R^{2}$

Figure 17: Boxplots of (a) $X_{1}$-specific, (b) $X_{2}$-specific, (c)
$X_{3}$-specific, and (d) the sum of all feature-specific $R^{2}$ in the three
models with $n=5000$, $p=500$. The dashed lines show the theoretical $R^{2}$.

###  .4 Plots of the mean absolute error (MAE)

Similar to Fig. 4, we show in Fig. 18 the mean absolute error (MAE) of
feature-specific $R^{2}$ for both signal and nuisance features averaged over
1,000 datasets when $p=500$. Note that the results of SAGE and SPVIM are
unavailable because none of them can complete the computation for $p=500$ with
the limited computational resources.

Figure 18: The mean of absolute error (MAE) of the feature-specific $R^{2}$ by
Q-SHAP averaged across 1,000 datasets with $p=500$