\section{Introduction}
Complex machine learning models like deep neural networks and ensemble methods like random forest and gradient boosting machines may well outperform simpler approaches such as linear regression or single decision trees, but are noticeably harder to interpret. This can raise practical, ethical, and legal issues, most notably when applied in critical systems, e.g., for medical diagnosis or autonomous driving. The field of explainable AI aims to address these issues by enhancing the interpretability of complex machine learning models.

The Shapley-value approach has quickly become one of the most popular model-agnostic methods within explainable AI. It can provide local explanations, attributing changes in predictions for individual data points to the model's features, that can be combined to obtain better global understanding of the model structure~\cite{lundberg2020local}. Shapley values are based on a principled mathematical foundation~\cite{shapley1953value} and satisfy various desiderata (see also Section~\ref{sec:interpretation}). They have been applied for explaining statistical and machine learning models for quite some time, see e.g.,~\cite{lipovetsky2001analysis,vstrumbelj2014explaining}. Recent interests have been triggered by Lundberg and Lee's breakthrough paper~\cite{lundberg2017unified} that introduces efficient computational procedures and unifies Shapley values and other popular local model-agnostic approaches such as LIME~\cite{ribeiro2016should}.

Humans have a strong tendency to reason about their environment in causal terms~\cite{sloman2005causal}, where explanation and causation are intimately related: explanations often appeal to causes, and causal claims often answer questions about why or how something occurred~\cite{lombrozo2017causal}. The specific domain of causal responsibility studies how people attribute an effect to one or more causes, all of which may have contributed to the observed effect~\cite{sober1988apportioning}. Causal attributions by humans strongly depend on a subject's understanding of the generative model that explains how different causes lead to the effect, for which the relations between these causes are essential~\cite{gerstenberg2012noisy}.

Most explanation methods, however, tend to ignore such relations and act as if features are independent. Even so-called counterfactual approaches, that strongly rely on a causal intuition, make this simplifying assumption (e.g.,~\cite{wachter2017counterfactual}) and ignore that, in the real world, a change in one input feature may cause a change in another. This independence assumption also underlies early Shapley-based approaches, such as~\cite{vstrumbelj2014explaining,datta2016algorithmic}, and is made explicit as an approximation for computational reasons in~\cite{lundberg2017unified}. We will refer to these as {\em marginal} Shapley values.

Aas et al.~\cite{aas2019explaining} argue and illustrate that marginal Shapley values may lead to incorrect explanations when features are highly correlated, motivating what we will refer to as {\em conditional} Shapley values. Janzing et al.~\cite{janzing2019feature}, following~\cite{datta2016algorithmic}, discuss a causal interpretation of Shapley values, in which they replace conventional conditioning by observation with conditioning by intervention, as in Pearl's {\em do}-calculus~\cite{pearl2012calculus}. They argue that, when the goal is to causally explain the prediction {\em algorithm}, the inputs of this algorithm can be formally distinguished from the features in the real world and `interventional' Shapley values simplify to marginal Shapley values. This argument is also picked up by~\cite{lundberg2020local} when implementing interventional Tree SHAP. Going in a different direction, Frye et al.~\cite{frye2019asymmetric} propose {\em asymmetric} Shapley values as a way to incorporate causal knowledge in the real world by restricting the possible permutations of the features when computing the Shapley values to those consistent with a (partial) causal ordering. In line with~\cite{aas2019explaining}, they then apply conventional conditioning by observation to make sure that the explanations respect the data manifold.

The main contributions of our paper are as follows.
(1)~We derive {\em causal} Shapley values that explain the total effect of features on the prediction, taking into account their causal relationships. This makes them principally different from marginal and conditional Shapley values. At the same time, compared to asymmetric Shapley values, causal Shapley values provide a more direct and robust way to incorporate causal knowledge. (2)~Our method allows for further insights into feature relevance by separating out the total causal effect into a direct and indirect contribution. (3)~Making use of causal chain graphs~\cite{lauritzen2002chain}, we propose a practical approach for computing causal Shapley values and illustrate this on a real-world example.

\section{Causal Shapley values}
\label{sec:interpretation}

In this section, we will introduce causal Shapley values and contrast them to other approaches. We assume that we are given a machine learning model $f(\cdot)$ that can generate predictions for any feature vector $\mathbf{x}$. Our goal is to provide an explanation for an individual prediction $f(\mathbf{x})$, that takes into account the causal relationships between the features in the real world.

Attribution methods, with Shapley values as their most prominent example, provide a local explanation of individual predictions by attributing the difference between $f(\mathbf{x})$ and a baseline $f_0$ to the different features $i \in {N}$ with ${N} = \{1,\ldots,n\}$ and $n$ the number of features:
\begin{equation}
f(\mathbf{x}) = f_0 + \sum_{i=1}^n {\phi}_i \: ,
\label{eq:efficiency}
\end{equation}
where ${\phi}_i$ is the contribution of feature $i$ to the prediction $f(\mathbf{x})$. For the baseline $f_0$ we will take the average prediction $f_0 = \mathbb{E} f(\mathbf{X})$ with expectation taken over the observed data distribution $P(\mathbf{X})$. Equation~(\ref{eq:efficiency}) is referred to as the {\em efficiency property}~\cite{shapley1953value}, which appears to be a sensible desideratum for any attribution method and we therefore take here as our starting point.

To go from knowing none of the feature values, as for $f_0$, to knowing all feature values, as for $f(\mathbf{x})$, we can add feature values one by one, actively setting the features to their values in a particular order $\pi$. We define the contribution of feature $i$ given permutation $\pi$ as the difference in value function before and after setting the feature to its value:
\begin{equation}
{\phi}_i(\pi) = {v}(\{j: j \preceq_\pi i\}) - {v}(\{j: j \prec_\pi i\}) \: ,
\label{eq:contperm}
\end{equation}
with $j \prec_\pi i$ if $j$ precedes $i$ in the permutation $\pi$
and where we choose the value function
\begin{equation}
{v}(S) = \mathbb{E} \left[f(\mathbf{X}) | \lvdo{S} \right] = \int d\mathbf{X}_{\bar{S}} \: P(\mathbf{X}_{\bar{S}}|\lvdo{S}) f(\mathbf{X}_{\bar{S}},\mathbf{x}_S) \: .
\label{eq:valuedef}
\end{equation}
Here $S$ is the subset of `in-coalition' indices with known feature values $\mathbf{x}_S$. To compute the expectation, we average over the `out-of-coalition' or dropped features $\mathbf{X}_{\bar{S}}$ with $\bar{S} = {N} \setminus S$, the complement of $S$. To explicitly take into account possible causal relationships between the `in-coalition' features and the `out-of-coalition' features, we propose to condition `by intervention' for which we resort to Pearl's \textit{do}-calculus~\cite{pearl1995causal}. In words, the contribution ${\phi}_i(\pi)$ now measures the relevance of feature $i$ through the (average) prediction obtained if we actively set feature $i$ to its value $x_i$ compared to (the counterfactual situation of) not knowing its value.

Since the sum over features $i$ in~(\ref{eq:contperm}) is telescoping, the efficiency property~(\ref{eq:efficiency}) holds for any permutation $\pi$. Therefore, for any distribution over permutations $w(\pi)$ with $\sum_{\pi} w(\pi) = 1$, the contributions
\begin{equation}
{\phi}_i = \sum_{\pi} w(\pi) {\phi}_i(\pi)
\label{eq:shapperm}
\end{equation}
still satisfy~(\ref{eq:efficiency}). An obvious choice would be to take a uniform distribution $w(\pi) = 1/n!$. We then arrive at (with shorthand $i$ for the singleton $\{i\}$):
\[
{\phi}_i = \sum_{S \subseteq {N}\setminus i} \frac{|S|! (n-|S|-1)!}{n!} \left[{v}(S \cup i) - {v}(S) \right] \: .
\]
Besides efficiency, the Shapley values uniquely satisfy three other desirable properties~\cite{shapley1953value}.
\begin{description}
	\item[Linearity:] for two value functions ${v}_1$ and ${v}_2$, we have ${\phi}_i(\alpha_1 {v}_1 + \alpha_2 {v}_2) = \alpha_1 {\phi}_i({v}_1) + \alpha_2 {\phi}_i({v}_2)$. This guarantees that the Shapley value of a linear ensemble of models is a linear combination of the individual models' Shapley values.
	\item[Null player (dummy):] if ${v}(S \cup i) = {v}(S)$ for all $S \subseteq {N} \setminus i$, then ${\phi}_i = 0$. A feature that never contributes to the prediction (directly nor indirectly, see below) receives zero Shapley value.
	\item[Symmetry:] if ${v}(S \cup i) = {v}(S \cup j)$ for all  $S \subseteq {N} \setminus \{i,j\}$, then ${\phi}_i = {\phi}_j$. Symmetry in this sense holds for marginal, conditional, and causal Shapley values alike.
\end{description}
Note that symmetry here is defined w.r.t.\ to the contributions ${\phi}_i$, not the function values $f(\mathbf{x})$. As discussed in~\cite{janzing2019feature} (Section 3), conditioning by observation or intervention then does not break the symmetry property. For a non-uniform distribution of permutations as in~\cite{frye2019asymmetric}, symmetry is lost, but efficiency, linearity, and null player still apply.

Replacing conditioning by intervention with conventional conditioning by observation, i.e., averaging over $P(\mathbf{X}_{\bar{S}}|\mathbf{x}_{S})$ instead of $P(\mathbf{X}_{\bar{S}}|\lvdo{S})$ in~(\ref{eq:valuedef}), we arrive at the conditional Shapley values of~\cite{aas2019explaining,lundberg2018consistent}. A third option is to ignore the feature values $\mathbf{x}_S$ and take the unconditional, marginal distribution $P(\mathbf{X}_{\bar{S}})$, which leads to the marginal Shapley values.

