Improving Local Explainability By Learning Causal Graphs From Data
Authors: Daan Roos
Sebastian Gerwinn
Jan-Willem van de Meent
Sara Magliacane
Venue: Accepted by TMLR ()

Abstract: Causal Shapley values take into account causal relations among dependent features to adjust the contributions of each feature to a prediction. A limitation of this approach is that it can only leverage known causal relations. In this work we combine the computation of causal Shapley values with causal discovery, i.e., learning causal graphs from data. In particular, we compute causal explanations across the Markov Equivalence Class (MEC), a set of candidate causal graphs learned from observational data, providing a list of causal Shapley values that explain the prediction. We propose two methods for estimating this list efficiently, drawing on the equivalences of the interventional distributions for a subset of the causal graphs. We evaluate our methods on synthetic and real-world data, showing that they provide explanations that are more consistent with the true causal effects compared to traditional Shapley value approaches that disregard causal relations. Our results show that even when the Markov Equivalence Class is learned incorrectly, in most settings the explanations of our framework are on average closer to true causal Shapley values than marginal and conditional Shapley values.

Published in Transactions on Machine Learning Research (03/2026)
Improving Local Explainability By Learning Causal Graphs From Data
Daan Roos MnA-Bosch Delta Lab, University of Amsterdam
d.f.c.roos@wo..l
Scbastian Gerwinn
Sebastian, Gerwinn Ode. bosch. com
Bosch Center for Artificial Intellience, Robert Bosch GmbH
Jan- Willem de Meent
j. w...onddemeennt@wv...l
YuA-Bosch Delta Lab, University of Amsterdam
Sara Magliacane
son.magliaonee@gaa. com
Cousal Machine Learning group, Saarland University, &
Amsterdam Machine Learning Lab, Universily of Amsterdam
Reviewed on OpenReview: https: // openreview. hormfo id= A1bX7TRQQQ
Abstract
Causal Shapley values take into account causal relations among dependent features to adjust the contributions of each feature to a prediction. A limitation of this approach is that can only leverage known causal relations. In this work we combine the computation of causal Shapley values with causal discovery, i.c., learning causal graphs from data. In particular, we compute causal explanations across the Markov Equivalence Class (MEC), set of candidate causal graphs learned from observational data, providing list of causal Shapley values that explain the prediction. We propose methods for estimating this list efficiently, drawing on the equivalences of the interventional distributions for subset of the causal graphs. We evaluate our methods on synthetic and real-world data, showing that they provide explanations that are more consistent with the true causal effects compared to traditional Shapley value approaches that disregard causal relations. Our results show that even when the Markov Equivalence Class is learned incorrectly, in most settings the explanations of our framework marginal and conditional Shapley values.
average closer to true causal Shapley values than
1 Introduction
The advancement of machine learning models has ushered in a new era of AT capabilities, surpassing simpler models in performance. However, this increase in complexity brings significant challenge: a marked docrease in interpretability with advanced models often resembling a black box (Zhang et al., 2020; Carvalho al., 2019). This opacity is not just technical concern, but it becomes matter of ethical and legal importance in high-stakes applications, e.g,, medical diagnostics autonomous vehicles, where the inability to explain a model's reasoning can lead to issues of trust and accountability (Doshi-Vekle Velez & Kim, 2017).
Shapley values are a notable method used to explain how machine learning models make their decisions (Zhou et al., 2021). They offer principled local approach to attributing the output of model its input features, breaking down the prediction into contributions from each feature. When using Shapley values for explainability, various methods assume independence between the model's input features (Strumbelj &

---

Published in Transactions on Machine Learning Research (03/2026)
Input datapoint x
Prediction model
For each subset features S:
For each group Gy estimate vS)::
G3
Group graphs interventional
Sampling method
@2
Causal Discovery
distributions for do(Xg)
Importance weighting
method
Data
Candidate causal graphs
List of groups of causal graphs
List of Causal Shapley values
Figure 1: Overview of our MEC Shapley framework that provides a range of causal Shapley values, describing the causal impact a feature on an input datapoint, even when the causal graph is unknown. We first apply causal discovery to learn a set of candidate causal graphs directly from the data. For cach subset of features S that we consider, we group these causal graphs based on sharing the same interventional distribution for an intervention on Xs in list of groups Gg. For each group of graphs, we can then estimate the value function v(S) for the input datapoint x" given the prodiction model with one of our two methods, sampling methods or importance weighting method, which allows to provide list of causal Shapley values that represent range of explanations ions for the inpul datapoint that reflect the different candidate causal relations between the features.
Kononenko, 2014; Lundberg & Lee, 2017). As a result, importance is only attributed to features that directly contribute to a prediction (Heskes et al., 2020). As a consequence, many inputs to the model will be off the data manifold (Frye et al., 2021) and may lead misleading explanations (Aas al., 2021). Relaxing this assumption lcads to explanations where correlated features share importance, regardless of their causal effect the prediction (Aas et al., 2021). As result, conditional Shapley values may altribute importance to features that are merely correlated with the prediction.
Heskes et al. (2020) integrate causal modeling with the Shapley value framework, leveraging the knowledge of causal relations among the features to adjust the contribution each feature. This ensures that features only contribute to explanations when they have a causal effect on the prediction, whether direct or indirect. This paper aims to address one of the key challenges in applying causal Shapley values: leveraging only known causal relations. As described in the overview in Fig. 1, we first use causal discovery methods, e.g. PC (Spirtes et al., 2000) and FGES (Chickering, 2003; Ramsey et al., 2017), to learn an equivalence class of causal graphs directly from data and then combine it with the causal Shaplcy value estimation in an efficient way. This process provides a set of explanations that reflect different possible causal structures, thereby representing causal uncertainty. Our contributions are as follows:
We introduce a framework for combining causal discovery with causal Shapley values (Heskes et al., 2020), extending their applicability beyond known causal relations. In particular, we address the challenge that causal discovery methods return set of candidate causal graphs and hence output is list of causal Shapley values.
We propose two methods to efficiently estimate causal Shapley values using only observational data, a sampling method based on simulating interventional data and an importance weighting method. Both of these methods leverage the ID algorithm (Shpitser & Pearl, 2006) to identify the estimand formula. for interventional distributions and hence reuse computation across the set of graphs returned by causal discovery methods. Additionally, the importance weighting method only requires the prediction model to be applicd the initial data, thus reducing the computational load.
We evaluate the performance of our methods using sinulated and real-world data. ingly, we find that even when the true causal graph or its equivalence class are not learned correctly, the explanations returned by our method tend to be closer on average to the true causal Shapley value than the marginal and conditional Shapley values.
2

---

Published in Transactions on Machine Learning Research (03/2026)
2 Background
In local explainability, the goal is to explain the prediction of model f for single input x" (Doshi-Velcz & Kim, 2017) that is an instantiation of the features X = (Xi...., Xp} with distribution P(X). A popular method for explaining the influence of each feature x on the model output f(x*) are Shapley values, (Shapley, 1953), introduced in game theory as a way of distributing value over a set of p 'players' The Shapley value of player for a value function is:
0i(1) = SC
1
E
(P-1) (--1)² |S|
- i [p],
(1)
where [p]= [[,.....)) p} is the set of all players and each S represents subset of [p] \ {i), called coalition. This value can be seen as weighted average of the marginal contribution of i to every possible coalition of other players S. It can be shown that the contribution of each player sums to the value gained when every player contributes, which ensures that the value produced by the whole group together is distributed over the player's individual Shapley values.
To apply Shapley values in local explainability, we consider each feature a "player". It is not obvious how to choose an appropriate value function and many options have been proposed (Sundararajan & Najmi, 2020). In particular, we need to choose what happens features that are in the coalition, denoted by the index set 5, and features that are not in the coalition, denoted by the index set S. A common choice for is the expectation of our model w.r.t. features that are not in the coalition. If we assume independent features, the value function is the expectation w.r.t. the marginal distribution of the features out of the coalition (Strumbelj & Kononenko, 2014; Datta et al., 2016; Lundberg & Lee, 2017). The value function is then:
(S) = Ex(x)))((X3 Xz)..
where xÿg are the features of the input that are in S, while X3 represents the features that are not in S over which we take the expectation. We refer to the Shapley values based on this value function as marginal Shapley values. Aas et al. (2021) take into account dependence between the features with the following value function, which we refer to as conditional Shapley values,
v(S) = -xxxxxxxxxxxxx
(3)
Heskes et al. (2020) propose causal Shapley values, which take into account causal relations between features, represented as a Structural Causal Model (Pearl, 2009; Spirtes et al., 2000), by replacing conditional distributions with interventional distributions. This results in choosing e as follows:
v(S) (Xxx Xz)].
(4)
Here P(X (Xo = x5)), represents the distribution after performing an intervention do(Xs = =x}) that forces Xg to take the value x5. This distribution is in general not the same as the conditional distribution P(X Xxs = x3) (Pearl, 2009). Heskes et al. (2020) show that with this value function, an indirect effect the contribution of feature only occurs when change in feature's value has causal effect the prediction, which leads to an explanation that more accurately reflects the effect of each feature. We provide a short introduction to causality in App. A.1 and an comparison between conditional and causal Shapley values, including a worked out example in our causal discovery setting in App. A.2.1.
3 MEC Shapley Values
In this paper, we exlend the causal Shapley value framework from leveraging known causal relations to the case in which we learn a set of possible causal graphs from observational data, ie. a Markov Equivalence Class (MEC), using causal discovery methods. This results a set of causal explanations, one for cach Markov- equivale causal graph that may have generated the data. Considering set of causal models compounds the computational challenges associated with computing Shapley values. After estimating the
3

---

