{'fraction_non_alphanumeric': 0.03371659725405018, 'fraction_numerical': 0.012454416834386146, 'mean_word_length': 4.663469134409209, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Partial order: Finding Consensus among Uncertain Feature Attributions

Post-hoc feature importance is progressively being employed to explain decisions of complex machine learning models. Yet in practice, reruns of the training algorithm and/or the explainer can result in contradicting statements of feature importance, henceforth reducing trust in those techniques. A possible avenue to address this issue is to develop strategies to aggregate diverse explanations about feature importance. While the arithmetic mean, which yields a total order, has been advanced, we introduce an alternative: the consensus among multiple models, which results in partial orders. The two aggregation strategies are compared using Integrated Gradients and Shapley values on two regression datasets, and we show that a large portion of the information provided by the mean aggregation is not supported by the consensus of each individual model, raising suspicion on the trustworthiness of this practice.


INTRODUCTION
The Machine Learning (ML) framework has proven to be an essential tool in many data intensive domains such as software engineering, medecine, and cybersecurity (Esteves et al., 2020;Kaieski et al., 2020;Salih et al., 2021). However, the lack of interpretability of complex models is still an important issue that limits potential trust in ML systems. For this reason, various model-agnostic techniques such as LIME (Ribeiro et al., 2016), SHAP (Lundberg and Lee, 2017), and Integrated 1 Main correspondance: gabriel.laberge@polymtl.ca Gradient (Sundararajan et al., 2017) have recently been developed to provide explanations of model decisions in the form of feature attributions. These attributions are meant to indicate the contribution (positive/negative) of individual features toward the model prediction, and their magnitudes (positive) can be used to rank features in order of importance.
As researchers and practitioners have started to apply these model-agnostic explanations to real-world settings, it has become apparent that they are subject to variability. Indeed, recomputing the explanation and/or retraining the model can induce a high variance in feature attributions (Visani et al., 2020;Shaikhina et al., 2021), and contradictions in feature importance rankings (Fel et al., 2021;Lyu et al., 2021;Zhou et al., 2021). This is a critical issue for ML explainability: if two practitioners who use the exact same data and learning algorithm cannot agree on the feature attributions/importance, how can they trust these explanations?
Information can be extracted from diverse explanations by aggregating them using various mechanisms. A common aggregation strategy is the arithmetic mean of feature attributions, which was previously used to combine explanations obtained from multiple reruns of the explainer (Molnar, 2020, Chapter 14), explanations of points in the neighborhood of the input of interest (Bhatt et al., 2020) and explanations of independently trained models (Shaikhina et al., 2021).
The alternative combination mechanism we introduce is two-fold. One the one hand, to mitigate variability in the explainer, we compute Confidence Intervals (CIs) that capture the ground-truth attributions with high probability. On the other hand, to address variability in explanations of diverse models, we take the consensus of feature importance among independently trained models. Both strategies are used in tandem, resulting in partial orders of feature importance (instead of the more common total orders). These partial orders have the ability of abstaining from considering a feature as more important than another, making them useful at ignoring misleading information. Moreover, they are constructed in such a way that the information they do contain is meaningful, with high probability.
On two regression datasets, the mean and consensus strategies are applied to the SHAP and Integrated Gradient explanations of Multi-Layered Perceptrons (MLP) and the following observations are made.
1. The partial orders resulting from the consensus aggregation are informative (structured), even when combining explanations of about 1K independently trained MLPs.
2. The total orders induced by the mean aggregation contain information that is not supported by the consensus, while information in the latter is supported by the former. This raises doubts on the trustworthiness of the mean strategy.
The remainder of this paper is organised as follows: Section 2 presents the necessary background knowledge, while Section 3 discusses the sources of variability (uncertainty) in feature attributions as well as aggregation strategies. Section 4 lays out the experimental details, results and discussions, and finally Section 5 concludes the paper.

BACKGROUND
Feature attribution methods allow to order features based on their importance in the prediction of a single black-box model. However the ML methodology inevitably leads to a diversity of models and our approach to feature importance in this context relies on the notion of partial order.

Additive Feature Attribution
The rapid growth of the field of eXplainable Artificial Intelligence (XAI) is primarily driven by the desire to explain black box models (Arrieta et al., 2020). Formally, given an input space X ⊆ R d , an output space Y, a black box model h : X → Y, and a specific input of interest x ∈ X , there is growing interest in capturing the reasons/logic behind the decision h(x).
One family of techniques that aim at providing such information are feature attributions, which are vectorvalued functionals 1 φ : H → R d whose vector output represents the contribution of each feature towards the prediction h(x).
This paper focuses on additive feature attributions which are further required to distribute the difference between the prediction h(x) and the average prediction with respect to a chosen background distribution B over X among the features, namely feature attributions which satisfy The component φ i (h) is interpreted as the increase/decrease of the model w.r.t the background we can attribute to feature i. Additive feature attributions are motivated by a contrastive question of the form: Why is h(x) so high/low compared to the average prediction on B?
The notion of feature importance refers to the magnitude of attributions |φ i (h)|, and can be used to define total orders ≤ φ(h) commonly used to rank features (Ribeiro et al., 2016;Bhatt et al., 2020;Molnar, 2020;Jiarpakdee et al., 2020) i We now recall two additive attribution methods: Shapley values and Integrated Gradients, that are both increasingly being used to explain black boxes.

Shapley Values
The Shapley values are a fundamental concept from cooperative game theory (Shapley, 2016). Letting [d] be the set of all d features, and given a subset P ⊆ [d] of features, we define the replace function r P : Moreover, let π be a permutation of [d], π(i) be the position of the feature i in π, and π :i = {j ∈ [d] : π(j) < π(i)}. The Shapley values, as defined in the library SHAP (Lundberg and Lee, 2017), are the average marginal contribution of specifying the ith feature from the background distribution across all coalitions where Ω is the uniform distribution over all d! permutations of the features.

Integrated Gradient
The Integrated Gradient (IG) originates from a different background: cost-sharing in economics. It is also known as the Aumann-Shapley value, and has been previously used to compute saliency maps of Convolutional Neural Networks (Sundararajan et al., 2017;Erion et al., 2021). The general definition of the IG is The main idea of this approach is to average the gradient along linear paths between reference inputs sampled from the background and the input x of interest.

Underspecification and Ensembling
Many ML algorithms consist in stochastic methods for learning a model from data and therefore potentially lead to a diversity of models. For simplicity, this paper focuses only on the stochasticity introduced by optimization procedures for MLPs. This includes the random initialization of the parameters, and the random order in which mini-batches are fed to the optimizer. So for given training data S and a learning algorithm A Φ with carefully selected hyperparameters Φ, we can describe the probability distribution A Φ (S) over MLPs via its sampling procedure h ∼ A Φ (S): randomly initialise the parameters of a model, perform stochastic gradient descent on the parameters to minimize the loss on S, and retrieve h.
The fact that practitioners, even with fine-tuned hyperparameters, can end up with a distribution over diverse MLPs with equivalent test performance is indicative of the underspecification of the task. The absence of evidence to choose one model over the others can be circumvented by the introduction of an aggregate predictor h E obtained by averaging predictions where E := {h (k) } M k=1 refers to an ensemble of models. A common justification for Eq. 6 is that averaging high complexity models acts as a regularizer, meaning h E often performs better on the test set than any individual model in E (Goodfellow et al., 2016, Section 7.11). However, a recent study suggests that in the presence of a mismatch between training/testing and deployment domains, there is a considerable variance of deployment performance between independently trained models, which could reduce the benefits of aggregating predictors (D'Amour et al., 2020). In light of those results, D'Amour et al. suggest that the whole set of good predictors that a learning algorithm can return carries valuable information about underspecification of the ML pipeline.
When feature attribution techniques like SHAP and IG are used to explain the predictions of different models trained on the same data with the same learning algorithm, they can produce different results (Shaikhina et al., 2021), hence suggesting that underspecification also affects feature importance. We argue this underspecification could be induced by a mismatch between the training/testing and explanation domains. Indeed, computing these feature attributions requires evaluating the models on distributions of inputs that differ from the dataset (see Appendix C). Note that this domain mismatch has previously been observed and exploited to attack/manipulate SHAP attributions (Slack et al., 2020a). While one could simply decide to consider the feature importance of the single aggregated model h E to deal with underspecification, the considerations of D'Amour et al. encourage us to map out explanations of all individual models. Combining the explanations of each individual model will naturally lead us to a partial order of feature importance, which we now briefly introduce.

Partial Orders
Definition 2.1. A partial order on the set of features [d] is a binary relation that is reflexive, antisymmetric and transitive i.e. for all i, j, k ∈ [d], we have: The subtlety about partial orders is that it is not required that any two features i and j be comparable. Indeed, if neither i j nor j i hold, i and j are incomparable, denoted by i ⊥ j. The associated strict relation, symbolized with i ≺ j, is defined by i j and i = j. Total orders ≤ (typically used for feature importance) can simply be viewed as special cases of of partial orders, those in which every pair of features is comparable.
The total number of ordered pairs of features that are strictly comparable is called the cardinality of the order | | := |{(i, j) ∈ [d] 2 : i ≺ j}|. This quantity takes values between 0 (all features are incomparable), and d(d − 1)/2 (total order).
With partial orders, we are interested in the two distinct concepts of maximum element and maximal element, which coincide for total orders. Feature i is called maximum if j ≺ i for all other features j, while feature i is called maximal if there exists no feature considered more important, i.e. j s.t. i ≺ j. Note that a maximum element must be unique and is always maximal, but when there is no maximum there can be several maximal elements.
A common way to picture a partial order on a finite set is to draw its Hasse diagram.
Definition 2.2. Given a set of features [d] and a partial order , the associated Hasse diagram G( ) is a directed graph with vertices V := [d], and directed edges Finally, one of the fundamental properties of partial orders is that they can be intersected to make new ones.
Definition 2.3. Given a set of partial orders { k } M k=1 , the intersection of these relations k k is the binary relation which is also a partial order.
The intersection operation formalizes the concept of consensus between various orders.

ATTRIBUTION UNCERTAINTY
We now discuss two sources of variability in feature attributions: the explainer and model uncertainties, whose existence motivates the development of mechanisms to aggregate diverse explanations

Explainer Uncertainty
In the majority of practical settings, the attributions presented in Equations 4 and 5 are intractable and must be approximated via the Monte Carlo (MC) method (details in Appendix B). Therefore, in experiments, an estimator φ of φ will be computed and used to order features as in Eq. 2. Because MC is random in essence, the feature attributions φ will fluctuate and result in instability of the total order on features.
Fortunately, MC allows for the construction of confidence intervals (CIs) I i around φ i to capture the true value φ i with high probability, i.e., P φ i (h) / ∈ I i ≤ δ. When CIs are considered, ordering features by relative importance naturally leads us to a first partial order.
Definition 3.1. Given the approximated feature attribution φ i (h), and their CIs I i , the partial ordering where |I i | refers to the image of the CI under the absolute value map.
The following proposition shows that, with high probability, this partial order is unlikely to contradict the true total order ≤ φ(h) from Eq. 2.

Model Uncertainty
The second source of variability is the randomness inherent to the learning algorithm A Φ (S), which is hard to characterize because it depends on the inner workings of the optimizer and the training data. For this reason, we discuss aggregation strategies to combine explanations of various models.

Mean Aggregation
It has recently been proposed to aggregate explanations of independently trained models by taking arithmetic mean of their feature attributions (Shaikhina et al., 2021), which is equivalent to explaining the aggregated model, as long as the feature attributions are linear w.r.t the model, which is the case for SHAP/IG attributions 2 .
Definition 3.2. Given a linear feature attribution φ, the mean aggregation on E is defined as the feature attribution of the aggregated model h E . It yields the total order denoted by ≤ φ(h E ) .
This aggregation strategy faces two main challenges. First of all, there can be high variability between attributions of independently trained models (Shaikhina et al., 2021), a useful information that is completely lost when averaging.
Secondly, the order ≤ φ(h E ) lacks theoretical guarantees that, as new models are added to the ensemble, no order relations will suddenly switch.

Consensus Aggregation
Given the identified shortcomings of the mean aggregation, we introduce another method, the consensus strategy, which takes into account the feature importance of each individual model M k=1 φ(h (k) ) , leading to the following key definition. Figure 1: Partial order computed via consensus among two models and by considering CIs.
Definition 3.3. Given an ensemble E, we define the consensus aggregation as φ(E) : An example of partial order computed via Definition 3.3 is illustrated in Figure 1. One of the crucial properties of the consensus aggregation is that relations cannot flip when new models are added to the ensemble.
Proposition 3.2. If E ⊆ E , the only possible disagreements between φ(E) and φ(E ) are relations In other words, when adding new models to an ensemble, an existing relation between two features cannot be inverted, it can only cease to hold, therefore making them incomparable. Moreover if a relation between two features does not hold, adding new models to the ensemble will not change this fact.
One parameter left to specify for the consensus strategy is the level of confidence δ of the CIs. Considering an ensemble of M models, if the true attributions were known, one would have a ground-truth order ≤ φ(E) . However, in practice, one obtains approximated orders φ(E) , which are unlikely to be identical to the groundtruth. Since relations in the partial order are meant to be informative, we at least want to ensure that every edge in the approximation is present in the groundtruth.

Comparing Aggregation Strategies
Comparing the mean and consensus aggregations is not an easy task because they result in partial and total orders. This difficulty is exacerbated by the lack of an oracle feature importance to which both methods can be contrasted. In face of this challenge, we decided to build a pseudo-oracle by labeling order relations supported by both methods as trustworthy. For an ensemble E of models we consider the number of trust- | as measures of trustworthiness of the mean and consensus orders.
) have the same sign for all models (and similarily for feature j), then Although the proposition does not guarantee that R C = 1, it still suggests that R C should be large in practice.

Convergence of Partial Orders
Determining the ultimate partial order of feature importance would virtually require all possible models that can be sampled, even though a finite number of them (which is bounded by the number of features) would then already define it. The question therefore arises as to how many models are required to stabilize the partial orders. Considering our inability to determine these ultimate partial orders, we decide to characterize convergence in a simplified setting. The ultimate orders are chosen to be the partial orders computed using an extremely large ensemble E. One can then collect smaller ensembles E M ⊂ E consisting of M models randomly picked without replacement from E, and count the number of disagreements between the two orders φ(E M ) and φ(E) (cf. Proposition 3.2). Repeating this process for increasing values of M will give an idea of the number of models required to stabilize the partial order.

EXPERIMENTS
The consensus mechanism proposed in this paper is compared to the mean aggregation using Shapley values and IG on two regression datasets, and it is shown that the total order induced by the mean contains untrustworthy information.

Experimental Setup
Given a dataset, we start by splitting the instances into train/tests sets with 90% and 10% ratios respectively. Afterward, a MLP architecture is chosen: for datasets with fewer than 10K instances, a single hidden layer with 100 units was used, while for larger datasets, four hidden layers of sizes [100,50,20,10] were considered. For simplicity, the Adam optimizer with default hyperparameters was used, and a grid-search over learning rates, number of epochs, and choice of learning rate scheduler was conducted. With the chosen hyperparameters, around 1000 models were trained and a pool of models E was formed with all MLPs that had a higher test-set performance than a Random Forest baseline. This choice was made because the learning algorithm would sometimes, but rarely, yield models with terrible performances (even worse than linear models).
After obtaining the pool E, we compute the additive feature attributions of the predictions h(x) ∀h ∈ E on some observations of interest x taken from the test set. Following the Formulate-Estimate-Explain (FAE) framework introduced for SHAP/IG explanations (Merrick and Taly, 2020), we apply the following steps: Formulate Given an instance x to explain, formulate a contrastive question, which directly specifies the background distribution B for the additive attribution.
Estimate Feature attributions φ are estimated with simple MC yielding φ and CIs I i . To assess whether or not enough MC samples are used, the relative attribution error , is reported for each model h and used as a proxy of MC error. During our experiments, when the attribution error for a model was above 2%, this attribution was rejected and not considered in the aggregation.
Explain Given MC estimates and CIs (with signifi- , the mean attribution of individual features φ i (h E ) are reported in the graph nodes, and used as color in a hot/cold heatmap. For quantitative comparisons between the orders, the trustworthiness ratios R M and R C are reported.
Finally, to observe the convergence of partial orders, the process discussed in Section 3.4 was repeated 10,000 times for every value of M ∈ {1, 2, 4, 8, . . . , 256, 512}, and each time recording the number of disagreements between φ(E M ) and φ(E) .
We now present results from the two regression datasets separately, followed by a discussion of the merits and limitation of our proposed approach. Additional results are presented in Appendix E.

House Price Prediction
The Houses regression dataset, available on Kaggle 3 , consists of predicting the selling price of 1446 houses based on 79 numerical and categorical features. For simplicity, only numerical features were selected and those with high Spearman correlation (> 0.8) were removed, leading to a total of 20 features. The IG attribution was considered for this dataset. After some data exploration, we became interested in explaining the price of the costliest house from the test set.

Formulate
The instance x to explain was the test set house with the largest selling price (430,000 USD predicted as 371,953 USD by h E ), and was compared to the training set B := {Training Set}.
Approximate For MC approximation of IG, a total of 50,000 samples where used leading to relative attribution errors of at most 0.5%. Results for SHAP are reported in Appendix E.1.1.
Explain Figure 3 illustrates the feature importance of each model in the pool of size 1115, as well as for the aggregate predictor h E , while Figure 2 presents the Hasse diagram of the corresponding consensus partial order G( φ(E) ) whose nodes are colored w.r.t the mean aggregation. Given the high trustworthiness ratio of the partial order (R C = 1), we can analyse the structure of the Hasse diagram with confidence.
First of all, while 1stFlrS=very large has the largest attribution variance (as evidenced by the spread of blue lines at the top of Figure 3), our approach brings to light that there is consensus among models that the feature has maximum importance. This provides strong evidence that the surface of the first floor explains most of the high price of the house. Secondly, feature OverallQual=8 (quality of materials and finish of the house from a scale of 1 to 10) is more important than 13 other features. Although this feature is of lesser importance compared to first floor surface, this observation still suggests that material/finish quality is primordial to increase this house selling price. The third most important feature according to mean aggregation is the surface of the second floor 2ndFlrSF=0, which is given an average negative attribution of −17, 324 USD. This can appear intuitive at first because one can imagine the house's price would increase if it had a second floor.    However, our approach reveals that the relative importance of this feature varies among models and as result this feature is not considered more important than any other feature. Therefore, to stay cautious, we would abstain from stating that the absence of a second floor explains a substantial drop in the house's price.
Figure 6 (left) shows the convergence of the partial order with up to 512 models considered out of 1115 models. We observe that it takes up to 128 models to guarantee less than 30/190 (∼ 16%) disagreements between the resulting diagram and its limit. Most strikingly, considering a single model leads to around 145 order relations that would disappear when adding more models.

Bikes Rental Forecast
The Bikes (aka Bike Sharing) dataset available on the UCI repository 4 aims at predicting the hourly count of bike rentals between years 2011 and 2012 in the  Washington state, based on time and weather features. Over all available features, only the felt-temperature was removed because of its high correlation with temperature, leading to 10 features total. During data exploration, it was reported that a lot of bikes were rented at 7h-8h on WorkingDay=True. However, there exists an instance in this time period with very low amounts of bike rentals, which we are going to explain.
Formulate Explain the test set instance x with the least amount of bike rentals (44 bikes predicted as 63 by h E ) between 7-8h on workingdays in 2011. The background B is chosen as the whole training set conditioned on yr=2011, workingday=True and 7 ≤ hr ≤ 8. It was decided to fix the year to 2011 in the background in order to find attributions that are more relevant to the contrastive question.
Approximate MC was conducted by independently running the Permutation explainer of the SHAP library 5 400 times, with each call using a subset of 100 random instances from the background. As explained in Appendix B.3, this procedure is equivalent to simple MC with 400 samples. The highest relative attribution error reported was 2%. Results for IG are discussed in Appendix E.2.1.
Explain The feature importance of all models in the pool of size 988 is shown in Figure 5 and the Hasse dia- gram is presented in Figure 4. Once again, we observe that all information contained in the partial order is supported by the total one (R C = 1), while only one fifth of the information in the total order could be considered trustworthy. More specifically, we note that feature humidity=very small is found to be given the most importance by h E to explain the drop in bike rentals. By itself this observation is already quite suspicious because prior knowledge suggests that people would rather be more inclined to rent bikes when humidity is low. Looking at Figure 5, we observe an important variability in the importance that different models attribute to humidity, and looking even closer, we can see that models further disagree on the relative importance of humidity compared to other features, which is why the feature is incomparable to any other in the partial order. It appears that presenting humidity as the most important feature, as done by the arithmetic mean, may be problematic and risky. Similar observations can be made for feature mnth=March. This example shows that using the arithmetic mean can be misleading as the relation of two of the top-three most important features to other features are not supported by the consensus strategy. Figure 6 (right) illustrates the convergence of the Hasse diagram. As previously, hundreds of models are needed to obtain stability, e.g. it takes 256 models to ensure that less than 7/45 (∼ 15%) of the order relations disagree with the limit. When only considering a single model, we see that about 17 order relations would disappear as the ensemble size is increased, only leaving 4 relations in the limit.

Discussion
The main takeaway from the case studies is that combining attributions via the arithmetic mean (or equivalently explaining the aggregate model h E ) often yields untrustworthy information about feature importance that is not supported by the consensus strategy. Moreover, we find that partial orders contain information (i.e. there are features which are deemed more important by every single model), even when considering up to ∼1K MLPs, which suggests that SHAP and IG feature attributions can still drive decision-making despite being subject to variability.
As presented in Figure 6, hundreds of models are required for the partial orders to stabilize, which is quite high compared to the size of ensembles commonly used by practitioners. However, we argue that reporting the number of relations present in φ(E M ) that are missing in the limit is a pessimistic measure of diagram stability. Indeed, some differences between φ(E M ) and φ(E) may be more critical than others for endusers. For example, disagreements of order relations between features deemed irrelevant by domain knowledge may have a lesser impact on interpretability than disagreements between features relevant to the user. It is therefore necessary to study more task-specific measures of stability for partial orders.
The performance of individual models is of high concern when combining explanations via models consensus because it suffices that a single model disagrees with all others to make two features incomparable. For this reason, there is a strong assumption that all models are equally as trustworthy. In this study we used the performance of a Random Forest baseline to fix an upper bound and select acceptable models, but the exact role of performance should be investigated further in the future.

Related Work
Confidence intervals for feature attributions have recently been applied for LIME (Visani et al., 2020) and SHAP (Merrick and Taly, 2020). To the best of our knowledge, we are the first to apply said confidence intervals to define partial orders of feature importance, via Definition 3.1. Moreover, confidence intervals with δ = 5% are systematically encouraged in the literature, although applying the union bound to Proposition 3.3 suggests that δ should decrease by a factor d(d − 1). Reducing δ so drastically will increase the width of the intervals, resulting in less relations in the partial order. Nonetheless, we see this as a natural price to pay in order to ensure remaining edges convey trustworthy information.
The explanation uncertainty induced by the diversity of good models has previously been investigated (Shaikhina et al., 2021;Fisher et al., 2019). While Shaikhina et al. characterize said uncertainty by computing the variance of feature attributions among independently trained models, Fisher et al. search for the maximum and minimum feature importance among all models with acceptable performance. Our approach is fundamentally different in that we are not measuring the variability of attributions of each single feature feature i, but focus instead on whether or not every single model agrees that feature i is more important than feature j. Moreover, as presented in the case study of Section 4.2, the feature with highest variance in attribution, can still be the most informative according to the partial order, suggesting attribution variance offers an incomplete view of explanation uncertainty.
An alternative to the frequentist approach to attribution uncertainty described in Section 3.1 is offered by BayesSHAP (Slack et al., 2020b), which leverages principles from Locally Weighted Bayesian Regression to compute posterior distributions over plausible Shapley values. However, while BayesSHAP can provide confidence intervals for attributions, the method does not consider model uncertainty and focuses on explaining a single model.

CONCLUSION
In this work, sources of uncertainty for feature attributions were described. A novel aggregation mechanism based on consensus among models was introduced and compared to the arithmetic mean. It was observed that the total orders resulting from the mean attributions can contain possibly misleading information that is absent from the more conservative partial orders provided by the consensus aggregation. Moreover, the partial orders were shown to be informative (structured), even when considering ∼1K MLPs, suggesting that there exists partial order in the seeming chaos of models produced by a learning algorithm and that IG/SHAP attributions are useful tools despite suffering from variability.
As future work, we intend to study other post-hoc expla-nations (LIME, SmoothGrad, Permutation Importance, Entropic Variable Projection), employ alternative models such as Gradient Boosted Trees, and apply our methodology to more practical settings on which the nuance introduced by partial orders will hopefully prove most beneficial.

References
From the union bound, the latter event occurs with probability at most 2δ.
Since the true order ≤ φ(h) is total, for each of the d(d − 1)/2 unordered pair of features {i, j}, either i < φ(h) j or j < φ(h) i (we assume that φ i (h) = φ j (h) is impossible considering the attributions are real numbers). Hence, using the union bound, the probability that there exists some unordered pair {i, j} such that j ≺φ

A.2 Proposition 3.2
Before proving the proposition, we wish to specify what are disagreements between orders. For two distinct features i and j and a partial order , there are exactly three mutually exclusive situations about relative feature importance: either i ≺ j, j ≺ i, or i ⊥ j.
We say that two partial orders 1 and 2 disagree on an unordered pair of features {i, j} if they fail to place {i, j} in the same configuration (among the three possible ones stated above). Thereof, six type of disagreements can occur between arbitrary orders.
Proposition A.2 (Proposition 3.2). Given two ensembles E and E such that E ⊆ E , the only possible disagreements between φ(E) and φ(E ) are relations i ≺ φ(E) j for which i ⊥ φ(E ) j.
Proof. It suffices to prove that four out of six possible disagreements between orders are impossible. Now let us show that i ≺ φ(E) j and j ≺ φ(E ) i is impossible. Indeed, suppose that j ≺ φ(E ) i and so by the above remark j ≺ φ(E) i. This means that j φ(E) i and j = i. Now if i ≺ φ(E) j also holds, i.e. i φ(E) j and i = j, we get by antisymmetry that i = j, a contradiction. Similarly, j ≺ φ(E) i and i ≺ φ(E ) j is impossible.
But then i ⊥ φ(E) j, which stands for i φ(E) j and j φ(E) i does not hold. Finally, i ⊥ φ(E) j and i ≺ φ(E ) j is also impossible for similar reasons.

A.3 Proposition 3.3 Proposition A.3 (Proposition 3.3). Given an ensemble
Proof. Since i ≮ φ(E) j, there exists at least one model (lets call it h (m) ) for which |φ i (h (m) )| ≥ |φ j (h (m) )|. Now note that the event i ≺ φ(E) j is included in the event i ≺ φ(h (m) ) j whose probability is bounded by 2δ (from Proposition 3.1).

A.4 Proposition 3.4 Proposition A.4 (Proposition 3.4). Given an ensemble
j , if the estimate attributions φ i (h (k) ) have the same sign for all models (and similarily for feature j), Proof. Because φ i (h (k) ) has the same sign for all models (and similarly for feature j), the absolute value function commutes with the arithmetic mean. Therefore assuming i ≺ φ(E) j, we have and so i < φ(h E ) j as desired.

B.1 Introduction
Monte-Carlo is fundamentally about approximating expectations of random variables. Letting ψ ∈ R n be a random vector, and (ψ (1) , ψ (2) , . . . , ψ (N ) ) be a random sequence of N iid observations of ψ, we wish to approximate the expectation µ = E[f (ψ)] for some function f : R n → R given our N iid observations. The frequentist interpretation of probabilities would suggest that the following is a good estimator of µ, the quality of which is assessed using the following theorem.
We now discuss how to apply MC to estimate the Integrated Gradient and Shapley feature attributions. To improve the clarity of the analysis, we shall make the dependence of the attribution on the input x of interest and the background distribution B explicit by writting φ(h, x, B), unlike in the main paper.

B.2 Integrated Gradient
As a reminder the IG attribution of feature i is the following By defining the random vector ψ = (z, α) ∼ B × U (0, 1), we can rewrite IG in the general form whose MC estimate is immediately apparent

B.3 Shapley Values
We once again show the definition of the Shapley attributions By letting ψ = (z, π) ∼ B × Ω be a random vector, sampling it N times iid, and taking the empirical average of the expression inside the brackets of Equation 14, we recover the IME algorithm (Štrumbelj and Kononenko, 2014).
However, this approach has a significant drawback: it evaluates the model h on redundant inputs. Indeed, several samples of r π:i (z, x) can yield the exact same point in R d , because some background samples may have similar feature values as x. Considering the cost of evaluating large ensembles at any given input, it is preferable to ignore duplicate inputs.
To this end, the SHAP library provides efficient estimates of Equation 14 via the PermutationExplainer. This algorithm is optimized for speed and therefore the SHAP source code discourages feeding the whole background distribution B to the explainer, but rather use a subset of 100 iid samples from B, which we call a background minibatch.
Before introducing the algorithm in detail, it is necessary to clarify and introduce notation. As stated previously, π is a permutation over [d], π(i) represents the position of feature i in π, and π −1 (m) is the feature at position m. Moreover π :i = {j ∈ [d] : π(j) < π(i)} and π i: = {j ∈ [d] : π(i) < π(j)} are the subsets of features that appear before and after i in π.
Algorithm 1 presents the pseudo-code of the procedures within the PermutationExplainer (vectors start at index 1).
Algorithm 1 PermutationExplainer given a budget of p permutation samples. Input: h, x Parameters: Background minibatch (z (1) , . . . , z (100) ) ∼ B 100 and Permutations (π (1) , . . . , π (p) ) ∼ Ω p Output: Note that the complex optimizations made by SHAP are hidden in the repeated calls of the function Efficiently evaluate, which is able to reuse model evaluations from its previous call.
The efficiency of this estimator comes at a cost, the loss of the iid assumption necessary to the application of the Central Limit Theorem. Indeed, once again thinking of ψ = (z, π) as a random vector, we observe that Algorithm 1 does not provide iid samples of ψ, but rather a 100 × 2p grid. This is because, the algorithm considers all pairings between permutations and background minibatch samples.
We observe that Algorithm 1 boils down to computing Taking the expectation of this expression w.r.t the random vector ψ = (z (1) , z (2) , . . . , z (100) , π (1) , π (2) , . . . , π (p) ) ∼ B 100 × Ω p , we obtain where to go from the first line to the second we used the following property based on the symmetry of the distribution Ω over permutations.
Since PermutationExplainer is an unbiased estimator, running it N times with iid permutations and background minibatches, and averaging results is a MC estimate of Shapley values PermutationExplainer(h, x; ψ ( ) ).

C DISTRIBUTIONAL SHIFT IN EXPLANABILITY
This section presents a preliminary investigation of why independently trained models can yield diverse (and even contradicting) feature attributions/importance. We provide empirical evidence of a distributional shift between the train/test domains and the explanation domain.
We denote probability distributions as D and their empirical counterparts by D := 1 N N =1 δ z ( ) which attribute uniform weights to each instance of a finite set {z (1) , . . . , z (N ) } sampled iid from D. Moreover, letting f be a (measurable) function, we define the push-forward distribution f • D of the distribution D mapped thoughf i.e. (20)

C.1 Investigating Explanation Distributions
All models in the ensemble are minimizing the empirical risk on the training data and the selection of hyperparameters and models is guided by the same performance metric. While this indirectly ensures that all models will give similar predictions on the train/test instances, they can disagree on other inputs. As Monte-Carlo (MC) estimates of the SHAP and IG attributions rely on the evaluation of the models on points which may not belong to the train/test datasets, we hypothesize that our models disagree more on MC points leading to variability in feature attribution. It can be argued the underspecification of our models (materialized by a stronger disagreement on MC points) translates to the underspecification of the feature importance, which our partial order based consensus approach aims to address.
To investigate this hypothesis, we now formalize the distributions of input generated during MC estimations of IG and SHAP, which we refer to as Explanation Distributions. We then study the distributional shift they represent from the empirical data distribution D Data , and more specifically, how the shift affects the models confidence as measured by the pointwise predictive variance

C.1.1 SHAP Distribution
When running Algorithm 1 to estimate Shapley values, the models will be evaluated on inputs r π:i (z, x) sampled by replacing subsets π :i of features of x with features of a background instance z ∼ B. Unless features are statistically independent (which they rarely are) the action of replacing certain features will generate unrealistic inputs, e.g. a very humid and hot day of January, or a house that was renovated before it was build. This observation has already been made and exploited to attack SHAP explanations by modifying the model on those unrealistic inputs, resulting in arbitrary explanations without affecting the test performance (Slack et al., 2020a). However, connections with distributional shift in ML remains to be fully explored.
By assuming that only training instances as are ever used in background samples while test set points are systematically being explained, we can formalize the distribution D Shap of all points that will ever be generated when computing SHAP attributions in practice, see Algorithm 2.

C.1.2 IG Distribution
When estimating IG with MC, the models will be fed random convex combinations of the input to explain x and background samples z ∼ B. Formally, all points that will ever be generated when using IG in a practical settings, assuming only train points are ever used for background and only test points are explained, are sampled from the distribution D IG described in Algorithm 3.

C.2 Toy Example
This toy regression dataset previously used by Botari et al. (2020) generates input points along a 2D spiral and labels then with the cumulative length of the spiral up to that point. Due to the manifold structure of the dataset, it is easy to visualize which inputs are close or far from the data, an information that is supported by the variance (disagreement) among model predictions according to Figure 7. We indeed observe that the pointwise variance between models is smaller on inputs close to the manifold compared to elsewhere by orders of magnitude.
Nonetheless, since the explanation distributions differ from the data-generating one, some inputs generated by SHAP and IG fall outside the manifold, and even into regions where model disagreements are ∼ 100 times larger than on the data points. for the Spiral toy dataset. We see that the 10 models disagree more, as measured by ∆(x), on some samples from the Explainer Distributions than on most training points. .

C.3 Real Datasets
When studying higher-dimensional datasets, we cannot as easily visualise the distributional shift. An alternative is visualize histograms of the empirical distributions log • ∆ • D with D = D Train , D Test , and D Explain , see Figure   8. In total, 20 models were used to compute the pointwise variance and the empirical explanation distributions D Explain were constructed by running Algorithms 2 and 3 as many times as there are instances in the datasets.  Looking at the top right plot, it is apparent that there is a considerable portion of the IG samples generated inputs on which the models disagree more than on train/test instances. The same behavior is also present on the top left plot (SHAP), although it is less pronounced. The difference between IG and SHAP could be explained by the fact that Bikes contain many features (hour, day of the week, month, and year) structured as a regular grid. Since linear interpolations between two vertices of a grid are unlikely to land on another vertex, models need not agree on IG points even though they are trained on the same data. On the other hand, replacing some feature values of a vertex with another one (as done in SHAP) necessarily yields a vertex of the grid. Focusing on the bottom row of Figure 8, the difference between the train/test/explain distributions are not as striking, although IG has a slight tendency to generate inputs on which the models agree more than on most train/test instances.
To refine our analysis, we computed the Mann-Whitney U-statistic as a principled mean to compare the explanation distributions and the dataset distribution D Data (train and test). Formally, given two empirical distributions D 1 , and D 2 , the U statistic in common language effect size is defined by and values close to 0 or 1 when one distribution is significantly shifted compared to the other. We employed the statistic in the one-sided Mann-Whitney test of SciPy 6 , whose null hypothesis H 0 : D 1 = D 2 is rejected in favor of the alternative H 1 : D 1 < D 2 when the U statistic is high enough compared to 1/2. Table 1 presents the U-statistics (and their respective p-values) resulting from comparing the datasets and explanations distributions across 10 reruns to account for the randomness of Algorithms 2 and 3. The main observation is that three out of four null hypotheses were rejected with high significance (all p-values way smaller than 1%), meaning that in these scenarios we can clearly measure a shift of model confidence between the dataset and explanation distributions. Taking a closer look at the U-statistics, we note as previously that, on the Bikes dataset, IG tend to generate inputs on which the models have a higher disagreement compared to SHAP, as evidenced by its higher U-statistic (0.756 versus 0.631).
Interestingly, IG yields a U-statistic below 0.5 on Houses, suggesting that models tend to agree more on interpolations between train/test instances than on the dataset as a whole. Since we used a one-sided test, the associated p-value of 1 means that we failed to provide strong evidence that models trained on Houses have larger disagreements on IG points than on the whole dataset. However, despite this fact, we still ended up with a partial order when taking the consensus of IG feature importance in Section 4.2, meaning that Table 1 presents an incomplete view of the underspecification of feature attributions. In truth, while we consider here all the inputs that occur in explaining any test instance with the whole training set as background, it could happen that the Explainer Distribution corresponding to a specific instance (and possibly different background) shows a different behavior, possibly suggesting that models have diverse explanations for that specific instance. Indeed, this Table  is principally meant as a preliminary mean to illustrate the general phenomenon of distributional shift occurring when applying a given explanation method to a specific dataset, and not as a definitive proof of why models have diverse explanations.

D EXPERIMENTAL DETAILS
All experiments involving MLPs were performed using a GPU (NIVIDIA GeForce RTX 2080 Ti) resulting in around 50 hours of runtime total. To allow for the fast training and explanations of the 1K MLPs, we made sure to systematically store 50 models on the graphics card and evaluate them in parallel. Otherwise, we would not have been able to run SHAP PermutationExplainer 400 times for 1000 models in ∼ 1h.
The online platform Weights & Biases 7 was used to help design experiments and to fine-tune hyperparameters.
Notably it performed grid searches over number of epochs, learning rate, and learning rate scheduler 8 , following Table 2. Given a specific set of hyperparameters Φ, the procedure shown in Algorithm 4 was employed with m = 5 and M = 10 to evaluate their fitness. The hyperparameters with highest fitness were then used to train ∼1K MLPs on the whole training set. Results are presented in Figure 9. Note that the Hasse diagram is identical to what we had with IG, which is a surprising result considering SHAP and IG are completely different mathematical quantities. However, considering that SHAP and IG attributions coincide on additive models, we hypothesize that feature interactions are very weak in Houses compared to the main effects, an hypothesis that still needs to be empirically confirmed.

E.1.2 IG on Test Instance B
Similarly to instance A, we aim at explaining extreme values of house prices from the test set, but now look at low prices.
Formulate The instance B to explain was the test set house with the lowest selling price (79,000 USD predicted as 95,041 USD by h E ), and was compared to all training instances.
Approximate For MC approximation of IG, a total of 50,000 samples where used leading to relative attribution errors of at most 0.5%.
Explain Figure 10 shows the results on a pool of 1115 selected models. Approximate We conducted MC by independently running the PermutationExplainer 200 times, with each call using a subset of 100 random instances from the background with the default number of permutations from SHAP. The highest relative attribution error reported was 1.875%.
Explain After accounting for performance and relative attribution errors, 990 models were considered for the consensus and the results are shown in Figure 11.

E.2.3 SHAP on Test Instance C
As we were exploring the dataset, it came to our attention that instances with the most bikes rented between 10 and 17h were on off-days, which is intuitive. However, there was still a lot of variance in bike rentals within this subset of the data and we attempt to explain one of the extreme values observed.
Formulate Explain the test set instance x with the most amount of bike rentals (715 bikes predicted as 722 by the aggregate model) between 10-17h on off-days in 2012. The background B is chosen as the whole training set conditioned on yr=2012, workingday=False and 10 ≤ hr ≤ 17.
Approximate MC was done as previously resulting in a maximal relative attribution error of 0.36%.
Explain A total of 990 were selected and their aggregated explanations are presented in Figure 12.