Up until here, our active interventional interpretation of Shapley values coincides with that in~\cite{datta2016algorithmic,janzing2019feature,lundberg2020local}. However, from here on Janzing et al.~\cite{janzing2019feature} choose to ignore any dependencies between the features in the real world, by formally distinguishing between true features (corresponding to one of the data points) and the features plugged as input into the model. This leads to the conclusion that, in our notation, $P(\mathbf{X}_{\bar{S}}|\lvdo{S}) = P(\mathbf{X}_{\bar{S}})$ for any subset $S$. As a result, any expectation under conditioning by intervention collapses to a marginal expectation and, in the interpretation of~\cite{datta2016algorithmic,janzing2019feature,lundberg2020local}, interventional Shapley values simplify to marginal Shapley values. As we will see below, marginal Shapley values can only represent direct effects, which makes that `root causes' with strong indirect effects (e.g.\ genetic markers) are ignored in the attribution, which is quite different from how humans tend to attribute causes~\cite{sober1988apportioning}. In this paper, we choose not to make this distinction between the features in the real world and the inputs of the prediction model, but to explicitly take into account the causal relationships between the data in the real world to enhance the explanations. Since the term `interventional' Shapley values has been coined for causal explanations of the prediction algorithm, ignoring causal relationships between the features in the real world, we will use the term `causal' Shapley values for those that do attempt to incorporate these relationships using Pearl's {\em do}-calculus.

The asymmetric Shapley values introduced in~\cite{frye2019asymmetric} have the same objective: enhancing the explanation of the Shapley values by incorporating causal knowledge about the features in the real world. In ~\cite{frye2019asymmetric}, this knowledge is incorporated by choosing $w(\pi) \neq 0$ in~(\ref{eq:shapperm}) only for those permutations $\pi$ that are consistent with the causal structure between the features, i.e., are such that a known causal ancestor always precedes its descendants. On top of this, Frey et al.~\cite{frye2019asymmetric} apply standard conditioning by observation. In this paper we show that there is no need to resort to asymmetric Shapley values to incorporate causal knowledge: applying conditioning by intervention instead of conditioning by observation is sufficient. Choosing asymmetric Shapley values instead of symmetric ones can be considered orthogonal to choosing conditioning by observation versus conditioning by intervention. We will therefore refer to the approach of~\cite{frye2019asymmetric} as {\em asymmetric conditional} Shapley values, to contrast them with {\em asymmetric causal} Shapley values that implement both ideas.

\section{Decomposing Shapley values into direct and indirect effects}

Our causal interpretation allows us to distinguish between direct and indirect effects of each feature on a model's prediction. This decomposition then also helps to understand the difference between marginal, symmetric, and asymmetric Shapley values. Going back to the contribution ${\phi}_i(\pi)$ for a permutation $\pi$ and feature $i$ in (\ref{eq:contperm}) and using shorthand notation $\underline{S} = \{j: j \prec_\pi i\}$ and $\bar{S} = \{j: j \succ_\pi i\}$, we can decompose the total effect for this permutation into a direct and an indirect effect:
\begin{eqnarray}
	{\phi}_i(\pi) \hspace*{-2.5mm} & = & \hspace*{-2.5mm}
	\mathbb{E}[f(\mathbf{X}_{\bar{S}},\mathbf{x}_{\underline{S} \cup i})|\lvdo{\underline{S} \cup i}] - \mathbb{E}[f(\mathbf{X}_{\bar{S} \cup i},\mathbf{x}_{\underline{S}})|\lvdo{\underline{S}}] ~~~~~~\mbox{(total effect)} \nonumber \\
	\hspace*{-2.5mm} & = & \hspace*{-2.5mm} \mathbb{E}[f(\mathbf{X}_{\bar{S}},\mathbf{x}_{\underline{S} \cup i})|\lvdo{\underline{S}}] - \mathbb{E}[f(\mathbf{X}_{\bar{S} \cup i},\mathbf{x}_{\underline{S}})|\lvdo{\underline{S}}] + ~~~~~~~\mbox{(direct effect)} \nonumber \\
	&& \! \! \! \! \mathbb{E}[f(\mathbf{X}_{\bar{S}},\mathbf{x}_{\underline{S} \cup i})|\lvdo{\underline{S} \cup i}] - \mathbb{E}[f(\mathbf{X}_{\bar{S}},\mathbf{x}_{\underline{S} \cup i})|\lvdo{\underline{S}}] ~~~~\mbox{(indirect effect)} \nonumber \\ && ~~\label{eq:decomposition}
\end{eqnarray}
The direct effect measures the expected change in prediction when the stochastic feature $X_i$ is replaced by its feature value $x_i$, without changing the distribution of the other `out-of-coalition' features. The indirect effect measures the difference in expectation when the distribution of the other `out-of-coalition' features changes due to the additional intervention $\ldo{i}$. The direct and indirect parts of Shapley values can then be computed as in~(\ref{eq:shapperm}): by taking a, possibly weighted, average over all permutations. Conditional Shapley values can be decomposed similarly by replacing conditioning by intervention with conditioning by observation in~(\ref{eq:decomposition}). For marginal Shapley values, there is no conditioning and hence no indirect effect: by construction marginal Shapley values can only represent direct effects. We will make use of this decomposition in the next section to clarify how different causal structures lead to different Shapley values.

\section{Shapley values for different causal structures}
\label{sec:toymodels}