Published in Transactions on Machine Learning Research (03/2026)
MEC, naive strategy is to iterate over each Directed Acyclic Graph (DAG) in the MEC and compute the corresponding causal Shapley value. However, the size of a MEC grows rapidly with the amount of nodes (He et al., 2015). We therefore want to reuse computation where possible by identifying which DAGs share the same interventional distribution for a subset of features S and by introducing two mothods to estimate the interventional expectation: (i) a sampling method and (ii) an importance weighting method. We describe each step of our pipeline the following.
3.1 Causal discovery and grouping of graphs
fk=1 for features [X] (X,,...,XX, our goal is to explain the prediction of a prediction model f for input data point x'. The first step is learning causal relations between the features from the data D with causal discovery. We assume that the causal Markov and faithfulness assumptions hold, T.e. that the conditional independences in the observational data correspond to d-separations in the true causal graph (Spirtes et al., 2000). We assume that we are in a causally sufficient setting, in which we do not have unobserved confounding or selection bias and the causal graph is a Directed Acyclic Graph (DAG). The output of most causal discovery methods is a Markov Equivalence Class (MEC), i.e. all graphs that induce the same conditional independences and dependences as in the observational data. We denote the MEC G and define it a set of candidate graphs G (Gi,...,Gal For many sets 8, multiple graphs in the MEC share the interventional distribution P(Xyldo(Xs x3)). By identifying which graphs share this distribution we only need to sample from one of them to estimate e(S) [E (X 3. x²) (ooXs = x3)]]- So for each intervention target S we rewrite the interventional distributions with the ID algorithm (Shpitser & Pearl, 2006) as an estimand of observational quantities. We group graphs with the same estimand for the intervention on in a list of set of graphs Gs, where for each group we estimate v(S) once.
Given a dataset D (x
3.2 Sampling method
After learning the MEC and grouping the graphs based on their interventional estimands, as described in Soc. 3.1, we can now use sampling mothod for estimating causal Shapley values, similar the method usod by Heskes al. (2020) for single, known causal graph. For each subset of features S and a graph G3 in each group of Gs, we follow three steps: (i) conditional distribution estimation, (ii) simulation of interventional data, and (iii) estimation of the value function.
First, we consider G; as a potential causal graph and follow the topological order of G; to estimate Pa. (Xp(Xaa,(0)) U)) for each variable Xj, where X () are the variables in X that are the causal parents of Xj in G.. We then generate simulated interventional samples from the interventional distribution Pc. (X|do(Xg = x3)) that assumes that G3 is the true causal graph. In particular, we all of the values of Xs x5S and use the truncated factorization formula (Eq. 11 in App. A.1) for which we sample each factor in topological order in G; as follows:
(K.G. 1mm
II Pc. (X,lXraa ())x(Xs= =x5) jES
= Pa. (Xldo(Xs x5)),
where ó is the Dirac delta function, forcing Xs to take the values x³.. We then use the nime simulated samples in a Monte-Carlo approximation of the value function o for the graph G. as follows:
= G.)
((((..))
NhE k=
Having estimated each interventional expectation, we combine them using Eq. I or, when sampling only a subset of all possible subsets, using Kernel SHAP (Lundberg & Lee, 2017; Aas et al., 2021).
3.3 Importance Weighting (IW) method
The sampling method from See. 3.2 may incur high computational cost, since it simulates new datapoints for each graph and intervention, each of which must be evaluated by potentially expensive model. To
4

---

Published in Transactions on Machine Learning Research (03/2026)
improve the computational efficiency, we introduce importance-weighhtngg method that re-uses single set of predicted samples for an intervention on X5 across all DAGs in the MEC. The pseudocode for the method is provided in Alg, 1. While this method avoids having to evaluate new simulated data points, this approach can be susceptible to high variance if the target interventional distribution differs significantly from the observational data, a known limitation of importance sampling that the sampling- based method avoids.
We start by defining the set of subsets of feature that we consider S similar KernelSHAP (Lundberg & Lee, 2017) (Lines 1-2). We learn the MEC the initial data D using standard causal discovery methods and group the graphs based
their interventional estimands, described in Sec. 3.1 (Lines 3-9). For each subset of features S E S we then compute the value function v(S) for each group of graphs (Lines 10-19), which we describe in detail in the following, and apply KernelSHAP (Lundberg & Lee, 2017) following Aas et al. (2021) to compute the final Shapley values (Lines 20-24).
The core of this algorithm is the computation of the value functions for cach group of graphs in an efficient way. The first step is to compute the predictions f(x (A) xz)) using the prediction model f and combination of the initial data D and the input datapoint x* (Line 11). We reweight these predictions based on the interventional distribution in each graph G; and estimate the expected value Ep,, (x,,dxoldo(xs=xx (Xx]do(x==x Gg), which is the value function used in causal Shapley values (Heskes et al., 2020). We rewrite vSS) for graph G, as weighted expectation over samples from the marginal distribution P(Xg), where 3 are the features not in S:
v (i)()) =Bpe,
(x33 X 3) |G.
(5)
= / [(xs ) (X xXs = =x3)) P(Xg)
5)
dXg
(6)
P(X 5)
Pa. (x sllo((xs = -x3))
= B
x(x3,
Xz)
(7)
P(Xg)
where in Eq. 6 we introduce P(X5) in the denominator and the numerator. To approximate this expectation, we first estimate the conditional densities Pa. (Xx(XPae;()) (i) from data D (Line 14). We then re-weight each
prediction f(x
(k) x3g) for each factorization in a Monte Carlo approximation of Eq. 7 as follows:
~ (S))
1
D
1
Σ
D
f(x
(k)
xx3)
Pa. (X (xgldo(xs (ooss =x))
(8)
D
P(Xs)
Po. (Xgldo x3)) = III Pa. pac, Pac, W x Xpao, (nn))
JEB
following the truncated factorization formula (Eq. 11), where we force all variables in S to take the value x5. including the parents of cach variable X, in G,, represented as Pac: (j) n.. We compute Eq. 8 for a single graph in each group and assign it to the other graphs in the same group (Lines 15-16). For each graph G; we use v(S); for each set of features S to compute the Kernel SHAP values (Lundberg & Lee, 2017) (Lines 19-22). We return a list of Kernel SHAP values, one for cach DAG in the MEC G (Line 23).
4 Related work
Connections between causality and explainability have been explored in various works (Beckers, 2022; Karimi et al., 2023). Many of these focus leveraging causal knowledge to improve explainability (Khademi & Honavar, 2020; Schwab & Karlen, 2019; Oesterle et al., 2023), but they do not use the framework of Shapley values, which provides desirable properties and computational tools in terms of interpretability and efficiency. Some approaches employ causal discovery for explainability outside the Shapley framework. For example, Sani et al. (2020) use causal discovery for explainability of black-bor methods. In contrast to our work, they consider unstructured data such as images and allow for potential latent confounding, but do not quantify causal effects. Similarly, Takahashi et al. (2024) also perform causal discovery for local explanations, but focus the context of counterfactual explanations.
5

---

Published in Transactions on Machine Learning Research (03/2026)
Algorithm 1 Importance weighting method for Estinating MBC Shapley Values
Input: Dataset D = {2=, prediction model f, target variable Y, data point to explain x", number of combinations consider Recomb I |2|0|| for p features.
Initialize set of subsets of features to consider S < {0, [p]} following (Lundberg & Lee, 2017).
Sample Neomb 2 subsets of the power set 2 \ {0, [p] with probability equal to kernel weights k(S)
(w-1) 1)/ (($) [S))s 1s1) and add them to S.
3: Learn the MEC G - [GI,...,Ge] by applying a causal discovery method to D
4: for each subset. 5 E S of features do
5: for each graph G2 E do
Identify estimand formula for pG; (Xxldo(Xx x3)) with ID algorithm. end for
8: Group graphs with same estimand formula for the intervention on S in a set of graphs G's.
end for
10: for each subset S S S of features do
1D 11: Compute predictions {f(x (k) xg)]k== with observational data (x5 1k=1 and input x3-
12:
for each group of graphs G's with the same estimand formula for P(Xgldo(Xs = xÿ)) do
13:
Choose any Gj in G'g-
14:
Estimate conditional probabilities Pa. (X, X Pac, (n) from observational data D. Calculate value of v for graph G, and set S with Eq. 8:
15:
Ihes P(x; (K) 1x (2)
Paa, (1) X (w)ns)
(S) D
1(x x's)
P(x
16: Store the value of (S) for the value functions of all other graphs in G's-
17: end for
18: end for
19: for each graph G2 E G do
vc. ((r))) 21: Compute the Kernel SHAP value SG for with Ga following (Lundberg & Lee, 2017).
20: Collect all value functions for each subset S C S a vector vG, <-
22: end for
23: return List of Kernel SHAP values f for each DAG in the MEC G.
Several works have combined Shapley values with causal reasoning. This work is closely related to that of Jung et al. (2022), who adopt the same value function in the context of causal contribution analysis, and propose a weighted estimator to estimate the value functions in the case of discrete data, Janzing al. (2020) argue that, since are explaining that interventional Shapley values become equivalent to marginal Shapley values. As shown by Heskes et al. (2020), this leads to explanations that do not consider indirect effects of the features on the predictions, e.g., when we use proxy of a feature instead of the feature itself, as discussed also by Frye et al. (2020) in the context of unresolved discrimination.
algorithm, the features should be seen independent, such
Frye et al. (2020) introduce asymmetric Shapley values, incorporating causal knowledge into explanations by weighting subsets according to their agreement with a partial causal ordering. They do not make use of the concepts of interventions, but noted by Heskes et al. (2020) they can be combined with causal Shapley values, specifically when making use of chain graphs. Budhathoki et al. (2022) and Strobl & Lasko (2022) apply causality to Shapley values, but in the context of root cause analysis, specifically to explain outliers instead of general predictions by any arbitrary model given a known causal graph. These works only leverage known causal relations, either as a partial causal ordering or a complete causal graph. An exception is the work by Strobl & Lasko (2022) that learn the causal graph, but assume a linear non-Gaussian model, a more specific setting than ours which leads to a single causal graph instead of MEC. Concurrently our work, Ng et al. (2025) propose Causal SHAP, which integrates the PC and IDA algorithms into SHAP to yield
6

---

Published in Transactions on Machine Learning Research (03/2026)
single attribution vector corrected by estimated causal strengths. In contrast, our framework is agnostic to the specific causal discovery method and explicitly considers the entire Markov equivalence class, producing a range of causal Shapley values that reflect multiple plausible causal structures.
5 Experiments
We evaluate our methods across synthetic datasets and two real-world datasets. One of the challenges when evaluating methods for explainability is that we lack agrood-upon metrics for asscssing the quality of explanations. We therefore focus on a tangible, well-defined objective: evaluating how well and how efficiently our methods approximate the theoretical causal Shapley values derived from ground-truth causal graph. These values provide principled framework for attributing model's prediction to its input features based on their causal effects.
In the synthetic datasets we have access to the true interventional distributions, allowing us to compare the output of our methods to the true causal Shapley values. We evaluate in a linear Gaussian and a non-linear setting. In addition to comparing with the ground truth causal Shapley value defined by the true interventional distributions, we also compare the outcome of our method with the marginal and conditional Shapley values, which do not take any causal relations into account. To implement our method we adapt the shapr package by Aas el. al. (2021). For causal discovery, we used the pcalg package (Kalisch el al., 2012) and the Tetrad toolbox (Ramsey et al., 2018). For identifying the interventional estimand formulas we use the causaleffect package (Tikka & Karvanen, 201 7a;b). For prediction we use XGBoost models (Chen & Guestrin, 2016). We provide the code to reproduce our experiments in the supplementary material.
5.1 Synthetic data
We generate data for two settings, linear Gaussian and non-linear causal models. For each setting, we consider Nsem 40 randomly generated causal models for which we generate synthetic data as follows.
For each model we first sample an Prdos-Renyt graph of size p 5, 10, 15 where we specify the expected number of neighbors d 2, for each node. We add our 'dependent" node Y, the target variable for our prediction model, for which we specify the minimum number of incoming edges (mmin 2, 5) and the expected number of incoming edges (m = 3, 6.
In the linear Gaussian setting, nodes are parametrized as X, = Ejera, WiyX, + e where Wy ~ U((-2, 2, -05) U (0. and & = (N, 1). For the non-linear setting we simulate data with the equa- tions X; who (EEps, w xx) +G where denotes the sigmoid function, w/y ~ U(-11.5, -0.5) (0.5, 1.5),, w² ~ ((-3, (1, 3)), and e; follows either a (0, 1) or a U--1, 1) distribution with equal probability. We apply causal discovery to standardized data. For each cansal model we simulate nest = 40 data points x" for which we compule the Shapley values.
Methods. For cach model, following Heskes et al. (2020) we train XGBoost model (Chen & Guestrin, 2016) subset of Mraia = 10K samples for 100 rounds neet = 40 data points x" which are unseen by the model. We compare the results of our methods with Marginal (Tundberg & Tee, 2017) and Conditional Shapley values (Aas et al., 2021). For the Conditional Shapley baseline, we estimate the value function v(S) = =[f(x) Xs = 231 by modeling the feature distribution P(X) as a multivariate Gaussian, where the conditional distributions P(Xs] | Xs) are derived analytically from the observational sample mean and covariance. This approach, implemented via the shapr package (Aas et al., 2021), serves as the standard baseline for dependent features but assumes linear relations between variables. We use the same conditional estimation for conditional and MEC Shapley values for a fair comparison. To evaluate whether our method provides benefits beyond improving conditional estimation, we additionally compare against conditional Shapley values estimated with Gaussian Copula and Conditional Inference Trees in App. A.3.3.
=zg]
As causal discovery methods, we use the PC algorithm (Spirtes et al., 2000) and the FGES algorithm (Ramsey et al., 2017). We call the combination of our method with PC the PC MEC Shapley method, while the combination with FGES is the FGES MEC Shapley method. We compare the results of these algorithms with causal discovery oracle, which we call the Oracle MEC method. Crucially, this is partial oracle: it is
7

---

Published in Transactions on Machine Learning Research (03/2026)
Method
M Marginal
Conditonal
Oaacle MEC
PC MEC
FGES MEC
Estimation
Sampling
IW
0.6
0.12
NL
0.4
0.09
0.2
0.06
10 Number of nodes p (a) Linear Gaussian
15
Number of nodes (b) Non-linear
Figure 2: Normalized L2 error between the ground truth causal Shapley value and the Shapley values returned by each method for different parameters: (a) linear Gaussian setting with varying number of features p, (b) nonlinear Gaussian sctting with varying number of features p, (c) lincar Gaussian sctting with varying number of expected neighbors d, and (d) linear Gaussian setting with varying expected number of nodes influencing Y m. For all MEC Shapley methods we plot the NL2E over the range of explanations, iee. the minimum and maximum of NL2E across all DAGs G, in the MEC G, averaged over all nc = 40 causal models. The marker represents the average NL2E across the DAGs in the MEC and models.
given the true underlying graph structure (and thus the true MEC), but must still estimate the required conditional distributions from the same nobs samples as the other methods. For all methods, we consider Sampling version, described in Sec. 3.2, and an IW version, described in Sec. 3.3.
We apply the PC algorithm to the linear data using partial correlation tests and for non-linear data we apply kernel-based conditional independence Lests (Zhang et al., 2011), both with significance threshold a = 0.05. We apply FGES in both settings using the BIC criterion. For the conditional and MEC Shapley methods, we model conditional distributions using Gaussian parametric assumption. In the non-linear setting, we additionally evaluate conditional Shapley values using two more flexible conditional estimators (Gaussian copula and conditional inference trees) in App. A.3.3. We approximate marginal Shapley values using samples from empirical marginal distribution p" (Xg). By default we Red 1000 observational samples for PC and FGES, Robs = 1000 samples to estimate conditional expectations and Nmme 1000 samples to estimate the Monte Carlo approximation. For 5, 10 we consider the number of subsets 7cmmb = 2", meaning that we use the full set of of combinations S. As the set of combinations increases exponentially with the amount of nodes, for p = 15 we set Acoomb = 8192 for the ground truth and Neomb = 4096 for the other methods. We provide ablations for these parameters in App. A.3.
Metrics. We evaluate how close are explanations from different methods in terms of agreement with the ground truth causal Shapley values &e which we compute using the ground truth graph G and sampling from the true interventional distributions with n 4000 Monte Carlo samples for each datapoint x".
We then compare on the L2 error between the Shapley value returned by each method and the ground truth causal Shapley value, normalized based on the its value for the datapoint to allow comparison across different. data points, defining the normalized L2 crror as follows:
-meethod l2
2.
(9)
10 ll2
where dhe refers to the ground truth causal Shapley value for x" and @method is the Shapley value returned by each method. Since MEC Shapley values methods provide a list of possible explanations, one for each causal graph in the MEC G [G).... G.), this results set of explanations (NI2EG,..... NL2Eg. }. Note that this is also obviously true for the Oracle MEC Shapley values method, since the ground truth MEC includes multiple DAGs. To visualize this compactly our results, we consider the minimum and maximum NL.2F across
8

---

Published in Transactions on Machine Learning Research (03/2026)
Method
Marginal
Conditional
Oracle MEC
PC MEC
FGES MEC
Estimation
Sampling
IW
0.8
0.75
NL2E
NL2E
0.6
0.50
0.4
0.25
0.2
Number of expected neighbors d
Number of expected parents of Y m
(a) Linear Gaussian
(b) Linear Gaussian
Figure 3: Normalized 1.2 error between the ground truth causal Shapley value and the Shapley values returned by each method for different parameters: (a) linear Gaussian setting with varying number of expected neighbors d, and (b) linear Gaussian setting with varying expected number of nodes influencing Y m. For all MEC Shapley methods we plot the NL2E over the range of explanations, i.e., the minimum and maximum of NL.E across all DAGs G, in the MEC G, averaged over all Nsa = 40 causal models. The marker represents the average NL2E across the DAGs in the MEC and models.
all DAGs in a MEC as the NL2E over the range of the explanations, i.e., [ming,ec NL.2FG,, maxG,eG NL22a,,)- Moreover, we report the average NL2E over the DAGs in the MEC, i.e., EG
Results. We report the results for all the synthetic data settings in the following. Fig. 2a shows the NL2E between each method and the ground truth causal Shapley as function of the number of nodes p 5, 10, 15 for linear Gaussian causal models for d = 2 and m = 3, 6,9.
For MEC Shapley values methods, the bars represent the NL2E over the range of explanations, so the minimum and maximum NL2E for each DAG a MEC, while the marker represents the average NT,2F over all DAGs in a MEC, averaged over all causal models. Since the baselines cntargiaal/condiiionn instead produce a single Shapley vector, we then report the average NL2E over all causal models. For completeness, we also report the standard error over the causal models for the same setting in A.3.2, showing that our results are stable across the different models.
Shapley)
We observe that as expected the estimated causal Shapley values resulting from the oracle MEC are closest to the ground truth causal Shapley values, but they still often do not include the ground truth causal Shapley values, because while the Oracle MEC has access to the ground truth MEC, it still uses the same nots finite sample data as other methods to estimate the conditional distributions. We can see that the average NL2E for all the MEC Shapley methods is still notably lower than the baselines, although the improvements are limited by the accuracy of causal discovery finite sample cases. For the oracle MEC and the FGES MEC can that even the worst explanations in the range can have comparable or better performances than the baselines. As expected, the Oracle MEC method consistently has the lowest. NT.2F, but as expected still has a relatively large spread of the NL2E due to the variability across the DAGs in the ground truth MEC. When comparing the performance of FGES and PC, we observe that the resulting explanation quality (NL2E) does mostly, but not always correlate directly with the global graph structure similarity (SHD). For instance, in the linear setting with 16 nodes (6=16), FGES yields explanations with lower error than PC (Fig. 2a), despite having a slightly higher average SHD (Appendix Table 3).
For non-linear causal models we see in Fig. 2b that all our methods provide better explanations average, and that in some cases, even the worst explanation in our range is better than the basclines. Finally, we see a slight difference between Sampling and IW, with Sampling outperforming IW. Surprisingly, here the marginal Shapley values perform better relative to the true causal Shapley values than the conditional ones.
Fig. 3a shows the NL2E between cach method and the ground truth causal Shapley as a function of the expected number of neighbors d 1, 2,3 3 for linear Gaussian models with p = 10 and m 6. As expected the oracle MFC Shapley is closest to the ground truth, and is not significantly affected by the expected
9

---

Published in Transactions on Machine Learning Research (03/2026)
Amnrtbbed Time (ma) Time (mis. leg ocall)
82
16
10
20
1.5
0.9
Number DAGs
MEC
FGES Sample Total (wl ID)
FOES Sample Compule-onty (wlo ID)
Conditinal
Marginal
Figure 4: Breakdown of computation time for 40 SCMs (scatter points) in the lincar-Gaussian regime. (top) The absolute time for our method (solid lines) compared to baselines (dashed/dotted) (bottom) The amortized time per explanation, showing that our method becomes more efficient than the conditional baseline as the number of DAGs in the MEC increases. The time measured is the computation time (total amortized) for computing the Shapley values for 40 datapoints. We plot LOESS-smoothed curves to highlight the trend. The plots for the Oracle and PC MECs look similar.
number of neighbors, which is as we can see for d = 3, not the case for PC and FGES MEC, but they on average still show comparable or closer agreement than the baselines. As the expected number of neighbors d increases, the baselines perform comparatively worse. The difference between Sampling and TW is again minimal. Fig. 3b shows that in the linear setting m, the expected number of parents, has no strong effect on performance.
In App. A.3 Fig. 8 we show how the methods compare when we consider only the cases where PC and FGES output an incorrect CPDAG. We can see that overall, performance not strongly affected. In almost all settings the average explanation closer than the basclines We show the mean and standard errors for these results in Fig. 9. We consider d and m for the non-linear case in Fig. 11, where we see that in most cases the worst explanation by both PC and FGES still outperforms the basclines. We provide ablations in terms of number of Monte Carlo samples tme in Fig. 12, which do not seem to influence the results substantially. We also vary the number of samples to estimate the conditional distributions node in Fig. 13, the number of combinations used to compute the Shapley values Reoms in Fig. 14 and the number of samples used to estimate the MEC ned in Fig. 15. As expected all mothods perform slightly better with more samples and higher numbers of combinations, but the differences bet ween the baselines and the MEC Shapley methods are qualitatively similar as what presented in the previous plots. In App.A.3.10 we show an analysis of the learned MECs in all scttings, reporting average sizes for PC and FGES and the quality of the CPDAGs compare to the oracle CPDAGs as measured by the Structural Hamming Distance (SHD). In Fig. 17 show how the SHD affects the wit range of NT.E. Additionally, in App. A.3.3 in the non-linear setting, we compare with conditional Shapley values with two more flexible conditional estimators (Gaussian copula and conditional inference trees) and show that the results of the baselines improve only marginally, while our methods, even with the conditional Gaussian estimators, stilll outperform the baselines.
Figure 4 evaluates the average time to compute the Shapley values for 40 data points, comparing the FGES Sample method against the marginal and conditional baselines as function of the number of DAGs in the learned MEC. We report per-SCM computation times (scatter points), and overlay LOESS-smoothed curves to highlight the overall trend. The top panel illustrates the absolute computation time. As expected, our full method is more computationally expensive than the baselines, as it must compute explanations for multiple DAGs. Furthermore, it incurs a significant, but fixed, one-time overhead for the ID algorithm identify opportunities for amortization (Total w/ ID). Crucially, however, we observe empirically that the total runtime scales favorably, growing sub-linearly with the MEC size, relative to a naive per-DAG iteration baseline, which would scale linearly with the amount of DAGs in the MEC. This suggests the
10

---

Published in Transactions on Machine Learning Research (03/2026)
ADNI Dataset
0.04
0.02
Marginal
Conditional
FGES MEC (W)
FGES MEC (sampling)
PC MEC (WW
PC MEC (samplinng)
Method
Figure 5: Average NL2E w.r.t the ground truth generated using the "gold standard" graph from Shen et al. (2020) on the ADNI dataset. Standard errors are too small to be visible.
Scaled Resture valuc
nigh
Scaled "cature
High
ccaled footure
408
FDG
ABETK
ABETA
774
AF0541
4P064]
AP0541
APOE42
AP0E42
AP0E42
-4:0
-0.23
0.25
0.2
-0.0
-0.2
0.02
225
-030
-0.25
350
Conditioal Shaplen value
Gold stardard cousel Shaple wile
FGES MEC causal Shaplel wabe
(a)
(b)
(c)
Figure 6: Sina plots for binary classification task using XGBoost the Alzheimer's dataset, showing the (a) conditional Shapley values, (b) causal Shapley values for the 'gold standard' DAG, (c) FGES MEC Shapley values. Sina plots for these and also the other methods are in Fig. 19.
amortization identified by the ID step can materially reduce average cost per explanation as the MEC size grows, though this is an empirical observation rather than complexity bound. The plot also decomposes the cost, revealing that while the ID algorithm is a major component, the core explanation cost (Compute-only w/o ID) scales even more elficiently. The bottom panel reveals the practical benefits of this structure by showing the amortized time per DAG. As the MEC size increases, this amortized cost continues to fall and begins to flatten, demonstrating that the marginal cost of explaining an additional causally equivalent graph becomes negligible. Even when including the initial ID overhead, the true average cost per explanation (Amortized Total w/ ID) consistently decreases, highlighting the power of amortizing computation.
In summary, our framework presents computational trade-off: higher, fixed initial cost to identify shared causal structures is exchanged for per-explanation efficiency and scalability as the causal uncertainty (i.c., MFC size) increases.
5.2 Real-world data
We apply our methods to two real-world settings: the Alzheimer's discase datasct obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (http://adit paaadi.lome .loni. .uc. edu) and to the bike rental dataset from (Fanaee-T & Gama, 2014), both used in Heskes et al. (2020). We provide an explanation of the datasets and detailed results in App. A.4. For the Alzheimer's datasct we consider causal Shapley values based on a "gold standard graph" from expert knowledge (Shen et al., 2020) as ground truth, where trained prediction model to predict the probability of Alzheimer's diagnosis.
11

---

Published in Transactions on Machine Learning Research (03/2026)
date
2012-10-09
2012-12-03
Marginal
Conditional
Causal (chain)
PC
Shapley value
cosyear
temp
cosyear
temp
cosyear
temp
cosyear
temp
Figure 7: Bar plot for the features cosyear and temp of the bike rental dataset for two days in October and December with approximately the same lemperature. We show the Causal chain Shapley values (Heskes et al., 2020), along with marginal, conditional, and PC Shapley values. FGES Shapley values are similar PC and have been omitted.
Fig. 5 shows how the outcomes the different algorithms compare to the ground truth Shapley values, as well as the baselines. The output both FGES and PC is a MEC with a single graph, leading to a single causal Shapley value for each feature. We see that the FGES MEC methods are closer the ground truth causal Shapley values than the baselines, especially when using the sampling method. The PC MEC methods are comparable to marginal and conditional Shapley values. This result highlights the critical interplay between the causal discovery and explanation stages of our pipeline. As shown in Appendix Fig. 18, the CPDAG learned by PC only learns correctly two edges, missing two edges and incorrectly orienting the edge between fideoxyglucosc (FDG) and amyloid beta (ABETA). On the other hand, the CPDAG learned by FGES correctly identifies four edges, missing only the connection wetween age and amyloid beta (ABETA). These errors in the graph structure likely explain why the PC MEC Shapley values are less aligned with the ground truth. This suggests that the benefit of our framework is directly tied to the ability to learn reasonably accurate causal structure. When the discovery step is successful (as with FGES), the resulting explanations show notable improvement. We also shwo the sina plots for ional, ground truth causal Shapley values and FGES MEC Shapley valuesin Fig. 6. The sina plots for all methods are in Fig. 19.
In Fig. 7 we look at the Shapley values for predictions of bike rental counts on two different dates with approximately equal temperature. Here the aim qualitative illustration, line with Heskes et al. (2020), showing how different Shapley formulations can yield distinct explanations for comparable prediction scenarios. While marginal Shapley values provide more less the same explanation for both days, focusing mostly on the temperature variable, temperature, resulting in different explanations for the two days. In App. A.4 we show additional sina plots for both datasets that show how the features get assigned different levels of importance by each method, in addition to providing more details on the real-world experiments.
see that causal and PC Shapley values assign credit to both season and
6 Conclusion and limitations
In this paper, we combined causal discovery with causal Shapley values to provide a more nuanced explanation of predictive models, particularly when the causal graph among features is not known. Our experimental results show that our methods provide explanations that are more closcly aligned with the true causal impacts of features on predictions, compared to traditional Shapley value approaches that do not incorporate causal reasoning at limited computational cost. Our results suggest that even when the causal graph is not in the set of learned graphs, the explanations our methods output are often closer to the ground truth causal Shapley values than other methods.
While our approach marks a step forward, several limitations need to be acknowledged. First, our claims are limited the accuracy and computational efficiency with which we can approximate the true underlying causal Shapley values. Whether causal Shapley values offer a better explanation than those provided by other methods requires interpretation, and will depend the particular analysis use under consideration.
12

---

Published in Transactions on Machine Learning Research (03/2026)
Since our method returns set of possible explanations (one per DAG in the MEC), practitioners should interpret these as reflecting causal uncertainty rather than a single definitive attribution.
A computational limitation is the scalability of our pipeline in terms of number of nodes and graph size. In particular, we are limited by the computational cost of causal discovery, which is generally an NP hard problem (Chickering et al., 2004), and even more so by Shapley estimation (which is exponential in the number of features without coalition subsampling). In terms of accuracy, we are dependent on the causal discovery methods that is employed, which can be inaccurate under wrong parametric assumptions for small data samples. Similarly, the accuracy of both the Sampling and IW method depend the correctness of the estimated conditional distributions. Our methods also assume causal sufficiency, i.c. no unmcasured confounders or selection bias, which may not always hold in real-world settings.
Acknowledgments
We thank SURF for the support in using the Dutch National Supercomputer Snellius. DR and JWvdM acknowledge the support of the Bosch Center for Artificial Intelligence. SM would like to acknowledge the VIDI grant CANES (VI.
References
Kjersti Aas, Martin Jullum, and Anders Loland. Explaining individual predictions when features are dependent: More accurate approximations to Shapley values. Artificial Intelligence, 298:103502, 2021.
ISSN 0004-3702. doi: 10.1016/jaartint.222.. 202..10350.
Sander Beckers. Causal explanations and XAI. In Bernhard Schölkopf, Caroline Uhler, and Kun Zhang (eds.), Proceedings of the First. Conference on Causal Learning and Reasoning, volume 177 of Proceedings of Machine Learning Research, pp. 90-109. PMLR, April 2022.
Kailash Budhathoki, Lenon Minorics, Patrick Bloebaum, and Dominik Janzing. Causal structure-based root cause analysis of outlicrs. In Kamalika Chaudhuri, Stefanic Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th International Conference on Machine Learning, volume PMLR, July 2022.
162 of Proceedings of Machine Learning Research, pp.
Diogo V. Carvalho, Eduardo M. Pereira, and Jaime S. Cardoso. Machine learning interpretability: A survey methods and metrics. Electronicsweek; 8(832), 2019. ISSN 2079-9292. doi: 10.3390/electronici elerroooscrrooss0000332222222
Tianqi Chen and Carlos Guestrin. XGBoost: A scalable tree boosting system. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Kdd '16, pp. 785 794, New York, NY, USA, 2016. ACM. ISBN 978-1-4503-4233-2. doi: 10.1145/293397722233377
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Mochine Learning Research, 3(mull):507-554 March 2003. ISSN 1532-4435. doi: 10.1662/11524449322188777
David Maxwell Chickering, David Heckerman, and Christopher Meek. Large-sample learning of bayesian networks np-hard. J. Mach. Learn. Res., 5:1287-1330, Docember 2004. ISSN 1532-4435.
Diego Colombo and Marloes H. Maathuis. Order-indepeudent constraint-based causal structure learning. Journal of Machine Learning Research, (921116))322-3922 2014.
Gabor Csardi and Tamas Nepusz. The igraph sore package for complex network research. InterJournal, Complex Systems:1695, 2006.
Anupam Datta, Shayak Sen, and Yair Zick. Algorithmic Transparency via Quantitative Input Influence: Theory and Experiments with Learning Systems. In 2016 IEEE Symposium on Security and Privacy (SP), pp. 598 617, San Jose, CA, May 2016. IEEE. ISBN 978-1-5000-0824-7. doi: 10.1109/SPP2016..22
Finale Doshi-Velez and Been Kim. Towards A Rigorous Science of Interpretable Machine Learning, March 2017.
13

---

Published in Transactions on Machine Learning Research (03/2026)
Hadi Fanace-T and João Gama. Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence, 2:113-127, June 2014. doi: 00.1007/s13748--13-000/44
Matteo Fasiolo. An Introduction to Menfast. R Package Version 0.2.8., 2014.
Christopher Frye, Colin Rowat, and llya Feige. Asymmetric Shapley values: Incorporating causal knowl- edge into model-agnostic explainability. In Proceedings of the 94th International Conference on Neural Information Processing Systems, NTPS'20, Red Hook, NY, USA, 2020. Curran Associates Inc. ISBN 978-1-71382-955--6
Christopher Frye, Damien de Mijolla, Tom Begley, Laurence Cowton, Megan Stanley, and Lya Feige. Shapley explainability on the data manifold. International Conference on Learning Representations, 2021.
Steven B. Gillispie and Michael D. Perlman. Enumerating markov equivalence classes of acyclic digraph models. UAI, 2001.
Yangbo He, Jinzhu Jia, and Bin Yu. Counting and exploring sizes of Markov equivalence classes of directed acyelic graphs. Journal of Machine Learning Research, 16(1):2589 2609,, January 2015. ISSN 1532-4435.
Tom Heskes, Toan Gabriel Bucur, Evi Sijben, and Tom Claassen. Causal shapley values: Exploiting causal knowledge to explain individual predictions of complex models. In Proceedings of the 34th International Conference on Newral Information Processing Systems, NIPS'20, Red Hook, NY, USA, 2020. Curran
Associates Inc. ISBN
Dominik Janzing, Lenon Minorics, and Patrick Bloebaun. Feature relevance quantification in explainable AI: A causal problem. In Silvia Chiappa and Roberto Calandra (eds.), Proceedings of the Twenty Third International Conference on Artificial Intelligence and Slatislics, volume 108 of Proceedings of Machine Learning Research, pp. 2907- 2916. PMLR, 2020.
Yonghan Jung, Shiva Kasiviswanathan, Jin Tian, Dominik Janzing, Patrick Blocbaum, and Elias Barcinboim. On measuring causal contributions via do-interventions. Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 10476 10501. PMLR, July 2022.
Markus Kalisch, Martin Mächler, Diego Colombo, Marloes H Maathuis, and Peter Bühlmann. Causal inference using graphical models with the R package pcalg, Journal of Statistical Software, 47(11):1-22) 26, 2012.
Moir-Hossein Karimi, Krikamol Muandet, Simon Kornblith, Bernhard Schölkopf, and Been Kim. On the relationship between explanation and prediction: A causal view. In Proceedings of the 40th International Conference on Machine Learning, ICML'23, Honolulu, Hawaii, USA,, 2023. JMLR.org.
Aria Khademi and Vasant Honavar. A Causal Lens for Pecking into Black Box Prodictive Models: Predictive Model Interpretation via Causal Attribution, August 2020.
Scott M. Lundberg and Su-In Lee. A unified approach to interpreting model predictions. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS'17, pp. 4768-4777, Rod Hook, NY, USA, 2017. Curran Associates Inc. ISBN 978-1-5108-6096-4
Woon Yee Ng, Li Rong Wang, Siyuan Liu, and Xiuyi Fan. Causal shap: Feature attribution with dependency awareness through causal discovery, 2025. URL http://axxiiv.oorraa/////// 00846.
Michael Oesterle, Patrick Blöbaum, Atalanti A. Mastakouri, and Elke Kirschbaum. Beyond Single-Featuur Importance with ICECREAM, July 2023.
Judea Pearl. Causality: Models, Reasoning, and Inference. 2009.
14

---

Published in Transactions on Machine Learning Research (03/2026)
Joseph Ramsey, Madelyn Glymour, Ruben Romerr and Clark Glymour. A million variables and more: The Fast Greedy Equivalence Search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images. International Journal of Data Science and doi: 1111007/41060-0166
Analytics, 3(2):121-129, March 2017. ISSN 2364-415X, 2364-4168.
Joseph D Ramsey, Kun Zhang, Madelyn Glymour, Ruben Sanchez Romero, Biwei Huang, Savini Samarasinghe, Elizabeth A Barnes, and Clark Glymour. TETRAD A TOOLBOX FOR CAUSAL DISCOVERY. 2018.
Numair Sani, Daniel Malinsky, and Llya Shpitser. Explaining the behavior of black-bar prediction algorithms with causal learning. CoRR, abs/ abss/2006.022 02482, 2020. URL https: / //axxivv////i//////// org/ wwgg/abs/2006..o 02482.
Patrick Schwab and Walter Karlen, CXPlain: Causal explanations for model interpretation under uncertainty. In Н. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché-Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Lloyd IS Shapley. A Value for n-Person Games. In Harold W. Kuhn and Albert W. Tucker (eds.), Contributions to the Theory of Games II, pp. 307 317. Princeton University Press, Princeton, 1953.
Xinpeng Shen, Sisi Ma, Prashanthi Vemuri, Gyorgy Simon, and Alzheimer's Disease Neuroimaging Initiative. Challenges and opportunities with cansal discovery algorithms: Application to Alzheimer's 's pathophysiology. Scientific Reporls, 10(1):2975, February 2020.
lyy Shpitser and Judea Pearl. Identification of conditional interventional distributions. In Proceedings of the Twenty-Second Conference on Uncertainty in Artificial Intelligence, UAI06, pp. 437-444, Arlington, Virginia, USA, 2006. AUAT Press. ISBN 0-9749039-2-2.
Peter Spirtes, Clark N. Glymour, and Richard Scheines. Causation, Prediction, and Search. Adaptive Computation and Machine Learning. MIT Press, Cambridge, Mass, 2nd ed edition, 2000. ISBN 978-0-262- 19440-2.
Eric V. Strobl and Thomas A. Lasko. Identifying patient-specific root causes of disease. In Proceedings of the 13th ACM International Conference on Bioinformatics, Computational Biology and Health Informaties, pp- 35445553.
1 10, Northbrook Illinois, August 2022. ACM. ISBN 978-1-4503-9386-7 doi:
Erik Strumbelj and Igor Kononenko. Explaining prediction models and individual predictions with feature contributions. Knowledge and Information Systems, 41()):447-665, December 2014. ISSN 0219-1377, 0219-3116. doi: 10.1007/110115-013-00777
Mukund Sundararajan and Amir Najmi. The many Shapley values for model explanation. In Proceedings of the 37th International Conference on Machine Learning, ICML JMLR.org, 2020.
Daisuke Takahashi, Shohei Shimizu, and Takuma Tanaka. Counterfacttal explanations of black-box machine learning models using causal discovery with applications to credit rating, 2024. URL attps://arxiv.orgg
0tt/22402.02677
Santtu Tikka and Juha Karvanen. Identifying causal effects with the R package causaleffect. Journal of Statistical Software, 76(12):1 30, 2017a. doi: 10.18637/jss.v07..12
Santtu Tikka and Juha Karvanen. Simplifying probabilistic expressions in causal inference. Journal of Machine Learning Research, 18(36):1-30, 2017b.
Hadley Wickham, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D' D'Agostino McGowan, Romain François, Garrett Grolemund, Alex Hayes, Lionel Henry, Jim Hester, Max Kuhn, Thomas Lin Pedersen, Evan Miller, Stephan Milton Bache, Kirill Müller, Jeroen Ooms, David Robinson, Dana Paige Seidel, Vitalie Spinu, Kohske Takahashi, Davis Vaughan, Claus Wilke, Kara Woo, and Hiroaki Yutani. Welcome to the tidyverse. Journal of Open Source Software, 4(43):1686, 2019. doi: 10.21105/josss.16688
Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kernel- based conditional independence test and application in causal discovery. In Proceedings of the Twenty-Seuenth Conference on Uncertainty in Artifiicial Intelligence, UAP11,, pp. 804-813, Arlington, Virginia, USA, 2011. AUAI Press. ISBN 978-0-9744039---22
15

---

Published in Transactions on Machine Learning Research (03/2026)
Yu Zhang, Peter Tiño, Ales Leonardis, and Ke Tang. A survey neural network interpretability. IEEE Transactions on Emerging Topics in Compulational Intelligence, 5:726-742. 772, 2020.
Jianlong Zhou, Amir H. Gandomi, Fang Chen, and Andreas Holzinger. Evaluating the quality of machine learning explanations: A survey on methods and metrics. Electronicsuree; 10(593), 2021. ISSN 2079-9292. doi: 10.3390 10.3390/electronnin /eee//ooooroic1100055333
A Appendix
A.1 Causality background
We introduce a short. background in causality. We assume that the dala-generating process can be modelled using a Structural Causal Model (SCM) (Pearl, 2009; Spirtes et al., 2000), such that each variable X₃ for E [p] has assignment,
X; (XP((()) Ut), i i.
(10)
where Xpa(t) are the parents, i.c. the variables that cause X; directly, and Uj is noise variable, with each U, being independent of the other noise variables. This SCM induces a joint probability distribution P(X) = I P(X,(Xra(t), as well as graph G. If the graph contains no cycles, we call it a directed acyelic graph (DAG). A do-intervention, denoted by d = xs), is a replacement of the causal mechanism for the variables Xs with constant Xs. This results in a new distribution P(X]do(Xg = =xs) that factorizes according to the truncated factorization formula:
P[X]do(Xs =x)) II (pxx)) =x)),
(11)
jÆS
where 8 is a Dirac delta function that forces the value of Xg to xs-
A.2 Relation between Causal and Conditional Shapley Values
In this section, discuss the relation between conditional (Aas et al., 2021) and causal Shapley values (Heskes et al., 2020), first showing theoretical analysis that they only coincide when there are no causal relations between the features and then providing a toy example illustrating more concretely how they differ.
A.2.1 Theoretical analysis: conditional = causal Shapley values if features are independent
We provide a theoretical analysis showing that, in general, conditional and causal Shapley values are only equivalent if features are independent of each other. We start by showing how their values functions relate.
The value function for coalition S for conditional and causal Shapley values will coincide iff P(X |xs- x5) = P(X do(Xs = x5)). As shown by Rule 2 in do-calculus (Pearl, 2009), in our causally sufficient setting this is true if Xs I X i.., they are d-separated (Pearl, 2009) by the empty separating set, in the causal graph Gxa: i.e., the true causal graph after we have remove the edges outgoing of X
In other words, in this setting, the interventional and observational conditional distribution coincide iff Xs only contains non-descendonts of Xs, i.e., variables that not descendants of any of the variables Xg. This is quite straightforward and well-known implication, but for completeness we provide short proof here:
Lemma. Let G be a causally sufficient DAG over features X (X1, (X1,...,4 Xp)- Under the causal Markon and faithfulness assumptions, P(Xg Xs = xXs) P(X3 do(Xg xxs)) for any value x's in the domain of Xs if set of variahles X C X only contains non-descendunts of the other variables 5 x X
Proof. Intuitively, if X1 E Xg is a descendant, of X2 Xg, then there exists directed path from X₂ X1, which cannot be blocked by the empty separating set, X V Xg in Gx and hence we cannot. apply Rule 2 in do-caleulus (Pearl, 2009) and these two distributions do not coincide general. If none of the X, Xx are a descendant of any X₂ Xg, then either: (i) there are paths between them, they are
16

---

Published in Transactions on Machine Learning Research (03/2026)
trivially d-separated; (ii) all paths that exist between any variable X1 E Xs and variable X₂ Xg that have edge out of 11 are cut in Gxg: so they are all trivially blocked; (iii) all paths that exist between any variable X1 E Xs and a variable X2 E Xg that have an edge into X1 have a collider, ie., X1 t- ................................ ........... .... X2, since otherwise they would be directed paths from X2 to X1. making X2 an ancestor of X1, or conversely, X1 a descendant of X2. These paths are also then trivially blocked, since they contain a collider. all paths between X5 and X5 are blocked in Gx5: then Xş and X5 are d-separated Gxg. Hence by the application of Rule 2 in do-calculus (Pearl, 2009), then these two distributions coincide.
This result shows that in general the observational and interventional distribution might not coincide, unless this condition holds for special choices of structural causal models and specifie interventional values x3s- We now use it to show that for a given S C [p], the value functions for conditional Shapley values, which we denote as (S), and causal Shapley values, which we denote as (S) coincide if Xs only contains non-descendants of the other variables Xg, or if the function f only depends on the variables Xg. More formally:
Corollary 1. Let G be causally sufficient DAG over features X = =x Xpt. Under the cousal Markov and failhfulness assumplions, for coalition S _ [p], the conditional and cousal value funclions coincide, i.e.,
Veona( EXXx,IX-==) x()] E 1 x3)] =: (eausat(S).
(12)
is one of the two following conditions hold:
1. the variables indered by this set Xs only contain non-descendants of the other variables Xg. = X Xs. This includes trivially S = 0 and S = [p] for which this always holds.
The function f depends only on the variables in Xs and not on any of the variables in X
Proof. The second condition is trivial prove, since the expected value of the function will collapse to of f(x3) for both conditional and causal Shapley values. We now focus on proving the first condition for functions f that also depend on the other variables Xg. The value functions are expectations over P(X Xs = x3s) and P(X do(Xs xt)), respectively, so it suffices to show these distributions are identical under this condition. As shown in Lemma A.2.1, this interventional and observational conditional distribution in general coincide if Xs only contains non- dendeescessssss of Xs.
In our setting this implies that for a feature sufficient condition for the conditional Shapley values (Aas et al., 2021), which we denote as écond,i, to be identical to the causal Shapley values (Heskes et al., 2020), which we denole as &eausal,;; is that the respective value functions coincide for all sets S C [p]. We show that this happens if the causal graph contains no edges. While in principle, there could be other cases in which the Shapley values might coincide, despite the individual value functions being different, we consider these as pathological cases. We formalize this result as follows:
Lemma. Let G be a causally sufficient DAG over features X = (Xi..., Xp). Under the causal Markov and faithfulness assumptions, for a fenture i, the conditional and causal Shapley values coincide, i.e.,
@oond, doausaliis
(13)
if G contoins no edges, i.e., the variables are independent.
Proof. A sufficient condition for the conditional Shapley values (Aas et al., 2021), which we denote as Ócond, to be identical to the causal Shapley values (Heskes et al., 2020), which we denote as Ócausal.i, is that the respective value functions coincide for all sets S C [p]. In this setting, the second condition of Cor. 1 cannot ever hold, e.g., since this would require f depends at the same time only on a variable X, E X and only on a different variable X, E X, which is logically impossible. So using the first condition of Cor. 1, these results would coincide if every variable is a non-descendant of all other variables, i.c., if the causal graph contains no edges.
17

---

Published in Transactions on Machine Learning Research (03/2026)
A.2.2 Example with three dependent features and linear predictor.
Heskes et al. (2020) (Section 4, Fig. 1) already provide examples of differences between marginal, conditional, asymmetric and causal Shapley values.
Here we provide slightly different example that makes the distinction between conditional and causal Shapley values more concrete also when we have Markov Equivalence Class. We focus well-known causal structure, the v-structure, for which the Markov Equivalence Class only contains one DAG, to show how this affects the MEC Shapley value methods. In this specific example, a marginal Shapley value would provide the correct output, since this example is similar to the fork example by Heskes et al. (2020), but with an additional variable. On the other hand, as discussed in detail by Heskes al. (2020), marginal Shapley values can only attribute direct effects, which might discard the contribution variable that has an indirect. effect, e.g,, in their chain example.
While in general our method allows for MECs with multiple DAGs, we focus this simple case to show how our output differs w.r.t. conditional Shapley values, even when are learning causal graph from data. Moreover, even in the case of multiple DAGs in the MEC, our approach will output a list of causal Shapley values, which would in the oracle case include the one the ground truth DAG, while conditional Shapley values will only output a single explanation.
Structural Causal Model. For simplicity, we consider the following linear Gaussian SCM with independent. standard Gaussian noise variables 6 Ez!
X1 C
(14)
X₂ C2
(15)
X3 = axx + 02X2 + E3
(16)
62, E3 (),
(17)
This SCM results in X2, (0, 1) and X3 N(O, aj + ag 1). For simplicity, we consider a linear prediction function that by design only uses Xy:
f(2,, Br1
(18)
We represent this SCM graphically in the figure below. Note that this SCM an instance of well-known
X1
X2
P
022
X₂
class of causal graphs, a v-structure (Pearl, 2009), one of the basic blocks of cansal discovery. In our sctting with causal sufficieney, the Markov Equivalence Class (MEC) of the v-structure only contains one DAG: the v-structure itself. So in this setting, there is only one MEC Shapley value, which coincides with the causal Shapley value with known graph.
Two value functions. We compare conditional and causal value functions:
Vound (S) = Ep(x,
(19)
Ucausal (S) =
(fX
(20)
Given one of these two value functions v(.), the Shapley values for features i can be written as:
1 2 -1 Σ
di(1)= = 3
(15))
eISU (i)) v(S).
(21)
SCPMM
18

---

Published in Transactions on Machine Learning Research (03/2026)
where can already simplify some terms. In particular, can show that regardless distribution P E (P() P(Jdo(X))) and corresponding value function v(.):
v00 [x X2, Xa)] v 2(1, = (flf( 23)] Bp [3x1] A
[p[3X1]
(22)
(23)
By using Cor. 1, we can further simplity the following:
Feond ((1}) = leausal (1))
= B
because
1 x Desc (2,3})
(24)
Ecoud ({2}) = Ucausal (2])
= 0
because
2 2 Desc( {1, 11,3)
(25)
Foond (1, 2}) Urausal [1, 2])
BZ
because
1,2 & Desc((3))
(26)
Ucond [1,33) Urausal (1,3))
Bai
because
I depends only on
(27)
We start by computing the part that is common to both types Shapley values by substituting the values we have alrcady computed:
1
1 1
1 + ((1, ((22,3))
(v)
[(((1)) v00
2 (1)
+ 2 2 ((1.3))
(28)
3
1
3e1 v(3))
[31
+
+B1 v((2,3))
(29)
3
1
v((3))
(38m1
v(2,3))
(30)
3
2
1
2 ((1]) 1 3)) (((33))+ o((1,3))
(2(v)=
[v((2) v00)
v ((1122)
(31)
3
[u((2,3)
v((3))
(32)
6
1
1
1
(1)
[v((3)) v0) ([L.33)
((1)) 2 ((2.3)) 3 + -
(33)
3
2
1 + (12,33))
1
[o((33)
(34)
3
We now compute the Shapley value functions that differ, v({3)) and v({2,3)). We start with the conditional Shapley value for feature 3:
Ucond ((3)) = EP(X [[XX, X2, x2)) =1 -XIxXXXX
(35)
= =[[SXi]+ + Cou(X,, Xz)
Var(Xs) (I' -((xs)
(36)
=0 a +0222 + 63)
53 +0
(37)
Var(Xs)
a1
3 aj az +1 33
(38)
where for the first equation, we first use the fact that f only depends on X1 and that X, and X₂ are independent, while then we apply the conditional expectation formula for two joint Gaussians and substitute the mean and variances structural equation for X3, where with X1. resulting in Coe(X,, Xs) = Coo(X1, X1) 01.
already computed. To compute the covariance, we substitute directly the
simplify the terms with X2 and e3 since they are both independent
We compute the causal Shapley value function, which simplifies substantially since intervening on a variable that is not an ancestor of X, has no effect on its distribution, hence simplifying to:
Vcausal ((3)) == X |f X2, xz)] (p(x)) 0.
(39)
which is substantially different from Weand (3)).
19

---

Published in Transactions on Machine Learning Research (03/2026)
We focus the conditional Shapley value for features 2. This requires a more complex computation,
based on partitioning the jointly Gaussian X1,X2 X3 in two groups: X1 and V := [X2, Xs]" with values
= [22, x3] so we can compute the conditional expectation using the standard formula for Gaussians:
Uoond (22,3)) == XXxxx|X2=222x X2 [X1, L2, 2z)] | xg
22,
[f(x)) [X
(40)
= =B + 8Ex, Ivv(o [vV[)-R[V) 0 3Ex, Evv(r 10,00²)
(41)
=B Evv
(42)
where the covariance matrix E for X1, V is composed as
Σxix, Ex,v
=
Σvx; Σvv:
We compute only the parts that we need for value function:
Σx, [Cov(X1, X2)
Xo)(XI1 = [0 a1]
(43)
Var(X2)
Cou(Xgg
a2
Evv Cov(X2, Xz)
(44)
Var(Xs)
Q2
of a +1
1
af
+ ag
+1
-02
(45)
of +1
1
By substituting all of these the value function formula, we get:
Voond ((2,3)) 3 [o
a1] af (x3 a2222),
1
a a2
+1
-0-
x22
(46)
-02
1
Ø3
3 a
(47)
We compute the causal Shapley value function for and 3, which again simplifies substantially since intervening on a variable that is not an ancestor of X2 has no effect on its distribution, hence simplifying to:
Veausad ( (2,
(oox,,do(XXo xx) [XX,,
(48)
which is again substantially different from Uoond (12,3)).
Difference in Shapley values. These changes in the value functions contribute to a change in the Shapley values. We use the simplified equations from before and we can now see clearly the difference in conditional (Aas et al., 2021) and causal Shapley values (Heskes al., 2020) for feature 1:
Urond ((3))
Φcoud1
[3300
Woond ((2,3))
8
a1
01
Bx1
D3 +
(13 2202)
3
2(a + a
+ 1)
af +
[ Urausal (33) Veausal (2,3)) 3x1,
ecansall
3
The causal Shapley value correctly attributes the explanation to the value 21 with the correct coefficient. On the other hand, the conditional Shapley value provides a biased output that can be quite far from the true causal effect, depending on the values of the datapoint and the coefficients in the SCM. Similarly, for features and get:
Øcond2
Weond (oann((1,2)) (12,3)) - Weond ((3))]
6
3 a
of
(03
(272)
of o +
a1
33
@rausal2
Veausal ([2,3))
Veausal ((3))
Øcond'3
Ucond (f3))
Dcond (2,3))
/
a1
F3
a1
(xz 2222)
3
?
+
+
2( + 1)
1 2
Oransal3
[oood ((3)) + Urond ((2,3)) -
20

---

Published in Transactions on Machine Learning Research (03/2026)
Neither of these two features has any effect the prediction, so the causal Shapley value provides the correct explanation, while the conditional Shapley value provides a biased explanation due to conditioning also on variables that are descendants of X1. In this specific example, a marginal Shapley value would provide the correct output, since this is a variation of the fork example in Heskes et al. (2020) with an additional variable. On the other hand, marginal Shapley values disregard indirect effects on the prediction, which might be an issuc many scttings which the variable that used the prediction is just a proxy of the actual variable that influences the prediction, e.g., ZIP code in credit card approval system, or other types of unresolved discrimination, as discussed by Frye et al. (2020).
A.3 Complete experiments synthetic data
Below we present the full set of results for all experiments performed. We investigate what happens when we only consider CPDAGs that were learned incorrectly, the range provided by each MEC Shapley, and ablations for the hyperparameters fmes nobsy Red: Reomb.
A.3.1 Performance on incorrectly learned CPDAGS
In Fig. 8 we report the same metrics in Fig. 2 in the main text, but now including only models for which neither PC and FGES learned the correct CPDAG, i.e. the oracle CPDAG. For each setting we report on how many models out of the original 40 are learned incorrectly: for Fig. Sa the number of models analyzed are 12, 18, and 26 for
5, 10, 15. For Fig. 8b, the numbers are 16, 32, and 36 for p = 5, 10, 15. For Fig. 8c, we have 5, 18, 30 models for d = 1,2,3, 3, Overall this ranges from a third for the simpler settings (linear, small number of nodes or sparse) to half for the default setting and almost all for the non-linear case. We see that in almost all cases results are not strongly affected by excluding the models that were learned perfectly. In almost all settings the average explanation is better than the marginal/ conditional nonginal/conditional Shapley values, and in
many cases even the worst explanation is comparable to, and sometimes better than, the baselines.
Method
Marginal
Conditionl
Oracle MEC
PC MEC
FGES MEC
Estimation Method
Sampling
0.75
0.5
NLZE
0.50
0.55
0.04
ID Number of Nades (p) (6) Linear Gauersian
10 Number of Nodes (p) (0) Non-Linesr
Number excected meichbors of (c) Lineer Gaussian
Figure 8: Normalized T.2 error between the ground truth causal Shapley value and the Shapley values returned by each method for different parameters, where only models are included where both PC and FGES learned incorrect CPDAG, i.e. a CPDAG not equivalent to the oracle CPDAG. The amount of models included is (a) 12, 18, 26 for p = 5, 10, 15 in the linear Gaussian setting, (b) 16, 32, 36 for p 5, 10, 15 in the non-lineat setting, (c) 5, 18, 30 for d = 1, 2,3 3 in the linear Gaussian setting.
A.3.2 Mean and standard error plots for Fig. 2
For completeness, in Fig. 9 we report the mean and standard error of the NL2E, averaged over the DAGs in the MEC and over the Nse = 40 causal models. In particular, the standard error is computed over Tsem = 40 causal models, after having averaged over DAGs. We report for (a) the lincar Gaussian setting where we vary p, (b) the non-linear setting with varying p, (c) the linear Gaussian setting with varying d. These results show the stability of the evaluation over the different causal models.
A.3.3 More flexible conditional estimators
To assess whether the performance gap between our method and the baselines is simply due to the misspecifi- cation of the Gaussian conditional estimator in the non-linear setting, we conduct an additional ablation study.
21

---

Published in Transactions on Machine Learning Research (03/2026)
Method
Marginal
Estimation
Sampling
0.
0.8
0.195
NL2E
0.100
0.6
15
15
Number nodes
Number of noces p (b)
Expected number of neighbors
Figure 9: Mean and standard error for NL2E between the ground truth causal Shapley value and the Shapley values returned by each method for different parameters and settings: (a) linear Gaussian setting with varying number of features p, (b) nonlinear Gaussian setting with varying number of features p, (c) linear Gaussian setting with varying number of expected neighbors d. For all MEC Shapley methods we plot the average NT.2E over all the explanations, with error bars representing the SE over neem = 40 causal models.
We extend the Conditional Shapley baseline by employing two more flexible, non-parametric conditional density estimators provided by the shapr package: Gaussian Copula and Conditional Inference Trees (ctree). Figure 10 reports these results. Note that while enhance the baselines with these flexible models to capture non-linear dependencies, the MEC Shapley variants continue to use the simpler Gaussian conditional modcling. This provides a conservative comparison, testing whether causal structure lcarning (even with misspecified parametric assumptions) yields better explanations than sophisticated conditional modeling that ignores causality.
Because the true data data-generating mechanism is non-lincan Gaussian conditionals are misspecified. Nevertheless, as shown in Fig. 10, MEC Shapley values remain systematically closer to the ground-truth causal Shapley values (computed from the true interventional distributions) than conditional Shapley values. Using more flexible conditional estimators narrows the gap to the causal target for the conditional baseline, but it does not remove the advantage of MEC Shapley: on average, MEC Shapley values are still closer to the causal Shapley values, even when MFC Shapley is computed under the misspecified Gaussian conditional model.
Marginal
Conditional (gaussian)
Conditional (copuls)
Conditional (ctree)
Method
Oracle MEC
PCC MEC
FGES
Estimation Method
Sampling
IW
0.12
0.09
0.06
10 Number of nodes Non-linear
15
Figure 10: Ablation study comparing Conditional Shapley values with flexible non-linear estimators against our MEC Shapley framework in the non-linear setting. We report the Normalized T.2 Error (NI.2F) for Conditional Shapley values estimated via Gaussian assumptions (standard), Gaussian Copula, and Conditional Inference Trees (ctree). Even when the Conditional baseline is enhanced with flexible estimators to capture non-linearities, the MEC Shapley methods (using standard Gaussian conditionals) consistently provide explanations closer to the ground-truth causal Shapley values.
22

---

Published in Transactions on Machine Learning Research (03/2026)
Method
Marginal
Conditional
Oracle MEC
PC MEC
FGES MEC
Estimation Method
Sampling
IW
0.12
0.12
0.09
0:09
0.06
0.06
2
6
Number of expected neighbors
Number of expected parents of Y
(a)
(b)
Figure 11: NL2E as function of (a) the expected number of neighbors of cach node d and (b) the expected
number of parents of Y in the non-linear setting.
A.3.4 Varying the amount of neighbors and parents
In Fig. 11 we vary (a) the expected number of neighbors for each node and (b) the number of expected parents of m, both in the non-linear setting. We can see that the number expected neighbors has a strong influence on causal discovery performance the non-linear setting, as the MRC Shapley values for PC get progressively worse. The difference between the range resulting from the oracle MEC and from the PC MEC shows that the true graph is average not contained in the MEC, and the best graph average leads to explanations deviating far from the one of the true graph. Increasing the expected number of parents of Y m seems to affect all methods equally, with the results getting worse as m increases. Here we fix the amount of nodes d = 10. We use the default settings Red causal Shapley values.
Tome
Robs 1000, with T' 4000 for the ground truth
A.3.5 Amount of Monte Carlo samples
Method
Marginal
Conditional
Oracle MEC
PC MEC
FGES MEC
Estimation Method
Sampling
TW
0.125
NL
0.4
0.100
0.075
0.050
250
500
1000
2000
4000
250
500
1000
2000
4000
Dmc
Dmc
(a) Linear Gaussian
(b) Non-linear
Figure 12: NL2E between the ground truth causal Shapley values and other methods in (a) the linear Gaussian setting and (b) non-linear setting as a function of number of Monte Carlo samples R'm used estimate cach expectation.
In Fig. 12 we perform an ablation where we vary Ruce the amount of Monte Carlo samples used to estimate expectations for (a) the linear Gaussian setting and (b) the non-linear setting. We can observe that the results in terms of error relative to the ground truth are, after averaging over all datapoints, not significantly influenced by the number of Monte Carlo samples ae- Here we fix nobs = 1000. For the ground truth causal Shapley values we used nue 4000. In this ablation the graphs contain p 10 nodes with d 2 and 6. We the full set of combinations S.
23

---

Published in Transactions on Machine Learning Research (03/2026)
A.3.6
Amount of data for estimation of distributions
Wethod
Marginal
Conditional
Oracle MEC
PC MEC
FGES MEC
Estimation M Method
Sampling
T
0.16
0.6
0.12
NL
NL2E
0.4
0.08
0.2
0.04
0.0
250
500
1000
2000
4000
250
500
1000
2000
4000
nobs
nobs
(a) Linear Gaussian
(b) Non-linear
Figure 13: NL2E between the ground truth causal Shapley values and other methods in (a) the lincar Gaussian setting and (b) non-linear setting as a function of number of observational samples Yoass used to estimate each marginal or conditional distribution.
Fig. 13 shows the average NL2E as function of Robs: the number of observational samples used to estimate marginal and conditional distributions, for (a) the linear Gaussian setting and (b) the non-lincar setting. Here we can see that the amount of samples used to estimate distributions affects all methods equally. A low amount of samples increases the error for all explanations, but we can see that as Robs increases, this effect starts to diminish. Interestingly, this holds also for the non-linear case, where the parametric assumptions for the conditional and MEC Shaplcy values fails hold.
A.3.7 Amount of combinations
In Fig. 14 we investigate what occurs when we use only a subset of [p]. Unsurprisingly, the error relative to the ground truth decreases as we use more combinations. We can see that the differences between methods stay consistent, with the average error decrcasing with the amount of combinations. Here we also include the ground truth computed using only a subset of the power set of [p], denoted Ground truth (limited), which uses the same number of combinations Beamb as the other methods. We can that the between the methods and the ground truth computed using a lower number of combinations stays consistent, which is important given that we do not use the full power set of subsets for the results on 15 nodes.
A.3.8 Amount of data for causal discovery
In Fig. 15 we vary the amount of samples available for causal discovery for the linear Gaussian setting (a) and the non-linear setting (b). We can that as the number of samples grows, the difference between the Oracle MEC and PC and FGES decreases, explained by the fact that the outcome of the causal discovery algorithms gets closer to the Oracle MEC.
A.3.9 Effect of number MC samples on computation time
In Fig. 16 we plot the computation time for single datapoint for each method. We show how Rom affects computation time, where report both the time taken including the identification of equivalent interventional distributions and the time taken just for computation of the Shapley values denoted as (w/o id). While the time taken to identify equivalent distributions is significant, once performed it markedly decreases computation time. The IW method also speeds up computation time. To compare MECs across graphs, we divide total computation time by the amount of graphs in cach MEC. As expected the higher the n the higher the average running time for the MEC Shapley methods and the Conditional Shapley values, with the running time for higher numbers of MC samples comparable between Conditional Shapley and the MEC methods.
24

---

Published in Transactions on Machine Learning Research (03/2026)
Marginal
Oracle MEC
FGES MEC
Estimation Method
Sampling
IW
Method
Conditionll
PC
Ground truth (limited)
0.15
0.6
NL2E
0.10
0.4
0.2
0.05
0.0
0.00
256
512
1024
128
1024
Number of combinations Ncomb (a) Linear Gaussian
Number of combinations ncamb
(b) Non-linear
Figure 14: NL2E between the ground truth causal Shapley values and other methods in (a) the linear Gaussian setting and (b) non-linear setting as a function of number of combinations nccomb used to estimate each Shapley value.
Method
Marginal
Conditional
Oracle MEC
PC MEC
FGES MEC
Estimation Method
Sampling
TW
0.6
0.125
NI ZE
NL2E
0.4
0.100
0.075
0.2
250
1000
2000
4000
250
500
1000
2000
4000
Red
Red
(a) Linear Gaussian
(b) Non-linear
Figure 15: NL2E between the ground truth causal Shapley values and other methods in (a) the linear Gaussian setting and (b) non-linear setting function of number of causal discovery samples ned used to estimate each MEC. Marginal, conditional and Oracle Shapley values are not affected by this parameter.
A.3.10 Analysis of MECs
In Table 1 we report the average number of DAGs in each MEC as well as the standard deviation for the linear Gaussian and non-linear setting, where we vary the number of nodes In Table we report the same, now varying d, the expected number of neighbors for each node. We do not report these tables for the expected number of parents of Y m, as is included in the causal discovery phase and as such this parameter does not influence the size of each MEC. We can see that the size of the MECs varies substantially across the different models within the same setting.
25

---

Published in Transactions on Machine Learning Research (03/2026)
Method
2.0
Marginal
Average timee (min)
Conditional
1.5
FGES
FGES (w/o id)
1.0
Estimation Method
0.5
Sampling
IW
1000
2000
4000
Dmc
Figure 16: The average runtime for the linear Gaussian setting as function of the number of MC samples Tyme
Nr. of Nodes
Type
MEC size
FGES
3.00 + 2.03
Oracle
4.9 + 5.33
6
PC
3.15 + 2.25
11
FGES
7.82 + 8.08
11
Oracle
10.2 + 13.0
11
PC
7.42 + 8.20
FGES
15.2 20.3
Oracle
17.9 14.8
16
PC
13.9 10.5
Nr. of Nodes p
Type
MEC size
FGES
4.55 + 5.10
Oracle
4.9 + 5.33
PC
2.7 = 2.40
FGES
5.5 4.11
Oracle
10.2 13.0
11
PC
3.18 4.25
16
FGES
11.2 = 12.0
16
Oracle
17.9 = 14.8
16
PC 4.97 + 5.48
Table 1: Tables reporting MEC sizes for the lincar Gaussian sctting (left) and non-lincar sctting (right) where the number of nodes p varies. We report mean and standard deviation for each setting for the Oracle, FGES, and PC MECs.
In Tables and we report the average Structural Hamming Distance (SHD) for each learned CPDAG with respect to the Oracle CPDAG well the standard deviation for the linear Gaussian and non-linear setting, In Table 3 we report these values for varying numbers of nodes d, and in Table 4 we vary the expected number of neighbors mn for cach node.
Fig. 17 relates structural discovery error to explanation error for the PC-based sampling variant for Nsmm = 40 causal models, On the x-axis we report the SHD between the learned CPDAG and the oracle CPDAG for each run; on the y-axis we report NL2E. The labels on cach point represent the number of CPDAGs with that SHD. For fixed learned CPDAG, our method produces one explanation per DAG in its Markov equivalence class (MRC). Accordingly, for each run compute NL.2 for every DAG in the learned MEC and summarize
26

---

Published in Transactions on Machine Learning Research (03/2026)
Nr. of Neighbors d
Type
MEC Size
Nr. of Neighbors d
Type
MEC Size
1
FGES
7.35 + 5.26
FGES
6.78 + 4.10
1
Oracle
7.15 + 5.30
Oracle
7.15 + 5.30
1
PC
6.92 + 5.23
PC
5.12 + 5.41
2
FGES
7.82 + 8.08
FGES
5.5 + 4.11
2
Oracle
10.2 13.0
Oracle
10.2 = 13.0
PC
7.42 8.20
PC
3.18 + 4.25
3
FGES
6.72 13.0
FGES
6.3 + 8.66
3
Oracle
8.82
14.2
Oracle
8.82 + 14.2
3
PC 3.92 = 2.67
PC 2.53 = 3.52
Table 2: Tables of MEC sizes for the linear Gaussian setting (left) and non-linear setting (right) where the number of neighbors d varies. We report mean and standard deviation for each setting for the Oracle, FGES, and PC MECs.
Nr. of Nodes
Type
SHD
Nr. of Nodes
Type
SHD
FGES
1.73 2.75
6
FGES
1.58 2.04
6
PC
1.58 + 2.10
6
PC
2.4 + 2.07
11
FGES
2.00 = 3.66
11
FGES
3.85 = 3.08
11
PC
4.03 + 3.50
11
PC
5.65 = 3.04
16
FGES
4.72 = 7.17
16
FGES
5.65 = 4.41
16
PC
4.62 + 4.26
16
PC 8.48 3.72
Table 3: Tables SHD (Structural Hamming Distance) with respect to the Oracle CPDAG for linear Gaussian setting (left) and non-linear setting (right). Mean and standard deviations are reported for FGES and PC CPDAGs across varying numbers of nodes p.
the resulting set by its Wit minimum (best-case), mean (average-case), and maximum (worst-case). The three curves plot these summaries after averaging over runs that atlain the same CPDAG-SHD value, and the shaded band spans the corresponding
worst envelope. Overall, the mean NL2E increases with CPDAG-SHD, indicating that larger discovery errors tend to degrade explanation quality. When the learned MEC contains only a single DAG, the best/mean worst summarris coincide and the envelope collapses.
Interestingly, the NL2E spread across the DAGs in each MEC does seem increase predictably with higher SHD. One explanation is that the spread is completely due to the size of the MEC, which is known to be asymptotically constant and empirically on average containing 4 DAGs (Gillispie & Perlman, 2001).
A.4 Complete experimental results real world data
A.4.1 Alzheimer's data set
Following Heskes et al. (2020), we consider the Alzheimer's disease data set obtained from the Alzheimer's Discase Neuroimaging Initiative (ADNT) database (http://adni, looii. usc. edu), which focuses on identifying biomarkers for early Alzheimer's detection and progression tracking. Our analysis incorporates the same features in (Shen et al,, 2020; Heskes et al., 2020), age (AGE), education level (EDU), gender (SEX), amyloid beta (ABETA), fudeoxyglucoss (FDG), phosphorylated tau (PTAU), and the number of apolipoprotein alleies (APOE41, APOE42), for total of 8 fcatures). As ground truth for the causal discovery we employ the gold standard causal graph from Shen et al. (2020), which includes dummy-encoded variables for the apolipoprotein alleles (APOEŞ1 and APOE42), and which we show in Fig. 18 (left). For causal discovery we discretized the continuous variables with a maximum of values for each variable using the Freodman-Diacom rule. Values that are encoded as being larger or smaller than cutoff value, such as ABETA 1700, are re-encoded their cutoff value. Like (Heskes et al., 2020), group together patients with mild cognitive
27

---

Published in Transactions on Machine Learning Research (03/2026)
Nr. of Neighbors
Type
SHD
Nr. of Neighbors
Type
SHD
1
FGES
0.6 + 0.98
1
FGES
1.27 + 1.65
1
PC
0.7 + 1.12
1
PC
2.6 + 2.53
2
FGES
2.0 + 3.66
2
FGES
3.85 + 3.08
2
PC
4. - 3.50
2
PC
5.65 - 3.04
3
FGES
9.5 + 9.46
3
FGES
7.38 = 5.29
3
PC
9.3 + 5.16
3
PC
9.1 + 3.74
Table 4: Tables of SHD (Structural Hamming Distance) with respect to the Oracle CPDAG for lincar Gaussian setting (left) and non-linear setting (right). Mean and standard deviations reported for FGES and PC CPDAGs across varying expected numbers of neighbors per node d.
NL2E vs SHD (PC, Sample) Label nr. of SCMs in bin
0.9
Curve
M
0.6
Min
NI
Mean
0.3
0.0
5
10
SHD(leamed CPDAG, oracle CPDAG)
Figure 17: NL2E vs. CPDAG-SHD (PC, sampling, Nsam 40) for p 10, 2, with lincar-Gaussial relations. On average the causal graphs have 20 edges. For each run, SHD is computed between the learned and oracle CPDAGs. We evaluate NL2E for all DAGs in the learned MEC and report the WiC min/mean/ hom/mean/maa /mx (green/bluu /b/ee/eele.... Curves average these summaries over learned CPDAGs with the same SHD (the labels on each point show how many cansal models have this SHD); the shaded band spans the averaged min max range.
impairment and Alzheimer's disease, dropping patients with the diagnosis of significant memory concerns. Rows containing missing values were removed. We are left with a total of N 1500 datapoints.
After these preprocessing steps, we apply the PC algorithm on a random subset of 80% of the data (0.8 1500 1200 datapoints) with the G-square test with 0.05, and the FGES algorithm using the discrete BIC score with penalty discount À 1. We prohibit edges from biomarkers into demographic variables following (Shen et al., 2020). The TETRAD library used to run both algorithms (Ramsey et al., 2018). In Fig. 18 on the right we show the Iearned CPDAGs, with blue edges identified by both FGES and PC, green edges exclusively by FGES, and red edges only by PC. The blue and green edges agree with the gold standard graph, whereas the red edge (FDG + ABETA) is oriented incorrectly. Compared to FGES, the CPDAG learned by PC is missing the edge from ABETA to PTAU and it inverts the edge from ABETA FDG. Neither of the methods include the edge AGE АBEТА. The two MECs represented by the CPDAGs each contain only single DAG.
28

---

Published in Transactions on Machine Learning Research (03/2026)
AP0241
AMO542
APOBAT
APOBAZ
PGRS
FDG
AGE
AGE
PGES
ABE
ARE
CDJ
Figure 18: Comparison of gold standard graph recreated from Shen et al. (2020) (left) with CPDAGs discovered by FGES and PC algorithms (right). Edges are color-coded to indicate discovery: blue for edges identified by both FGES and PC, green exclusively by FGES, and red only by PC. The blue and green edges agree with the gold standard graph, whereas the red edge is oriented incorrectly.
Like (Heskes et al., 2020), we consider a binary classification task, grouping together patients with mild cognitive impairment and Alzheimer's discase, dropping patients with the diagnosis of significant memory concerns. We train an XGBoost model (Chen & Guestrin, 2016) for 100 rounds on the subset used for causal discovery (0.8 x 1500 = 1200 dalapoints), with the other 20% unseen (0.2 x 1500 = 300 dalapoints) by the prediction model and used for evaluation. In the sina plots in Fig. 19 we show for each feature the marginal, conditional, gold-standard causal, FGES, and PC Shapley values, respectivcly. We can see that FGES and ground truth Causal Shapley values both capture the indirect ellect of APOE41 on the prediction, while the marginal Shapley values do not, as well as attributing a larger importance to APOF42.
Ackuowledgments. ADNI data collection and sharing was funded by the Alhhemer's Disease Neuroimag- ing Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W8IXWI1-12-2-0012); See http://adni .oni woooiiusc.eeeeooooooooeoooeeeeee www.uuc.ccc.... fett/wp-content, Datpr//rrAAADDDaaaaaaaaaAAaAAae butions and http://adni. lonii usc www.uuu..uuuuuuuuuuuu List.pdf for a complete listing of ADNI investigators.
Data_ Use Uawe_A_rrementtpppeeee pdf, section 12, for further contri- to Apppa//AAADAAAACoollleeeeeee
A.4.2 Bike rental dataset
We consider the bike rental dataset of Fanace-T & Gama (2014), with the same features and preprocessing as used by Heskes et al. (2020): the number of days since January 2011 (trend), two cyclical variables to represent scason (cosyear, sinyear), the temperature (temp), feeling temperature (atemp), wind speed (windspeed), and humidity (lum), for a total of 7 continuous fral We use the partial order ({Irend), &cosyear, sinyear), ) (all weather variables)) ) from Heskes et al. (2020) as background information for the PC and FGES algorithm, both ran using Tetrad. For PC we apply Fisher's Z test with a 0.05, and for FGES the BIC criterion with x = 1. The resulting CPDAGs are shown in Fig. 20. In this dataset we do not have the ground truth causal graph and moreover, there are likely many latent confounders, so we can only reason qualitatively. The CPDAG discovered by PC is sparser than the one by FGES, missing the link between trend and temperature, as well as temperature and windspeed. Moreover the two CPDAGs disagree on the direction of the edge between windspeed and the feeling temperature alemp. The CPDAG discovered by FGES contains a single undirected odge between humidity and windspeed, while this edge is oriented windspeed -> hum in the CPDAG discovered by PC. In both graphs, some of the edges that we intutively expect are missing, e.g., trend + cosyear, possibly because sinyear and cosyear are deterministic functions of the trend.
We split the data (total number of N = 730 datapoints) into 80% (0.8 x 730 = 584 datapoints) training data and 20% (0.2 x 730 = 146 datapoints) test data, and trained an XGBoost model for 100 rounds, predicting the amount of bike rentals based on the previously mentioned features. In Fig 21 we show sina plots for the marginal Shapley values, the conditional Shapley values and the PC Shapley values (right). We see that the PC Shapley values better capture the contribution of cosyear than marginal and conditional Shapley values, which represents the seasonality of the bike rental (Le. higher bike rental in summer than in winter).
29

---

Published in Transactions on Machine Learning Research (03/2026)
Scaled feature value Low
High
Sca led feature value Low
High
AGE
AGE
SEX
SEX
EDU
EDU
FDG
FDG
ABETA
ABETA
PTAU
PTAU
APDE41
APOE41
APOE42
APOE42
-0.50
-0.25
0.00
0.25
0.50
-0.50
-0.25
0.00
0.25
0.50
Marginal Shapley value
Conditional Shapley value
(a)
(b)
Scaled festure wabe
1899
Scaled feature value
High
Scaled feature vahre in
462
30E
323
EDU
ED0
100
AB:TH
K0ETS
AB1TA
P7U0
P1AU
AP0141
AP0E4.
AP0E4L
APO2/2
APDE42
AP
-450
-0:25
1.00
0.25
0.22
-0.5
-022
0.22
025
1.0
-0.0
-0.25
0:00
0.2
350
Gold standard causal Shapley value
FGES MBC causal Shapley value
PC MEC causal Shapley vabe
(c)
(d)
Figure 19: Sina plots for a binary classification task using XGBoost on the Alzheimer's dataset. We plot (a) marginal, (b) conditional, (c) the causal Shapley values for the 'gold standard' DAG, (d) FGES MEC Shapley values and (e) PC MEC Shapley values.
FGES
PC
trend
trend
cosyear
sinyear
cosyear
sinyear
temp
temp
atemp
hum
atemp
hum
windspeed
windspeed
Figure 20: The learned CPDAGs on the bike rental dataset for FGES (left) and the PC algorithm (right).
The figure also shows that there is little difference in the altribution between lemperal ure and perceived temperature when using conditional Shapley values, unlike with PC Shapley values where the difference is more noticeable. In Fig. 22 we show the sina plots for the two DAGs in the CPDAG discovered by FGES. We can observe that when windspeed causes humidity (a), the Shapley values for windspeed are slightly more spread out than when humidity causes windspeed (b). Compared to both FGES plots, in the sina plot of PC
30

---

Published in Transactions on Machine Learning Research (03/2026)
Scaled feature vOMC
High
Scaled fealure VOlU Low
High
Scaled feature volup Low
100d
200ye1
32y031
tomp
tamp
3010
w0p0poo
widggood
win0000d
-000
-1000
1000
2000
-200
-1000
1000
2000
-2000
-1000
1000
2000
Marghal Shaplyy value
Conditonal Shapley valuo
PC Shaplyy valuo
Figure 21: Sina plots for regression task using XGBoost on the bike rental dataset. We plot marginal (left), conditional (middle) and PC Shapley values (right).
Scaled feature value Low
High
Scaled feature value
gh
sond
bond
cosyear
cosyear
sinyea
temp
atomp
atomp
windspeed
windseeed
hum
hum
-2000
-1000
1000
2000
-2000
-1000
1000
2000
FGES S hapley value (windspeed ->i hum)
FGES Shapley value (hum - windspeed)
(a)
(b)
Figure 22: Sina plots for a regression task using XGBoost on the bike rental dataset. We plot FGES shapley values for the two DAGs in the CPDAG discovered by FGES (shown in Fig. 20 on the left).
we can observe slightly more spread in temperature, which is consistent with trend being an additional cause in the FGES CPDAG. In Fig. 7 we show the bar plot of the Shapley values for two days with temperature about 13 degrees Celsius in October and December. While this temperature represents an average day in October, it is a warm day for December, hence the contribution of the temperature to the prediction of bike rental is substantially more positive for December for both Causal chain and PC Shapley values. Similarly, there is a significant seasonal drop in bike rentals from October to December, which can be seen in both Causal Chain and PC Shapley values, while this does not hold true for marginal Shapley values. Finally, similarly to what is shown in the sina plot, PC Shapley values attribute a bigger impact to temperature than Causal chain Shapley values.
A.5 Implementation and compute
To implement our method, we adapted the 'shapr' package by Aas et al. (2021) (https://github commgh.hh.o.//// NorssReggessntrall/hhaap relcased under the MIT license. The causal discovery methods we used were implemented in the 'pealg' package (Kalisch et al., 2012) released under GPL-2, the 'kpealg' package (Zhang et al., 2011) released under GPL,2, and the Tetrad package (Ramsey et al., 2018) released under the GPI.2
31

---

Published in Transactions on Machine Learning Research (03/2026)
date
2012-10-09
2012-12-03
Causal (chain)
Conditional
Marginal
PC
Shapley value
-500-
cosyear
temp
cosyear
temp
cosyear
temp
cosyear
temp
Figure 23: Bar plot of Shapley values for the features cosyear and temp of the bike rental dataset for two days October and December with the same temperature. We show the Causal chain Shapley values from Heskes al. (2020), along with marginal, conditional, and PC Shapley values. Both FGES Shapley values are similar to PC and have been omitted.
license. We make use the 'igraph' package (Csardi Nepusz, 2006), released under the GPL license, to represent graphs. To identify causal effects we make use of the 'causaleffect' package (Tikka & Karvanen, 2017a) released under GPT.-2 The prediction models trained were xgboost models (Chen & Guestrin, 2016) released under the Apache License v2.0. For estimating conditional Gaussian distributions of the 'mvnfast' package (Fasiolo, 2014), released under GPL-2. For analysis of results and the creation of plots we make use the 'tidyr', 'dplyr' and 'ggplot2' packages, all part of the tidyverse (Wickham et al., 2019), released under the MIT license. We provide the code to reproduce all our experiments in the supplementary material.
make use
A.5.1 Hyperparameters
In this section we list each hyperparameter used in the experiments, along with an explanation, how we vary them, and their default value in bold.
When sampling random DAGs, specify p the amount of nodes excluding the dependent variable (Y), d the expected number of neighbors for each node X; (representing features), and m, the expected amount of parents of Y, i.e. the expected number of features directly affecting the dependent variable. We use -= 5, 10, 15, d = 1, 2,3, and m = 3, 6,9 for p= 10.
We consider two parameterizations for each graph: linear Gaussian and non-linear. In the linear Gaussian setting, nodes are parametrized as X, Ejera, WiyXi + 6, where wij U((-2, -0.5) u (0.5, 2)) and 6, 10, 1).
For the non-linear setting consider parameteriiattons of the form X; w (E EjePa, w ) + 4 where a denotes the sigmoid function, w ~ U((11.5, -0.5) U (0.5, 1.5)), wh ~ U((33, 1) (1,3)), and e; follows either a (0, 1) or a U(-1,1) 1) distribution with equal probability. As such, for each DAG we get two different. SCMs, one linear Gaussian and one non-linear. For each and each parameterization we sample 7sem = 40 random SCMs.
For each setting of hyperparameters generate Tsse 40 random causal models.
32

---

Published in Transactions on Machine Learning Research (03/2026)
For causal discovery, use the PC method implemented in the 'pcalg' package (Kalisch et al., 2012). We set a default significance level of = 0.05. We use the default settings, which uses PC-stable by Colombo & Maathuis (2014). For the linear Gaussian parameterizations we use partial correlation conditional independence tests ('gaussCT) and for the non-linear setting kernel-based conditional independence tests implemented by the 'kpcalg' package (Zhang et. al., 2011); specifically the Hilbert-Schmide Independence Criterion gamma test. For FGES we use the TETRAD toolbox (Ramsey et al., 2018) using the BIC score with default penalty discount 1. We use the following values for number of observational samples that we apply PC and FGES to ned 250, 500, 1000, 2000, 4000.
Prediction models were trained using XGBoost (Chen & Guestrin, 2016) using 7 0.1, and 'subsample 0.8º for 100 rounds. We used Thrain = 10K observational samples to train each model.
For computing Shapley values we consider the hyperparameters Noomb: Rme: Robs ccmmb specifies the number of combinations used by all methods, where we use the full set of combinations for p 5, 10, and for p Peomb 8192 for the ground truth and Neooh 4096 for the other methods. Nmm specifies the number of Monte Carlo samples used to estimate each expectation. We use the values Rmc = 250, 500, 1000, 2000, 4000, with thee = 4000 fixed for the ground truth. Robs specifies the amount of observational datapoints used to estimate marginal and conditional distributions, which we vary Robs = 250, 500, 1000, 2000, 4000.
For each SCM we evaluate performance on Aetstt = 40 data points unseen by the prediction model.
A.5.2 Computing infrastructure
The full experiment suite was run in parallel on cluster of AMD Rome CPUs with 16 cores each. Total compute time for all experiments listed was approximately 88 days. Individual runs in the default setting per 40 data points varied from 30 seconds (marginal, conditional), up to 4 hours for MECs of the largest size. Runs for 15 nodes could take up to 8 hours for 40 datapoints.
33


Reviews:
[{"review_id": "N2eUA1mA2H", "reviewer": "Reviewer_EMdU", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "First of all, as a reviewer within the realm of causal AI, it is important to thank the authors for their cool project. Any time invested into the philosophy, science and engineering to make systems explainable in a causally consistent manner that can actually be deployed widely is awesome and much appreciated.\n\nThe authors present in their manuscript a complete pipeline that goes from data (purely associational) to learning different possible causal structures to then give causal explanations based on those all while respecting the fact that this requires re-use of computation to efficiently get said explanations (since we observe [super]exponential behavior for both DAG possibilities and MEC size). A derivation for the estimation methods is provided. Extensive experimentation is provided covering different data sets, settings and ablating over important parameters. An appendix with additional discussion is provided. Further code for reproduction is provided (not evaluated).", "claims_and_evidence_yn": "Yes", "claims_explanation": "A key challenge for causal Shapely values, which poses the key method the authors build upon, is the requirement of the (true) causal graph. This is solved by simply attaching a causal discovery process a priori to the pipeline that learns from the same data the predictor-to-be-explained was trained on. This naturally comes with challenges of both design choices regarding performance (accuracy, time and such). The authors focus on a robust pipeline that re-uses computation to be efficient. The authors present two estimation methods: 1) using sampling where interventional data is generated using a candidate graph given the TF-formula, and 2) using IW where ID is used with the graph to get an estimand applicable to the available data. The key technical contribution is the grouping of graphs into sets that share interventional distributions and said estimation methods for the interventional expectation in Eq.4. While said contribution might be seen as only slightly incremental, the execution both technically and presentation-wise are solid. Furthermore, a proper study is presented that gives a reproducible account of our topic at hand. Reproducibility, precise language, interesting findings, important problem, mechanistic explanations - in terms of scientific endeavour everything checks out.", "audience_yn": "Yes", "audience_explanation": "My judgment on this is that likely anyone interested in XAI or causal AI will find this paper interesting and might even consider using it in their paper.", "requested_changes": "IMHO the main changes to this particular submission are of polishing nature. One could go back on any segment of the paper and make sure that they are to be understood by as many readers as possible for instance being self-contained e.g. in Alg.1 when reading something like $p$ one might have forgotten by that point that it was the symbol to refer to the total number of features or the \"final\" feature of our data. Or polish the key figure such that $S_i$ becomes more clear when trying to understand (particularly for the first read) the grouping of the graphs. Things of that nature essentially, which there are many of that one can come up with to make the paper \"perfect\".", "broader_impact_concerns": "Explanations and sound reasoning are crucial in our current age of AI. IMHO this type of research will benefit all of us greatly down the road. Other than that the paper does not raise any concerns. My hope lies in the continuation of such research.", "additional_comments": "Ideally the authors can extend their experimentation. More data sets, more CD methods.", "extra_scores": {}, "extra_text": {}}, {"review_id": "ZkdlaOqVAn", "reviewer": "Reviewer_z4iG", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This manuscript looks at the potential of using causal discovery to improve explainability analysis through shapley values.  In particular, they focus on the case where the causal graph is not known a priori and needs to be estimated from data, and how that can be used to modify the structure put into estimating the shapley values.\n\nOverall, the manuscript is clearly written with a clear idea.  The algorithm is straightforward to implement although computationally challenging in larger dimensions.  There is a clear advantage over the chosen baselines in non-linear scenarios, although advantages in linear scenarios are less clear.\n\nThe biggest weakness to me is that it is not clear when this would provide substantially different or better results than correcting for the relationships between the variables; e.g., do we need the causal viewpoint, or do we just need to estimate relationships between variables better?  The chosen baselines are fairly minimal, considering standard (independent) shapley analysis and correcting for a linear structure--or so I assume, as the algorithmic details on the baseline are not particularly precise.  I would be more convinced by the approach is there was clear theoretical analysis about when the results will diverge (e.g., are they essentially the same in the linear gaussian case?) and if there was greater attempts to capture non-linear relationships between variables in a baseline model without causality.  This, to me, is addressable.\n\nThe other weakness is the challenges of scalability; the primary shown results are based on relatively small and not particularly connected graphs with a large number of measurements.  While that certainly covers a subset of problems, the majority of ongoing scientific challenges that I interact with have fewer measurements and more variables, making it unclear how broadly applicable this method is.  I would encourage the authors to consider scaling to a greater extent, which is challenging with the causal discovery approaches.", "claims_and_evidence_yn": "No", "claims_explanation": "The first gap in the evidence that I see is that that there is not a clear enough evidence about when this approach would outperform a conditional method in a linear-gaussian case, as the results are largely overlapping.  I would appreciate this being examined from both a theoretical and practical perspective.  In particular, is it just that the two-step method (causal discovery + interventional effects) has a more robust estimation of the correlation structure?\n\nThe second gap is that there is not a clear enough description of the baseline conditional approach and how it is implemented.  There are several different versions of this approach, so details are necessary.  It is somewhat implied that it is the linear formulation, which would make sense as for why it does not capture non-linear relationships.  If so, then it should be considered what happens when you capture non-linear conditional relationships.\n\nAdditionally, it would be helpful to visualize the Sina plots in the main draft as to show how the estimates change.  For example, Figure 16 is helpful for understanding the differences between the methods; however, there the primary difference is that it captures indirect effects.  Given that the marginal method cannot capture indirect effects, it would be informative to see what the conditional results show differently there.  That would help to give a sense of the practical importance of the method.", "audience_yn": "Yes", "audience_explanation": "Shapley analysis is a commonly used method, and improvements in its efficacy would be appreciated.", "requested_changes": "I would like the authors to:\n(a) evaluate mathematically when the proposed approach differs from a conditional estimation.\n(b) provide more details on the baseline algorithms (how were the conditional relationships estimated)?\n(c) evaluate how different conditional estimation procedures impacts results\n(d) visualize shapley values to compare against all methods\n(e) ensure all details are included on data sets (e.g., sample numbers, etc., on the real-world data)\n\nI think these are all necessary.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "dXZGAZr6Kb", "reviewer": "Reviewer_Su8y", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper focuses on overcoming the limitation of existing causal shapley value methods that leverage only known causal relations. The main contributions of the paper are as follows:\n- The authors propose an efficient framework that combines causal discovery methods with the shapley value methods, to yield a set of explanations for the datapoints. \n- The paper presents two variants of the framework: sampling based and importance weighting based.\n- The experiments on both simulated and real-world datasets are comprehensive covering variations like varying degree in the graph, number of parents of the prediction variable and ablation studies for the hyper-parameters. The results demonstrate the effectiveness of the method in estimating the shapley value as compared to the conditional and marginal shapley value baselines.\n\n**Key Strengths:**\n\n- The paper is well-written and clearly lays down the algorithm and mathematics behind the proposed algorithm.\n- The experiments are technically sound and comprehensive. The time analysis for the different methods is interesting.\n\n**Key Weaknesses:**\n- The experiment methodology is unclear at certain points. It is unclear what the bars in the graphs for different methods indicate. Please see below for more details.", "claims_and_evidence_yn": "No", "claims_explanation": "I have concerns about the claim that the proposed method can still output largely close feature shapley values to the ground truth/oracle values despite the incorrectly learned graph.\n\nIn synthetic cases, this claim is supported by the numbers of the proposed method. However, on the real world dataset, the PC MEC method performs quite worse and this behaviour is attributed to the incorrectly learned edge between two variables. The two observations don’t align together as if the model were robust to minor mistakes in the graph, it would be reflected in the real world case as well.\n\nTheoretically, I would expect the model to give abnormal values to the incorrect graphs depending on the SHD rather than covering the underlying mistake with close shapley values. These abnormal values could then help in identifying incorrect graphs and even localize to what relations might be oriented wrong. A discussion from the authors on this would be helpful.", "audience_yn": "Yes", "audience_explanation": "The idea of combining shapley value estimation with causal discovery is very natural and helpful. The authors also present initial ideas for tackling the computational complexity of the problem as well. Therefore, it could be a good starting point for the people working on explainability in diverse real world disciplines.", "requested_changes": "There are 2 adjustments which I would like to propose:\n- A description of the graphs is necessary, especially where the spread is coming from. If the spread is coming from the different graphs in the MEC, then why would the oracle MEC have the spread. A clear description of what the plotted values are and what the bars are will make it more clear.\n- The claim about robustness to incorrect graphs needs to be justified with regards to why it is required and what it intends to achieve. A plot with NL2E on y and SHD on x axis could be interesting to see what the tolerance of the algorithm is relative to the size of the problem. Consistency across both synthetic and real world datasets is a must.  \n\nWhile the second adjustment is critical to securing my recommendation for acceptance, the first one is slightly less critical but still important for clear messaging and hence acceptance.", "broader_impact_concerns": "Most concerns are already addressed.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)