\newcommand{{\em D}}{{\em D}}
\newcommand{{\em E}}{{\em E}}
\newcommand{{\em R}}{{\em R}}
\newcommand{\stupid}[1]{\color{red}\underline{#1}}
\begin{figure}
	\centering
	\begin{tabular}{c|cc|cc|cc}
		& \multicolumn{2}{c|}{{\em D}} & \multicolumn{2}{c|}{{\em E}} & \multicolumn{2}{c}{{\em R}} \\[0.3em] 
		& direct & indirect & direct & indirect & direct & indirect \\ \midrule
		$\phi_1$ & 0 & 0 & 0 & $\frac{1}{2} \beta \alpha x_1$ & 0 & $\beta \alpha x_1$ \\
		$\phi_2$ & $\beta x_2$ & 0 & $\beta x_2 - \frac{1}{2} \beta \alpha x_1$ & 0 & $\beta x_2 - \beta \alpha x_1$ & 0 \\ \bottomrule
	\end{tabular}\\
	\vspace*{2em}
	\begin{minipage}{0.45\textwidth}
		\centering
		\tikzstyle{arrow} = [thick,->,>=stealth]
		\tikzstyle{dashedarrow} = [thick,->,>=stealth,dashed]
		\tikz{
		
			\node[latent] at (0,0) (y1) {$Y$};
			\node[latent] at (0,1.5) (x12) {$X_2$};
			\node[latent] at (0,3) (x11) {$X_1$};
			\node[text height=1em, anchor=north] at (0,4.2) {Chain};
			\draw[arrow] (x12) -- node[anchor=west]{$\beta$} (y1);
			\draw[arrow] (x11) -- node[anchor=west]{$\alpha$} (x12);
			\draw[dashedarrow] (x11) to[bend right] (y1);
		
			\node[latent] at (1.5,0) (y2) {$Y$};
			\node[latent] at (1.5,1.5) (x22) {$X_2$};
			\node[latent] at (1.5,3) (x21) {$X_1$};
			\node[text height=1em, anchor=north] at (1.5,4.2) {Fork};
			\draw[arrow] (x22) -- node[anchor=west]{$\beta$} (y2);
			\draw[arrow] (x22) -- (x21);
			\draw[dashedarrow] (x21) to[bend right] (y2);
		
			\node[latent] at (3,0) (y3) {$Y$};
			\node[latent] at (3,1.5) (x32) {$X_2$};
			\node[latent] at (3,3) (x31) {$X_1$};
			\node[latent] at (3.6,2.25) (z)  {$Z$};
			\node[text height=1em, anchor=north] at (3,4.2) {Confounder};
			\draw[arrow] (x32) -- node[anchor=west]{$\beta$} (y3);
			\draw[arrow] (z) -- (x31);
			\draw[arrow] (z) -- (x32);
			\draw[dashedarrow] (x31) to[bend right] (y3);
		
			\node[latent] at (4.5,0) (y4) {$Y$};
			\node[latent] at (4.5,1.5) (x42) {$X_2$};
			\node[latent] at (4.5,3) (x41) {$X_1$};
			\node[text height=1em, anchor=north] at (4.5,4.2) {Cycle};
			\draw[arrow] (x42) -- node[anchor=west]{$\beta$} (y4);
			\draw[arrow] (x41) to[bend left=15] (x42);
			\draw[arrow] (x42) to[bend left=15] (x41);
			\draw[dashedarrow] (x41) to[bend right] (y4);			
		}
	\end{minipage}
	\hfill
	\begin{minipage}{0.53\textwidth}
		\begin{tabular}{r|c|cc|cc}
			& marginal & \multicolumn{2}{c|}{conditional} & \multicolumn{2}{c}{causal} \\[0.3em] 
			& & \rotatebox{90}{symmetric} & \rotatebox{90}{asymmetric} & \rotatebox{90}{symmetric} & \rotatebox{90}{asymmetric} \\ \midrule
			Chain		& \stupid{{\em D}} & {\em E} & {\em R} & {\em E} & {\em R} \\
			Fork		& {\em D} & \stupid{{\em E}} & {\em D} & {\em D} & {\em D} \\
			Confounder 	& {\em D} & \stupid{{\em E}} & \stupid{{\em E}} & {\em D} & {\em D} \\
			Cycle		& \stupid{{\em D}} & {\em E} & {\em E} & {\em E} & {\em E} \\
			\bottomrule
		\end{tabular}
	\end{minipage}
	\caption{Direct and indirect Shapley values for four causal models with the same observational distribution over features (such that $\mathbb{E}[X_1] = \mathbb{E}[X_2] = 0$ and $\mathbb{E}[X_2|x_1] = \alpha x_1$), yet a different causal structure. We assume a linear model that happens to ignore the first feature: $f(x_1,x_2) = \beta x_2$. The bottom table gives for each of the four causal models on the left the marginal, conditional, and causal Shapley values, where the latter two are further split up in symmetric and asymmetric. Each letter in the bottom table corresponds to one of the patterns of direct and indirect effects detailed in the top table: `direct' ({\em D}, only direct effects), `evenly split' ({\em E}, credit for an indirect effect split evenly between the features), and `root cause' ({\em R}, all credit for the indirect effect goes to the root cause). Shapley values that, we argue, do not provide a proper causal explanation, are underlined and indicated in red.}
	\label{fig:fourmodels}
\end{figure}

To illustrate the difference between the various Shapley values, we consider four causal models on two features. They are constructed such that they have the same $P(\mathbf{X})$, with $\mathbb{E}[X_2|x_1] = \alpha x_1$ and $\mathbb{E}[X_1] = \mathbb{E}[X_2] = 0$, but with different causal explanations for the dependency between $X_1$ and $X_2$. In the causal chain $X_1$ could, for example, represent season, $X_2$ temperature, and $Y$ bike rental. The fork inverts the arrow between $X_1$ and $X_2$, where now $Y$ may represent hotel occupation, $X_2$ season, and $X_1$ temperature. In the chain and the fork, different data points correspond to different days. For the confounder and the cycle, $X_1$ and $X_2$ may represent obesity and sleep apnea, respectively, and $Y$ hours of sleep. The confounder model implements the assumption that obesity and sleep apnea have a common confounder $Z$, e.g., some genetic predisposition. The cycle, on the other hand, represents the more common assumption that there is a reciprocal effect, with obesity affecting sleep apnea and vice versa~\cite{ong2013reciprocal}. In the confounder and the cycle, different data points correspond to different subjects. We assume to have trained a linear model $f(x_1,x_2)$ that happens to largely, or even completely to simplify the formulas, ignore the first feature, and boils down to the prediction function $f(x_1,x_2) = \beta x_2$. Figure~\ref{fig:fourmodels} shows the explanations provided by the various Shapley values for each of the causal models in this extreme situation. Derivations can be found in the supplement.

Since in all cases there is no direct link between $X_1$ and the prediction, the direct effect of $X_1$ is always equal to zero. Similarly, any indirect effect of $X_2$ can only go through $X_1$ and hence must also be equal to zero. So, all we can expect is a direct effect from $X_2$, proportional to $\beta$, and an indirect effect from $X_1$ through $X_1$, proportional to $\alpha$ times $\beta$. Because of the sufficiency property, the direct and indirect effect always add up to the output $\beta x_2$. This makes that, for all different combinations of causal structures and types of Shapley values, we end up with just three different explanation patterns, referred to as {\em D}, {\em E}, and {\em R} in Figure~\ref{fig:fourmodels}.

To argue which explanations make sense, we call upon classical norm theory~\cite{kahneman1986norm}. It states that humans, when asked for an explanation of an effect, contrast the actual observation with a counterfactual, more normal alternative. What is considered normal, depends on the context. Shapley values can be given the same interpretation~\cite{merrick2019explanation}: they measure the difference in prediction between knowing and not knowing the value of a particular feature, where the choice of what's normal translates to the choice of the reference distribution to average over when the feature value is still unknown.

In this perspective, marginal Shapley values as in~\cite{datta2016algorithmic,janzing2019feature,lundberg2020local} correspond to a very simplistic, counterintuitive interpretation of what's normal. Consider for example the case of the chain, with $X_1$ representing season, $X_2$ temperature, and $Y$ bike rental, and two days with the same temperature of 13 degrees Celsius, one in fall and another in winter. Marginal Shapley values end up with the same explanation for the predicted bike rental on both days, ignoring that the temperature in winter is higher than normal for the time of year and in fall lower. Just like marginal Shapley values, symmetric conditional Shapley values as in~\cite{aas2019explaining} do not distinguish between any of the four causal structures. They do take into account the dependency between the two features, but then fail to acknowledge that an {\em intervention} on feature $X_1$ in the fork and the confounder, does not change the distribution of $X_2$.

For the confounder and the cycle, asymmetric Shapley values put $X_1$ and $X_2$ on an equal footing and then coincide with their symmetric counterparts. Asymmetric {\em conditional} Shapley values from~\cite{frye2019asymmetric} have no means to distinguish between the cycle and the confounder, unrealistically assigning credit to $X_1$ in the latter case. Asymmetric and symmetric {\em causal} Shapley values do correctly treat the cycle and confounder cases.

In the case of a chain, asymmetric and symmetric causal Shapley values provide different explanations. Which explanation is to be preferred may well depend on the context. In our bike rental example, asymmetric Shapley values first give full credit to season for the indirect effect (here $\alpha \beta x_1$), subtracting this from the direct effect of the temperature to fulfill the sufficiency property ($\beta x_2 - \alpha \beta x_1$). Symmetric causal Shapley values consider both contexts -- one in which season is intervened upon before temperature, and one in which temperature is intervened upon before season -- and then average over the results in these two contexts. This symmetric strategy appears to better appeal to the theory dating back to~\cite{lewis1974causation}, that humans sample over different possible scenarios (here: different orderings of the features) to judge causation. However, when dealing with a temporal chain of events, alternative theories (see e.g.~\cite{spellman1997crediting}) suggest that humans have a tendency to attribute credit or blame foremost to the root cause, which seems closer in spirit to the explanation provided by asymmetric causal Shapley values.

By dropping the symmetry property, asymmetric Shapley values do pay a price: they are sensitive to the insertion of causal links with zero strength. As an example, consider a neural network trained to perfectly predict the XOR function on two binary variables $X_1$ and $X_2$. With a uniform distribution over all features and no further assumption w.r.t.\ the causal ordering of $X_1$ and $X_2$, the Shapley values are $\phi_1 = \phi_2 = 1/4$ when the prediction $f$ equals $1$, and $\phi_1 = \phi_2 = -1/4$ for $f=0$: completely symmetric. If we now assume that $X_1$ preceeds $X_2$ (and a causal strength of 0 to maintain the uniform distribution over features), all Shapley values stay the same, except for the asymmetric ones: these suddenly jump to $\phi_1 = 0$ and $\phi_2 = 1/2$ for $f=1$, and $\phi_1 = 0$ and $\phi_2 = -1/2$ for $f=0$. More details on this instability of asymmetric Shapley values can be found in the supplement, where we compare Shapley values of trained neural networks for varying causal strengths.

To summarize, unlike marginal and (both symmetric and asymmetric) conditional Shapley values, causal Shapley values provide sensible explanations that incorporate causal relationships in the real world. Asymmetric causal Shapley values may be preferable over symmetric ones when causality derives from a clear temporal order, whereas symmetric Shapley values have the advantage of being much less sensitive to model misspecifications.

\section{A practical implementation with causal chain graphs}

In the ideal situation, a practitioner has access to a fully specified causal model that can be plugged in~(\ref{eq:valuedef}) to compute or sample from every interventional probability of interest. In practice, such a requirement is hardly realistic. In fact, even if a practitioner could specify a complete causal structure (including potential confounding) and has full access to the observational probability $P(\mathbf{X})$, not every causal query need be identifiable (see e.g., \cite{pearl2012calculus}). Furthermore, requiring so much prior knowledge could be detrimental to the method's general applicability. In this section, we describe a pragmatic approach that is applicable when we have access to a (partial) causal ordering plus a bit of additional information to distinguish confounders from mutual interactions, and a training set to estimate (relevant parameters of) $P(\mathbf{X})$. Our approach is inspired by~\cite{frye2019asymmetric}, but extends it in various aspects: it provides a formalization in terms of causal chain graphs, applies to both symmetric and asymmetric Shapley values, and correctly distinguishes between dependencies that are due to confounding and mutual interactions.

In the special case that a complete causal ordering of the features can be given and that all causal relationships are unconfounded, $P(\mathbf{X})$ satisfies the Markov properties associated with a directed acyclic graph (DAG) and can be written in the form
\[
P(\mathbf{X}) = \prod_{j \in {N}} P(X_j|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(j)}) \: ,
\]
with $\mathop{\textit{pa}}(j)$ the parents of node $j$. With no further conditional independences, the parents of $j$ are all nodes that precede $j$ in the causal ordering. For causal DAGs, we have the interventional formula~\cite{lauritzen2002chain}:
\begin{equation}
P(\mathbf{X}_{\bar{S}}|\lvdo{S}) = \prod_{j \in \bar{S}} P(X_j|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(j)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(j) \cap S}) \: ,
\label{eq:interventional}
\end{equation}
with $\mathop{\textit{pa}}(j) \cap T$ the parents of $j$ that are also part of subset $T$. The interventional formula can be used to answer any causal query of interest.

\begin{figure}
	\centering
	\tikz{
	
		\node (t1) at (-3,2) [align=left]{partial causal ordering:\\[1em]
		$\left\{(1,2),(3,4,5),(6,7)\right\}$};
		\node (t2) at (-0.5,2) {\Huge $\Rightarrow$};	
		\node (t3) at (4.5,2) {\Huge $\Rightarrow$};		
		\node (1) at (0,4) [circle,draw,inner sep=0.5mm] {\scriptsize $X_1$};
		\node (2) at (1,4) [circle,draw,inner sep=0.5mm] {\scriptsize $X_2$};
		\node (3) at (1.9,2.8) [circle,draw,inner sep=0.5mm] {\scriptsize $X_3$};
		\node (4) at (3.1,2.8) [circle,draw,inner sep=0.5mm] {\scriptsize $X_4$};
		\node (5) at (2.5,1.8) [circle,draw,inner sep=0.5mm] {\scriptsize $X_5$};
		\node (6) at (1,0.5) [circle,draw,inner sep=0.5mm] {\scriptsize $X_6$};
		\node (7) at (2,0.5) [circle,draw,inner sep=0.5mm] {\scriptsize $X_7$};
		\node (8) [draw,label=above:$\tau_1$,inner sep=0.5mm,fit=(1) (2)] {};
		\node (9) [draw,label=above:$\tau_2$,inner sep=0.5mm,fit=(3) (4) (5)] {};
		\node (10) [draw,label=above:$~~\tau_3$,inner sep=0.5mm,fit=(6) (7)] {};
		\draw (1) -- (2);
		\draw (3) -- (4) -- (5) -- (3);
		\draw (6) -- (7);
		\edge {8} {9};
		\edge {8,9} {10};
		\node (11) at (5,4) [circle,draw,inner sep=0.5mm] {\scriptsize $X_1$};
		\node (12) at (6,4) [circle,draw,inner sep=0.5mm] {\scriptsize $X_2$};
		\node (13) at (6.9,2.8) [circle,draw,inner sep=0.5mm] {\scriptsize $X_3$};
		\node (14) at (8.1,2.8) [circle,draw,inner sep=0.5mm] {\scriptsize $X_4$};
		\node (15) at (7.5,1.8) [circle,draw,inner sep=0.5mm] {\scriptsize $X_5$};
		\node (16) at (6,0.5) [circle,draw,inner sep=0.5mm] {\scriptsize $X_6$};
		\node (17) at (7,0.5) [circle,draw,inner sep=0.5mm] {\scriptsize $X_7$};
		\node (21) at (5.5,4.8) [circle,draw,inner sep=0.7mm]
		{\scriptsize $Z_1$};
		\node (22) at (6.5,-0.3) [circle,draw,inner sep=0.7mm]
		{\scriptsize $Z_2$};
		\node (18) [draw,label=above:$\tau_1$,inner sep=0.5mm,fit=(11) (12) (21)] {};
		\node (19) [draw,label=above:$\tau_2$,inner sep=0.5mm,fit=(13) (14) (15)] {};
		\node (20) [draw,label=above:$~~\tau_3$,inner sep=0.5mm,fit=(16) (17) (22)] {};
		\edge {18} {19};
		\edge {18,19} {20};
		\draw[-Latex] (13) to[bend right=15] (14);
		\draw[-Latex] (14) to[bend right=15] (13);
		\draw[-Latex] (13) to[bend right=15] (15);
		\draw[-Latex] (15) to[bend right=15] (13);
		\draw[-Latex] (15) to[bend right=15] (14);
		\draw[-Latex] (14) to[bend right=15] (15);
		\draw[-Latex] (22) to (16);
		\draw[-Latex] (22) to (17);
		\draw[-Latex] (21) to (11);
		\draw[-Latex] (21) to (12);
	}	
	\caption{From partial ordering to causal chain graph. Features on an equal footing are combined into a fully connected chain component. How to handle interventions within each component depends on the generative process that best explains the (surplus) dependencies. In this example, the dependencies in chain components $\tau_1$ and $\tau_3$ are assumed to be the result of a common confounder, and those in $\tau_2$ of mutual interactions.}
	\label{fig:chaingraph}
\end{figure}

When we cannot give a complete ordering between the individual variables, but still a partial ordering, causal chain graphs~\cite{lauritzen2002chain} come to the rescue. A causal chain graph has directed and undirected edges. All features that are treated on an equal footing are linked together with undirected edges and become part of the same chain component. Edges between chain components are directed and represent causal relationships. See Figure~\ref{fig:chaingraph} for an illustration of the procedure. The probability distribution $P(\mathbf{X})$ in a chain graph factorizes as a ``DAG of chain components'':
\[
P(\mathbf{X}) = \prod_{\tau \in {\cal T}} P(\mathbf{X}_\tau|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)}) \: ,
\]
with each $\tau$ a chain component, consisting of all features that are treated on an equal footing.

How to compute the effect of an intervention depends on the interpretation of the generative process leading to the (surplus) dependencies between features within each component. If we assume that these are the consequence of marginalizing out a common confounder, intervention on a particular feature will break the dependency with the other features. We will refer to the set of chain components for which this applies as $\onder{{\cal T}}{confounding}$. The undirected part can also correspond to the equilibrium distribution of a dynamic process resulting from interactions between the variables within a component~\cite{lauritzen2002chain}. In this case, setting the value of a feature does affect the distribution of the variables within the same component. We refer to these sets of components as ${\cal T}_{\overline{\textrm{\upshape \scriptsize confounding}}}$.

Any expectation by intervention needed to compute the causal Shapley values can be translated to an expectation by observation, by making use of the following theorem (see the supplement for a more detailed proof and some corollaries linking back to other types of Shapley values as special cases).
\begin{theorem}
For causal chain graphs, we have the interventional formula 
\begin{eqnarray}
P(\mathbf{X}_{\bar{S}}|\lvdo{S}) \hspace*{-2.5mm} & = & \hspace*{-2.5mm} \prod_{\tau \in {\cal T}_{\textrm{\upshape \scriptsize confounding}}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S}) \times \nonumber \\
&& \prod_{\tau \in {\cal T}_{\overline{\textrm{\upshape \scriptsize confounding}}}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\mathbf{x}_{\tau \cap S}) \: .
\label{eq:chaininterventional}
\end{eqnarray}
\end{theorem}

\begin{proof}
\begin{eqnarray*}
	P(\mathbf{X}_{\bar{S}}|\lvdo{S}) \isequaldo{1} \prod_{\tau \in {\cal T}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\lvdo{S}) \hfill \\
	\isequaldo{3}  \prod_{\tau \in {\cal T}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\lvdo{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\lvdo{\tau \cap S}) \\
	\isequaldo{2}  \prod_{\tau \in {\cal T}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\lvdo{\tau \cap S}) \: ,
\end{eqnarray*}
where the number above each equal sign refers to the standard \textit{do}-calculus rule from~\cite{pearl2012calculus} that is applied. For a chain component with dependencies induced by a common confounder, rule (3) applies once more and yields $P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S})$,
whereas for a chain component with dependencies induced by mutual interactions, rule (2) again applies and gives $P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{ \mathop{\textit{\scriptsize pa}}(\tau) \cap S},\mathbf{x}_{\tau \cap S}))$.
\end{proof}

To compute these observational expectations, we can rely on the various methods that have been proposed to compute conditional Shapley values~\cite{aas2019explaining,frye2019asymmetric}. Following~\cite{aas2019explaining}, we will assume a multivariate Gaussian distribution for $P(\mathbf{X})$ that we estimate from the training data. Alternative proposals include assuming a Gaussian copula distribution, estimating from the empirical (conditional) distribution (both from~\cite{aas2019explaining}) and a variational autoencoder~\cite{frye2019asymmetric}.

\section{Illustration on real-world data}

\begin{figure}[t]
	\centering
	\begin{minipage}{.49\linewidth}
		\includegraphics[width=\textwidth]{trend_plot.pdf}
		\includegraphics[width=\textwidth]{corr_plots_top.pdf}
		\includegraphics[width=\textwidth]{corr_plots_bottom.pdf}
	\end{minipage}
	\begin{minipage}{.5\linewidth}
		\vfill
		\includegraphics[width=\textwidth]{sina_plot.pdf}
	\end{minipage}
	\caption{Bike shares in Washington, D.C.\ in 2011-2012 (top left; colorbar with temperature in degrees Celsius). Sina plot of causal Shapley values for a trained XGBoost model, where the top three date-related variables are considered to be a potential cause of the four weather-related variables (right). Scatter plots of marginal (MSV) versus causal Shapley values (CSV) for temperature ({\em temp}) and one of the seasonal variables ({\em cosyear}) show that MSVs almost purely explain the predictions based on temperature, whereas CSVs also give credit to season (bottom left).}
	\label{fig:trendplot}
\end{figure}

To illustrate the difference between marginal and causal Shapley values, we consider the bike rental dataset from~\cite{fanaee2013bikerental}, where we take as features the number of days since January 2011 ({\em trend}), two cyclical variables to represent season ({\em cosyear}, {\em sinyear}), the temperature ({\em temp}), feeling temperature ({\em atemp}), wind speed ({\em windspeed}), and humidity ({\em hum}). As can be seen from the time series itself (top left plot in Figure~\ref{fig:trendplot}), the bike rental is strongly seasonal and shows an upward trend. Data was randomly split in 80\% training and 20\% test set. We trained an XGBoost model for 100 rounds. 

We adapted the R package SHAPR from~\cite{aas2019explaining} to compute causal Shapley values, which essentially boiled down to an adaptation of the sampling procedure so that it draws samples from the interventional conditional distribution~(\ref{eq:chaininterventional}) instead of from a conventional observational conditional distribution. The sina plot on the righthand side of Figure~\ref{fig:trendplot} shows the causal Shapley values calculated for the trained XGBoost model on the test data. For this simulation, we chose the partial order $(\{\textit{trend}\},\{\textit{cosyear},\textit{sinyear}\},\{\textrm{all weather variables}\})$, with confounding for the second component and no confounding for the third, to represent that season has an effect on weather, but that we have no clue how to represent the intricate relations between the various weather variables. The sina plot clearly shows the relevance of the trend and the season (in particular cosine of the year, which is -1 on January 1 and +1 on July 1). The scatter plots on the left zoom in on the causal (CSV) and marginal Shapley values (MSV) for {\em cosyear} and {\em temp.} The marginal Shapley values for {\em cosyear} vary over a much smaller range than the causal Shapley values for {\em cosyear}, and vice versa for the Shapley values for {\em temp}: where the marginal Shapley values explain the predictions predominantly based on temperature, the causal Shapley values give season much more credit for the higher bike rental in summer and the lower bike rental in winter. Sina plots for marginal and asymmetric conditional Shapley values can be found in the supplement.

\begin{wrapfigure}[21]{r}{0.5\textwidth}
	\begin{center}
		\includegraphics[width=0.48\textwidth]{final_bar_plots.pdf}
	\end{center}
	\caption{Asymmetric (conditional), (symmetric) causal and marginal Shapley values for two different days, one in October (brown) and one in December (gray) with more or less the same temperature of 13 degrees Celsius. Asymmetric Shapley values focus on the root cause, marginal Shapley values on the more direct effect, and symmetric causal Shapley consider both for the more natural explanation.}
	\label{fig:barplots}
\end{wrapfigure}
The difference between asymmetric (conditional, from~\cite{frye2019asymmetric}), (symmetric) causal, and marginal Shapley values clearly shows when we consider two days, October 10 and December 3, 2012, with more or less the same temperature of 13 and 13.27 degrees Celsius, and predicted bike counts of 6117 and 6241, respectively. The temperature and predicted bike counts are relatively low for October, yet high for December. The various Shapley values for {\em cosyear} and {\em temp} are shown in Figure~\ref{fig:barplots}. The marginal Shapley values provide more or less the same explanation for both days, essentially only considering the more direct effect {\em temp}. The asymmetric conditional Shapley values, which are almost indistinguishable from the asymmetric causal Shapley values in this case, put a huge emphasis on the `root' cause {\em cosyear}. The (symmetric) causal Shapley values nicely balance the two extremes, giving credit to both season and temperature, to provide a sensible, but still different explanation for the two days.

\section{Discussion}

In real-world systems, understanding {\em why} things happen typically implies a causal perspective. It means distinguishing between important, contributing factors and irrelevant side effects. Similarly, understanding why a certain instance leads to a given output by a complex algorithm asks for those features that carry a significant amount of information contributing to the final outcome. Our insight was to recognize the need to properly account for the underlying causal structure between the features in order to derive meaningful and relevant attributive properties in the context of a complex algorithm.

For that, this paper introduced causal Shapley values, a model-agnostic approach to split a model's prediction of the target variable for an individual data point into contributions of the features that are used as input to the model, where each contribution aims to estimate the total effect of that feature on the target and can be decomposed into a direct and an indirect effect. We contrasted causal Shapley values with (interventional interpretations of) marginal and (asymmetric variants of) conditional Shapley values. We proposed a novel algorithm to compute these causal Shapley values, based on causal chain graphs. All that a practitioner needs to provide is a partial causal order (as for asymmetric Shapley values) and a way to interpret dependencies between features that are on an equal footing. Existing code for computing conditional Shapley values is easily generalized to causal Shapley values, without additional computational complexity. Computing conditional and causal Shapley values can be considerably more expensive than computing marginal Shapley values due to the need to sample from conditional instead of marginal distributions, even when integrated with computationally efficient approaches such as KernelSHAP~\cite{lundberg2017unified} and TreeExplainer~\cite{lundberg2020local}.

Our approach should be a promising step in providing clear and intuitive explanations for predictions made by a wide variety of complex algorithms, that fits well with natural human understanding and expectations. Additional user studies should confirm to what extent explanations provided by causal Shapley values align with the needs and requirements of practitioners in real-world settings. Similar ideas may also be applied to improve current approaches for (interactive) counterfactual explanations~\cite{wachter2017counterfactual} and properly distinguish between direct and total effects of features on a model's prediction. If successful, causal approaches that better match human intuition may help to build much needed trust in the decisions and recommendations of powerful modern-day algorithms. 


\section*{Broader Impact}

Our research, which aims to provide an explanation for complex machine learning models that can be understood by humans, falls within the scope of explainable AI (XAI). XAI methods like ours can help to open up the infamous ``black box'' of complicated machine learning models like deep neural networks and decision tree ensembles. A better understanding of the predictions generated by such models may provide higher trust~\cite{ribeiro2016should}, detect flaws and biases~\cite{kusner2017counterfactual}, higher accuracy~\cite{bhatt2020explainable}, and even address the legal ``right for an explanation'' as formulated in the GDPR~\cite{gdpr2017}.

Despite their good intentions, explanation methods do come with associated risks. Almost by definition, any sensible explanation of a complex machine learning system involves some simplification and hence must sacrifice some accuracy. It is important to better understand what these limitations are~\cite{kumar2020problems}. Model-agnostic general purpose explanation tools are often applied without properly understanding their limitations and over-trusted~\cite{kaur2020interpreting}. They could possibly even be misused just to check a mark in internal or external audits. Automated explanations can further give an unjust sense of transparency, sometimes referred to as the `transparency fallacy'~\cite{edwards2017slave}: overestimating one's actual understanding of the system. Last but not least, tools for explainable AI are still mostly used as an internal resource by engineers and developers to identify and reconcile errors~\cite{bhatt2020explainable}.

Causality is essential to understanding any process and system, including complex machine learning models. Humans have a strong tendency to reason about their environment and to frame explanations in causal terms~\cite{sloman2005causal,lombrozo2017causal} and causal-model theories fit well to how humans, for example, classify objects~\cite{rehder2003causal}. In that sense, explanation approaches like ours, that appeal to a human's capability for causal reasoning should represent a step in the right direction~\cite{mittelstadt2019explaining}.

\begin{ack}
This research has been partially financed by the Netherlands Organisation for Scientific Research (NWO), under project 617.001.451. 
\end{ack}


\section{\textit{Do}-calculus for cyclic graphs}

For completeness, we here repeat the rules of {\em do}-calculus for cyclic graphs, in the notation of the generalized ID algorithm of ~\cite{forre2019causal}, which generalizes~\cite{pearl2012calculus}. We are given a causal graph $G$. To each node $X_i$ which is intervened upon, we add an `intervention node' $I_{X_i}$, with a directed edge from $I_{X_i}$ to $X_i$ that we clamp to the value $x_i$. The corresponding graph is called $\hat{G}^+$. $\hat{G}_{\mathit{do}(\mathbf{W})}$ is now obtained by removing from $\hat{G}^+$ all incoming edges to variables that are part of $\mathbf{W}$, except those from the corresponding intervention nodes $I_{\mathbf{W}}$. We use shorthand
\[
\mathbf{Y} \stackrel[G]{\sigma}{\Perp} \mathbf{X} ~|~ \mathbf{Z}, \mathit{do}(\mathbf{W})
\]
to indicate that $\mathbf{Y}$ and $\mathbf{X}$ are $\sigma$-separated by $\mathbf{Z}$ in the graph $\hat{G}_{\mathit{do}(\mathbf{W})}$.  $\sigma$-separation is a generalization of standard d-separation (see~\cite{forre2019causal} for details).

\textit{Do}-calculus now consists of the following three inference rules that can be used to map interventional and observational distributions.
\begin{enumerate}
	\item Insertion/deletion of observation:
	\[
	 P(\mathbf{Y}|\mathbf{X},\mathbf{Z},\mathit{do}(\mathbf{W})) = P(\mathbf{Y}|\mathbf{Z},\mathit{do}(\mathbf{W})) \mbox{~~if~~} \mathbf{Y} \stackrel[G]{\sigma}{\Perp} \mathbf{X} ~|~ \mathbf{Z},\mathit{do}(\mathbf{W}) \: .
	\]
	\item Action/observation exchange:
	\[
	P(\mathbf{Y}|\mathit{do}(\mathbf{X}),\mathbf{Z},\mathit{do}(\mathbf{W})) = P(\mathbf{Y}|\mathbf{X},\mathbf{Z},\mathit{do}(\mathbf{W})) \mbox{~~if~~} \mathbf{Y} \stackrel[G]{\sigma}{\Perp} I_{\mathbf{X}} ~|~ \mathbf{X}, \mathbf{Z}, \mathit{do}(\mathbf{W})  \: .
	\]
	\item Insertion/deletion of actions:
	\[
	P(\mathbf{Y}|\mathit{do}(\mathbf{X}),\mathbf{Z},\mathit{do}(\mathbf{W})) = P(\mathbf{Y}|\mathbf{Z},\mathit{do}(\mathbf{W})) \mbox{~~if~~} \mathbf{Y} \stackrel[G]{\sigma}{\Perp} I_{\mathbf{X}} ~|~ \mathbf{Z},\mathit{do}(\mathbf{W}) \: .
	\]	
\end{enumerate}
Through consecutive application of these rules, we can try to turn any interventional probability of interest into an observational probability.

\section{Shapley values for linear models}

\begin{figure}
		\centering
		\tikzstyle{arrow} = [thick,->,>=stealth]
		\tikzstyle{dashedarrow} = [thick,->,>=stealth,dashed]
		\tikz{
		
			\node[latent] at (0,0) (y1) {$Y$};%
			\node[latent] at (0,1.5) (x12) {$X_2$};
			\node[latent] at (0,3) (x11) {$X_1$};
			\node[text height=1em, anchor=north] at (0,4.2) {Chain and fork};
			\draw[arrow] (x12) -- node[anchor=west]{$\beta_2$} (y1);
			\draw[arrow] (x11) -- (x12);
			\draw[arrow] (x11) to[bend right] node[anchor=east]{$\beta_1$}(y1);
		
			\node[latent] at (3,0) (y3) {$Y$};
			\node[latent] at (3,1.5) (x32) {$X_2$};
			\node[latent] at (3,3) (x31) {$X_1$};
			\node[latent] at (3.6,2.25) (z)  {$Z$};
			\node[text height=1em, anchor=north] at (3,4.2) {Confounder};
			\draw[arrow] (x32) -- node[anchor=west]{$\beta_2$} (y3);
			\draw[arrow] (z) -- (x31);
			\draw[arrow] (z) -- (x32);
			\draw[arrow] (x31) to[bend right] node[anchor=east]{$\beta_1$}(y3);
		
			\node[latent] at (6,0) (y4) {$Y$};%
			\node[latent] at (6,1.5) (x42) {$X_2$};
			\node[latent] at (6,3) (x41) {$X_1$};
			\node[text height=1em, anchor=north] at (6,4.2) {Cycle};
			\draw[arrow] (x42) -- node[anchor=west]{$\beta_2$} (y4);
			\draw[arrow] (x41) to[bend left=15] (x42);
			\draw[arrow] (x42) to[bend left=15] (x41);
			\draw[arrow] (x41) to[bend right] node[anchor=east]{$\beta_1$}(y4);			
		}
	\caption{Three causal models with the same observational distribution over features, yet a different causal structure. To connect to the models in the main text, we set $\beta_1 = 0$ and $\beta_2 = \beta$, except that for the `fork' we set $\beta_2 = 0$, $\beta_1 = \beta$, and then swap the indices.}
	\label{fig:threemodels}
\end{figure}

We will make use of the {\em do}-calculus rules above to derive the causal Shapley values for the four different models in Figure~1 in the main text. To this end, we consider the three models in Figure~\ref{fig:threemodels} that predict $f(x_1,x_2) = \beta_1 x_1 + \beta_2 x_2$ for general values of $\beta_1$ and $\beta_2$. All three models have the same observational probability distribution, with $\mathbb{E}[X_i] = \bar{x}_i$ and $\mathbb{E}[X_{3-i}|X_i=x_i] = \alpha_i x_i$, for $i = 1,2$, yet different causal structures. We will arrive at the main text's results for the `chain', `confounder', and `cycle' by setting $\beta_1 = 0$, whereas for the `fork' we set $\beta_2 = 0$ and swap the two indices. We then further need to take $\bar{x}_1 = \bar{x}_2 = 0$, and $\alpha = \alpha_2$.

Following the definitions in the main text, the contribution of feature $i$ given permutation $\pi$ is the difference in value function before and after setting the feature to its value:
\begin{equation}
{\phi}_i(\pi) = {v}(\{j: j \preceq_\pi i\}) - {v}(\{j: j \prec_\pi i\}) \: ,
\label{eq:contperm}
\end{equation}
with value function
\begin{equation}
{v}(S) = \mathbb{E} \left[f(\mathbf{X}) | \lvdo{S} \right] = \int d\mathbf{X}_{\bar{S}} \: P(\mathbf{X}_{\bar{S}}|\svdo{S}) f(\mathbf{X}_{\bar{S}},\mathbf{x}_S) \: ,
\label{eq:valuedef}
\end{equation}
where we use shorthand $\svdo{}$ for $\lvdo{}$. Combining these two definitions and substituting $f(\mathbf{x}) = \sum_i \beta_i x_i$, we obtain
\[
{\phi}_i(\pi) =
\beta_i \left(x_i - \mathbb{E} [X_i | \svdo{j: j \prec_\pi i}]\right) + \sum_{k \succ_\pi i} \beta_k \left( \mathbb{E} [X_k | \svdo{j: j \preceq_\pi i}] - \mathbb{E} [X_k | \svdo{j: j \prec_\pi i}] \right) \: .
\]
The first term corresponds to the direct effect, the second one to the indirect effect. Symmetric causal Shapley values will follow by averaging these contributions for the two possible permutations $\pi = (1,2)$ and $\pi = (2,1)$. Conditional Shapley values result when replacing conditioning by intervention with conventional conditioning by observation, marginal Shapley values by not conditioning at all.

To start with the latter, we immediately see that for {\em marginal Shapley values} the indirect effect vanishes and the direct effect simplifies to
\[
\phi_i = \phi_i(\pi) = \beta_i (x_i - \mathbb{E}[X_i])  = \beta_i (x_i - \bar{x}_i) \: ,
\]
as also derived in~\cite{aas2019explaining}.

For symmetric conditional Shapley values, we do get different contributions for the two different permutations, but by definition still the same for the three different models:
\begin{eqnarray}
\phi_1(1,2) \hspace*{-2.5mm} & = & \hspace*{-2.5mm} \beta_1 (x_1 - \mathbb{E}[X_1]) + \beta_2 (\mathbb{E}[X_2|x_1] - \mathbb{E}[X_2]) = \beta_1 (x_1 - \bar{x}_1) + \beta_2 \alpha_1 (x_1 - \bar{x}_1) \nonumber \\
\phi_2(1,2) \hspace*{-2.5mm} & = & \hspace*{-2.5mm} \beta_2 (x_2 - \mathbb{E}[X_2|x_1]) = \beta_2 (x_2 - \bar{x}_2) - \beta_2 \alpha_1 (x_1 - \bar{x}_1) \: .
\label{eq:asymmetric}
\end{eqnarray}
Here the first term in the contribution for the first feature corresponds to the direct effect and the second term to the indirect effect. The contribution for the second feature only consists of a direct effect. The contributions for the other permutation follow by swapping the indices and the final Shapley values by averaging to arrive at the {\em symmetric conditional Shapley values}
\begin{eqnarray}
\phi_1 \hspace*{-2.5mm} & = & \hspace*{-2.5mm} \beta_1 (x_1 - \bar{x}_1) - {1 \over 2} \beta_1 \alpha_2 (x_2 - \bar{x}_2) + {1 \over 2} \beta_2 \alpha_1 (x_1 - \bar{x}_1) \nonumber \\
\phi_2 \hspace*{-2.5mm} & = & \hspace*{-2.5mm} \beta_2 (x_2 - \bar{x}_2) - {1 \over 2} \beta_2 \alpha_1 (x_1 - \bar{x}_1) + {1 \over 2} \beta_1 \alpha_2 (x_2 - \bar{x}_2) \: ,
\label{eq:symmetric}
\end{eqnarray}
where now the first two terms constitute the direct effect and the third term the indirect effect.

\begin{table}
\begin{center}
\begin{tabular}{r|ccc} \toprule
expectation & chain & confounder & cycle \\ \midrule
$\mathbb{E}[X_1|\sdo{2}]$ & $\mathbb{E}[X_1]$ & $\mathbb{E}[X_1]$ & $\mathbb{E}[X_1|x_2]$ \\
$\mathbb{E}[X_2|\sdo{1}]$ & $\mathbb{E}[X_2|x_1]$ & $\mathbb{E}[X_2]$ & $\mathbb{E}[X_2|x_1]$ \\
\bottomrule
\end{tabular}
\end{center}
\caption{Turning expectations under conditioning by intervention into expectations under conventional conditioning by observation for the three models in Figure~\ref{fig:threemodels}.}
\label{tab:rewriting}
\end{table}

The {\em asymmetric conditional Shapley values} consider both permutations for the confounder and the cycle, and hence are equivalent to the symmetric Shapley values for those models. Yet for the chain, they only consider the permutation $\pi(1,2)$ and thus yield $\bm{\phi} = \bm{\phi}(1,2)$ from~(\ref{eq:asymmetric}).

To go from the symmetric conditional Shapley values to the causal symmetric Shapley values, we follow the same line of reasoning, but have to replace $\mathbb{E}[X_2|x_1]$ by $\mathbb{E}[X_2|\sdo{1}]$ and $\mathbb{E}[X_1|x_2]$ by $\mathbb{E}[X_1|\sdo{2}]$. Table~\ref{tab:rewriting} tells whether
the expectations under conditioning by intervention reduce to expectations under conditioning by observation (because of the second rule of \textit{do}-calculus above) or to marginal expectations (because of the third rule). For the chain we have
\[
P(X_2 | \sdo{1}) = P(X_2 | x_1) \mbox{~~since~~} X_2 \stackrel[G]{\sigma}{\Perp} I_{X_1} ~|~ X_1 \mbox{~~(rule 2),~yet~~} P(X_1 | \sdo{2}) = P(X_1) \mbox{~~since~~} X_1 \stackrel[G]{\sigma}{\Perp} I_{X_2} \mbox{~~(rule 3),}
\]
for the confounder
\[
P(X_2 | \sdo{1}) = P(X_2) \mbox{~~since~~} X_2 \stackrel[G]{\sigma}{\Perp} I_{X_1}  \mbox{~~and~~} P(X_1 | \sdo{2}) = P(X_1) \mbox{~~since~~} X_1 \stackrel[G]{\sigma}{\Perp} I_{X_2} \mbox{~~(rule 3),}
\]
and for the cycle
\[
P(X_2 | \sdo{1}) = P(X_2 | x_1) \mbox{~~since~~} X_2 \stackrel[G]{\sigma}{\Perp} I_{X_1} ~|~ X_1 \mbox{~~and~~} P(X_1 | \sdo{2}) = P(X_1|x_2) \mbox{~~since~~} X_1 \stackrel[G]{\sigma}{\Perp} I_{X_2} ~|~ X_2 \mbox{~~(rule 2).}
\]
Consequently, for the confounder the {\em symmetric} and {\em asymmetric causal Shapley values} coincide with the marginal Shapley values (consistent with~\cite{janzing2019feature}) and for the cycle with the symmetric conditional Shapley values from~(\ref{eq:symmetric}). For the chain, the causal symmetric Shapley values become
\begin{eqnarray}
\phi_1(1,2) \hspace*{-2.5mm} & = & \hspace*{-2.5mm} \beta_1 (x_1 - \bar{x}_1) + {1 \over 2} \beta_2 \alpha_1 (x_1 - \bar{x}_1) \nonumber \\
\phi_2(1,2) \hspace*{-2.5mm} & = & \hspace*{-2.5mm} \beta_2 (x_2 -\bar{x}_2) - {1 \over 2} \beta_2 \alpha_1 (x_1 - \bar{x}_1) \: ,
\label{eq:asymmetric2}
\end{eqnarray}
where the asymmetric causal Shapley values coincides with the asymmetric conditional Shapley values from~(\ref{eq:asymmetric2}).

Collecting all results and setting $\bar{x}_1 = \bar{x}_2 = \beta_1 = 0$, $\beta_2 = \beta$, and $\alpha_1 = \alpha$ (after swapping the indices for the `fork'), we arrive at the Shapley values reported in Figure~1 in the main text. Note that for most Shapley values, the indirect effect for the second feature vanishes because we chose to set $\beta_1 = 0$. The exceptions, apart from the marginal Shapley values, are the causal Shapley values for the chain and the confounder, as well as the asymmetric conditional Shapley values for the chain: these show no indirect effect for the second feature even for nonzero $\beta_1$.

\section{Proofs and corollaries on causal chain graphs}

In this section we expand on the proof of Theorem~1 in the main text and add some corollaries to link back to other approaches for computing Shapley values.

The probability distribution for a causal chain graph boils down to a directed acyclic graph of chain components:
\begin{equation}
P(\mathbf{X}) = \prod_{\tau \in {\cal T}} P(\mathbf{X}_\tau|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)}) \: .
\label{eq:chaingraph}
\end{equation}
For each (fully connected) chain component, we further need to specify whether (surplus) dependencies within the component are due to confounding or due to mutual interactions. Given this information, we can turn any causal query into an observational distribution with the following interventional formula.

\begin{theorem}
	For causal chain graphs, we have the interventional formula 
	\begin{eqnarray}
	P(\mathbf{X}_{\bar{S}}|\lvdo{S}) \hspace*{-2.5mm} & = & \hspace*{-2.5mm} \prod_{\tau \in {\cal T}_{\textrm{\upshape \scriptsize confounding}}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S}) \times \nonumber \\
	&& \prod_{\tau \in {\cal T}_{\overline{\textrm{\upshape \scriptsize confounding}}}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\mathbf{x}_{\tau \cap S}) \: .
	\label{eq:chaininterventional}
	\end{eqnarray}
\end{theorem}

\begin{proof}
Plugging in~(\ref{eq:chaingraph}) and using shorthand $\svdo{} = \lvdo{}$, we obtain
\[
P(\mathbf{X}_{\bar{S}}|\svdo{S}) = P(\mathbf{X}|\svdo{S}) = \prod_{\tau \in {\cal T}} P(\mathbf{X}_{\tau}|\mathbf{X}_{\tau' \prec_G \tau},\svdo{S})
\isequaldono{1} \prod_{\tau \in {\cal T}} P(\mathbf{X}_\tau|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)},\svdo{S}) = \prod_{\tau \in {\cal T}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\svdo{S}) \: ,
\]
where in the second step we made use of {\em do}-calculus rule~(1): the conditional independencies in the causal chain graph $G$ are preserved when we intervene on some of the variables. Now rule~(3) tells us that we can ignore any interventions from nodes in components further down the causal chain graph as well as those from higher up that are shielded by the direct parents:
\[
P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\svdo{S})
\isequaldono{3} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\svdo{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\svdo{\tau \cap S}) \: .
\]
Rule~(2) then states that conditioning by intervention upon variables higher up in the causal chain graph is equivalent to conditioning by observation:
\[
P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\svdo{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\svdo{\tau \cap S}) \isequaldono{2}
P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\svdo{\tau \cap S}) \: .
\]
For a chain component with dependencies induced by a common confounder, rule (3) applies once more and makes that we can ignore the interventions:
\[
P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\svdo{\tau \cap S}) = P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S}) \: .
\]
For a chain component with dependencies induced by mutual interactions, rule (2) again applies and allows us to replace conditioning by intervention with conditioning by observation:
\[
P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\svdo{\tau \cap S}) = P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{ \mathop{\textit{\scriptsize pa}}(\tau) \cap S},\mathbf{x}_{\tau \cap S})) \: .
\]
\end{proof}

Algorithm~\ref{alg:sampling} provides pseudocode on how to estimate the value function ${v}(S)$ by drawing samples from the interventional probability~(\ref{eq:chaininterventional}). It assumes that a user has specified a partial causal ordering of the features, which is translated to a chain graph with components ${\cal T}$, and for each (non-singleton) component $\tau$ whether or not surplus dependencies are the result of confounding. Other prerequisites include access to the model $f()$, the feature vector $\mathbf{x}$, (a procedure to sample from) the observational probability distribution $P(\mathbf{X})$, and the number of samples $\onder{N}{samples}$.

\begin{algorithm}[t]
	\caption{Compute the value function ${v}(S)$ under conditioning by intervention.}
	\label{alg:sampling}
	\begin{algorithmic}[1]
		\Function{ValueFunction}{$S$}
		\For{$k \gets 1 \mbox{~to~} \onder{N}{samples}$}
		\ForAll{$j \gets 1 \mbox{~to~} |{\cal T}|$} \Comment{run over all chain components in their causal order}
			\If{confounding$(\tau_j)$} 
				\ForAll{$i \in \tau_j \cap \bar{S}$}
					\State Sample $\tilde{x}_i^{(k)} \sim P(X_i|\tilde{\mathbf{x}}_{\mathop{\textit{\scriptsize pa}}(\tau_j) \cap \bar{S}}^{(k)},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau_j) \cap \bar{S}})$ \Comment{can be drawn independently}
				\EndFor
			\Else
				\State Sample $\tilde{\mathbf{x}}_{\tau_j \cap \bar{S}}^{(k)} \sim P(\mathbf{X}_{\tau_j \cap \bar{S}}|\tilde{\mathbf{x}}_{\mathop{\textit{\scriptsize pa}}(\tau_j) \cap \bar{S}}^{(k)},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau_j) \cap \bar{S}},\mathbf{x}_{\tau_j \cap S})$ \Comment{e.g., Gibbs sampling}
			\EndIf
		\EndFor
		\EndFor
		\State ${v} \gets {\displaystyle {1 \over \onder{N}{samples}} \sum_{k=1}^{\onder{N}{samples}} f(\mathbf{x}_S,\tilde{\mathbf{x}}_{\bar{S}}^{(k)})}$
		\State \Return ${v}$
		\EndFunction
	\end{algorithmic}
\end{algorithm}

Theorem~1 connects to observations made and algorithms proposed in recent papers.
\begin{corollary}
	With all features combined in a single component and all dependencies induced by confounding, as in~\cite{janzing2019feature}, causal Shapley values are equivalent to marginal Shapley values. 
\end{corollary}
\begin{proof}
	With just a single confounded component $\tau$, $\mathop{\textit{pa}}(\tau) = \emptyset$ and~(\ref{eq:chaininterventional}) reduces to $P(\mathbf{X}_{\bar{S}})$.
\end{proof}
\begin{corollary}
	With all features combined in a single component and all dependencies induced by mutual interactions, causal Shapley values are equivalent to conditional Shapley values as proposed in~\cite{aas2019explaining}.
\end{corollary}
\begin{proof}
	With just a single non-confounded component $\tau$, $\mathop{\textit{pa}}(\tau) = \emptyset$ and~(\ref{eq:chaininterventional}) reduces to $P(\mathbf{X}_{\bar{S}}|\mathbf{x}_{S})$.
\end{proof}
\begin{corollary}
	When we only take into account permutations that match the causal ordering and assume that all dependencies within chain components are induced by mutual interactions, the resulting asymmetric causal Shapley values are equivalent to the asymmetric conditional Shapley values as defined in~\cite{frye2019asymmetric}.
\end{corollary}
\begin{proof}
	Following~\cite{frye2019asymmetric}, asymmetric Shapley values only include those permutations $\pi$ for which $i \prec_\pi j$ for all known ancestors $i$ of descendants $j$ in the causal graph. For those permutations, we are guaranteed to have $\tau \prec_G \tau'$ for all $\tau, \tau' \in \mathcal{T}$ such that $\tau \cap S \neq \emptyset, \tau' \cap \bar{S} \neq \emptyset$. That is, the chain components that contain features from $S$ are never later in the causal ordering of the chain graph $G$ than those that contain features from $\bar{S}$. We then have
	\[
	P(\mathbf{X}_{\bar{S}}|\mathbf{x}_S) = \prod_{\tau \in {\cal T}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{S}) = \prod_{\tau \in {\cal T}} P(\mathbf{X}_{\tau \cap \bar{S}}|\mathbf{X}_{\mathop{\textit{\scriptsize pa}}(\tau)  \cap \bar{S}},\mathbf{x}_{\mathop{\textit{\scriptsize pa}}(\tau) \cap S},\mathbf{x}_{\tau \cap S}) = P(\mathbf{X}_{\bar{S}}|\svdo{S}) \: ,
	\]
	where in the last step we used interventional formula~(\ref{eq:chaininterventional}) in combination with the fact that $\onder{{\cal T}}{confounding} = \emptyset$.
\end{proof}

\section{Additional illustrations on the bike rental data}

Figure~\ref{fig:sinaplots} shows sina plots for asymmetric conditional Shapley values (left) and marginal Shapley values (right), to be compared with the sina plots for symmetric causal Shapley values in Figure~3 of the main text. In this case, the sina plots for asymmetric causal Shapley values are virtually indistinguishable from those for the asymmetric conditional Shapley values.

It can be seen that the marginal Shapley values strongly focus on temperature, largely ignoring the seasonal variables. The asymmetric Shapley values do the opposite: they focus on the seasonal variables, in particular {\em cosyear} and put much less emphasis on the temperature variables.

\begin{figure}[t]
	\centering
	\begin{minipage}{.49\linewidth}
	\includegraphics[width=\textwidth]{sina_plot_asymmetric.pdf}
	\end{minipage}
	\begin{minipage}{.49\linewidth}
	\includegraphics[width=\textwidth]{sina_plot_marginal.pdf}
	\end{minipage}
	\caption{Sina plots of asymmetric (conditional) Shapley values (left) and marginal Shapley values (right) on the bike rental data. See Figure~3 in the main text for further details.}
	\label{fig:sinaplots}
\end{figure}

\section{Comparing symmetric and asymmetric Shapley values on the XOR problem}

We consider the standard XOR problem with binary features $X_1$ and $X_2$ and binary output $Y$:
\begin{center}
\begin{tabular}{cc|c}
	$X_1$ & $X_2$ & Y \\ \hline
	0     & 0     & 0 \\
	0     & 1     & 1 \\
	1     & 0     & 1 \\
	1     & 1     & 0
\end{tabular}
\end{center}
We generate a dataset of $n$ samples by drawing features and corresponding outputs with probabilities $p_{ij} = P(X_1=i,X_2=j)$. We will choose $p_{00} = p_{11} = \frac{1}{4} (1 + \epsilon)$ and $p_{01} = p_{10} = \frac{1}{4} (1 - \epsilon)$. With $\epsilon > 0$, the probability of the two features having the same values is larger than the probability of them having different values. $\hat{p}_{ij}$ is the same probability estimated from the data, e.g., by computing the frequencies of the four input combinations. We train a neural network on the generated data, which yields a function $\hat{f}(X_1,X_2)$ hopefully closely approximating the correct XOR function. The parameter $\epsilon$ captures the dependency between the two features and can be interpreted as a measure of the causal strength when the two features are causally related.

We will now compute the various Shapley values for the data point $(X_1,X_2) = (0,0)$. The value functions with all features either `out-of-coalition' or `in-coalition' are the same for all Shapley values:
\begin{eqnarray*}
\nu(\{\}) & = & \mathbb{E}\left[f(\mathbf{X})\right] = \sum_{i,j} \hat{p}_{ij} \hat{f}(i,j) \approx \frac{1}{2} (1 - \epsilon) \\
\nu(\{1,2\}) & = & \hat{f}(0,0) \approx 0 \: ,
\end{eqnarray*}
where we use the convention that the Shapley values computed from the fitted probabilities and learned neural network appear before the $\approx$-sign, and those that we obtain when the fitted probabilities equal the probabilities used to generate the data and when the learned neural network equals the XOR function after the $\approx$-sign.

The value functions for the case that one input is `in-coalition' and the other `out-of-coalition' does depend on the type of Shapley value under consideration. For the marginal Shapley values we get
\begin{eqnarray}
\nu(\{1\}) & = & \mathbb{E}\left[f(0,X_2)\right] = \sum_{j} \left(\sum_i \hat{p}_{ij}\right) \hat{f}(0,j) \nonumber \approx {1 \over 2} \\
\nu(\{2\}) & = & \mathbb{E}\left[f(X_1,0)\right] = \sum_{i} \left(\sum_j \hat{p}_{ij}\right) \hat{f}(i,0) \approx \frac{1}{2} \: ,
\label{eq:marginal}
\end{eqnarray}
yet for the conditional Shapley values
\begin{eqnarray}
\nu(\{1\}) & = & \mathbb{E}\left[f(0,X_2)|X_1=0\right] = \sum_{j} \frac{\hat{p}_{0j}}{\sum_{i} \hat{p}_{ij}} \hat{f}(0,j) \approx \frac{1}{2}(1-\epsilon) \nonumber \\
\nu(\{2\}) & = & \mathbb{E}\left[f(X_1,0)|X_2=0\right] = \sum_{i} \frac{\hat{p}_{i0}}{\sum_{j} \hat{p}_{ij}} \hat{f}(i,0) \approx \frac{1}{2}(1-\epsilon) \: .
\label{eq:conditional}
\end{eqnarray}
The value functions for the causal Shapley values depend on the presumed causal model that generates the dependencies. In case the dependencies are assumed to be the result of confounding, we get the value functions in~(\ref{eq:marginal}) as for the marginal Shapley values and when the dependencies are assumed to be the result of mutual interaction the value functions in~(\ref{eq:conditional}) as for the conditional Shapley values. The more interesting case is when we assume a causal chain, e.g., $X_1 \rightarrow X_2$:
\begin{eqnarray}
\nu(\{1\}) & = & \mathbb{E}\left[f(0,X_2)|\mathit{do}(X_1=0)\right] = \mathbb{E}\left[f(0,X_2)|X_1=0\right] = \sum_{j} \frac{\hat{p}_{0j}}{\sum_{i} \hat{p}_{ij}} \hat{f}(0,j) \approx \frac{1}{2}(1-\epsilon) \nonumber \\
\nu(\{2\}) & = & \mathbb{E}\left[f(X_1,0)|\mathit{do}(X_2=0)\right] =\mathbb{E}\left[f(X_1,0)\right] = \sum_{i} \left(\sum_j \hat{p}_{ij}\right) \hat{f}(i,0) \approx \frac{1}{2} \: ,
\label{eq:causal}
\end{eqnarray}
and the same with indices 1 and 2 interchanged for the causal chain $X_2 \rightarrow X_1$.


\begin{figure}[t]
	\centering
	\includegraphics[width=\textwidth]{XOR_plot.pdf}
	\caption{The conditional symmetric, causal symmetric and (causal/conditional) asymmetric Shapley values of data point $(X_1,X_2) = (0,0)$ under the assumption of a causal chain $X_1 \rightarrow X_2$ for different causal strengths $\epsilon$. The bars indicate the mean Shapley value and standard deviation of 100 runs with a neural network trained on 100 data points generated according to a particular causal strength $\epsilon$. Red lines give the theoretical Shapley values, to be expected for an infinite amount of samples and when the neural networks would have perfectly learned to represent the XOR function. The third row gives the difference between the Shapley values for $X_1$ and $X_2$ and clearly shows the discontinuity of asymmetric Shapley values for $\epsilon = 0$.
	}
	\label{fig:xor_plot}
\end{figure}


Given these value functions, we can now compute the various Shapley values. For marginal and symmetric Shapley values we have
\begin{eqnarray}
\phi_1 & = & \frac{1}{2}[\nu(\{1\}) - \nu(\{\})] + \frac{1}{2}[\nu(\{1,2\}) - \nu(\{2\})] \nonumber \\
\phi_2 & = & \frac{1}{2}[\nu(\{2\}) - \nu(\{\})] + \frac{1}{2}[\nu(\{1,2\}) - \nu(\{1\}]) \: , \nonumber
\end{eqnarray}
whereas for asymmetric Shapley values, assuming the causal chain $X_1 \rightarrow X_2$,
\begin{eqnarray}
\phi_1 & = & \nu(\{1\}) - \nu(\{\}) \nonumber \\
\phi_2 & = & \nu(\{1,2\}) - \nu(\{1\}) \: , \nonumber
\end{eqnarray}
and the same with indices 1 and 2 interchanged for the causal chain $X_2 \rightarrow X_1$.

With the expressions above, we can compute the various Shapley values based on a learned neural network and the actual frequencies of the generated feature combinations and compare those with the theoretical values obtained when the estimated frequencies equal the probabilities used to generate the data and the neural network indeed managed to learn the XOR function. For the latter we distinguish the following cases.
\begin{description}
	\item[identical:] $\phi_1 = \phi_2 \approx \frac{1}{4} \epsilon - \frac{1}{4}$. This applies to marginal, symmetric conditional, symmetric causal assuming confounding, symmetric causal assuming mutual interaction.
	\item[symmetric causal:] $\phi_1 \approx -\frac{1}{4}$ and $\phi_2 \approx \frac{1}{2}\epsilon - \frac{1}{4}$ assuming the causal chain $X_1 \rightarrow X_2$ and vice versa for $X_1 \rightarrow X_2$.
	\item[asymmetric:] $\phi_1 \approx 0$ and $\phi_2 \approx \frac{1}{2}\epsilon - \frac{1}{2}$ assuming the causal chain $X_1 \rightarrow X_2$ and vice versa for $X_1 \rightarrow X_2$. These apply both to asymmetric conditional and asymmetric causal.
\end{description}

In this example, symmetric causal Shapley values are clearly to be preferred over asymmetric causal Shapley values for small causal strengths. Inserting a causal link with zero strength ($\epsilon = 0$), asymmetric Shapley values jump from the symmetric $\phi_1 = \phi_2 \approx -\frac{1}{4}$ to the completely asymmetric $\phi_1 \approx 0$ and $\phi_2 \approx - \frac{1}{2}$, assigning all credit to the second feature, even though the first feature in reality does not affect the second feature at all. Symmetric Shapley values, on the other hand, are insensitive to the insertion of a causal link with zero strength: in the limit $\epsilon \rightarrow 0$ symmetric causal Shapley values correctly converge to marginal Shapley values.

Figure~\ref{fig:xor_plot} shows the results of a series of simulations, computing different Shapley values for trained neural networks and comparing these to the theoretical values. The discontinuity of asymmetric Shapley values (conditional and causal asymmetric Shapley values are identical in this example) is most clearly seen in the third row, showing the difference between the Shapley values for $X_1$ and $X_2$. Symmetric conditional Shapley values do not distinguish between the Shapley values for $X_1$ and $X_2$ for any causal strength $\epsilon$, whereas the symmetric causal Shapley values are identical for $\epsilon = 0$ and then slowly start to deviate for larger values of $\epsilon$.

\section{Shapley values for predicting dementia}

We consider the Alzheimer's disease data set obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (\url{http://adni.loni.usc.edu}). The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment and early Alzheimer's disease. As features we consider age ({\em age}), gender ({\em gender}), education level ({\em pteducat}), fudeoxyglucose ({\em FDG}), amyloid beta ({\em ABETA}), phosphorylated
tau ({\em PTAU}), and the number of apolipoprotein alleles ({\em APOE4}). Data was normalised and randomly split in 80\% training and 20\% test set. We consider a binary classification problem, where we grouped together patients with mild cognitive impairment and early Alzheimer's disease, to distinguish these from the healthy ``cognitive normal'' subjects. We trained a multi-layered perceptron with five hidden units.

To compute the causal Shapley values, we chose the partial order (\{{\em gender}, {\em APOE4}, {\em age}, {\em pteducat}\}, \{{\em ABETA}\}, \{{\em FDG}, {\em PTAU}\}), in line with the ``gold standard'' causal graph from~\cite{shen2020challenges}. We assume confounding in the first and third component. Since interventional expectations over variables in the first component simplify to marginal expectations, we can sample the discrete variables {\em gender} and {\em APOE4} from their empirical distributions. All other variables are sampled from conditional Gaussian distributions.

\begin{figure}[t]
	\includegraphics[width=\textwidth]{sina_plot_adni.pdf}
	\caption{Sina plots of marginal (left), symmetric causal (middle), and asymmetric causal Shapley values (right) for a multi-layered perceptron trained on the ADNI dataset. See the text for further explanations.}
	\label{fig:sina_adni}
\end{figure}

The sina plots in Figure~\ref{fig:sina_adni} show the marginal, symmetric causal, and asymmetric causal Shapley values, respectively, for the predicted probability of having dementia or mild cognitive impairment. As the dependencies among the features are relatively weak, the marginal and symmetric causal Shapley values are quite similar. The asymmetric causal Shapley values for different values of {\em APOE4}, a known risk factor of {\em ABETA}, are more clearly separated than those for the marginal and symmetric causal Shapley values. The asymmetric Shapley values provide an explanation with a relatively strong focus on the root cause {\em APOE4}, which first gets all the credit for the indirect effect through {\em ABETA}.

\begin{ack}
ADNI data collection and sharing was funded by the Alzheimer's Disease
Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). See \url{http://adni.loni.usc.edu/wp-content/themes/freshnews-dev-v2/documents/policy/ADNI_Data_Use_Agreement.pdf}, section 12, for further contributions and \url{http://adni.loni.usc.edu/wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf} for a complete listing of ADNI investigators.
\end{ack}

