MCCE: Monte Carlo sampling of valid and realistic counterfactual explanations for tabular data

MCCE: Monte Carlo sampling of valid and realistic counterfactual explanations for tabular data Annabelle Redelmeier, Martin Jullum, Kjersti Aas and Anders Løland Norwegian Computing Center, P.O. Box 114, Oslo, NO-0314, Norway. *Corresponding author(s). E-mail(s): [email protected]; Contributing authors: [email protected];[email protected]; [email protected]; Abstract We introduce MCCE: Monte Carlo sampling of valid and realistic Counterfactual Explanations for tabular data, a novel counterfactual explanation method that generates on-manifold, actionable and valid counterfactuals by modeling the joint distribution of the mutable features given the immutable features and the decision. Unlike other on-manifold methods that tend to rely on variational autoencoders and have strict prediction model and data requirements, MCCE handles any type of prediction model and categorical features with more than two levels. MCCE first models the joint distribution of the features and the decision with an autoregressive generative model where the conditionals are estimated using decision trees. Then, it samples a large set of observations from this model, and finally, it removes the samples that do not obey certain criteria. We compare MCCE with a range of state-of-theart on-manifold counterfactual methods using four well-known data sets and show that MCCE outperforms these methods on all common performance metrics and speed. In particular, including the decision in the modeling process improves the efficiency of the method substantially. Keywords: explainable AI, counterfactual explanations, Gower distance, conditional distribution 1 arXiv:2111.09790v2 [stat.ML] 25 Jan 2024 2Monte Carlo sampling of valid and realistic counterfactual explanations 1 Introduction It is exceedingly apparent that complex, black-box AI, machine learning, and statistical models deployed in industry need sound and efficient explanations. In this paper, we discuss counterfactual explanations, a type of local prediction explanation method. Counterfactual explanations (CEs)1explain predictions, f(x), from a fitted prediction model f(·) or some other deterministic ML/AI system. It does so by providing examples that yield a different decision than the feature vector, x, one is trying to explain2. These examples, ˜ x, are called counterfactuals, because they give an idea of what features could be changed to obtain a different decision. Generating a feature vector with a different decision is sometimes referred to as generating recourse (Karimi et al,2022; Guidotti,2022). As a simple example, imagine that a bank utilizes a black-box machine learning model to decide whether or not an individual should receive a loan. The model takes the four features age,sex,salary, and defaulted last year as inputs, and outputs the probability of defaulting on the loan ∈[0,1]. If the probability is in the desired decision interval c= [0,0.1), the bank customer is granted the loan. Now suppose that a 36-year-old female customer with a salary of $68,000 who defaulted on a loan last year, is denied a loan. To explain to the customer why her loan was denied, a set of counterfactuals could be provided. One counterfactual may show that if she had a salary of $80,000 rather than $68,000 she would have received the loan. Another may show that only increasing her salary to $72,000, but without having defaulted on her loan last year would also make her eligible for the loan. For counterfactuals to be useful, they should be of low cost, i.e., limit the number and magnitude of feature changes. Low-cost counterfactuals are easier for individuals to implement and/or case handlers to understand. Counterfactuals should also be actionable, i.e., not change features that are fixed by the individual. Otherwise, they will be impossible to implement in practice. Further, the counterfactuals should lie on the data manifold, since this ensures a realistic and plausible combination of features. Finally, counterfactuals should be valid, i.e., obtain the desired decision. Without this, the counterfactual would not be a counterfactual at all. 1.1 Related work Recently there has been an explosion of papers proposing counterfactual explainers. Surveys are given by e.g., Verma et al (2021), Stepin et al (2021) and Guidotti (2022). In the most recent survey (Guidotti,2022), the different approaches are categorized according to (i) the strategy used to generate the counterfactuals, (ii) whether the approach is model-agnostic, (iii) 1We use ‘counterfactual explanation’ or ‘CE’ to refer to the literature or explanation type and ‘counterfactual’ or ‘example’ to refer to the instance produced. 2Adecision is derived from a prediction,f(x), using a pre-defined cutoff value or interval c, characterizing the desired decision. For example, if f(x)=0.39 and c= (0.5,1], then since f(x)/∈c, we give instance xa decision of 0 and say xhas received an undesirable decision. Monte Carlo sampling of valid and realistic counterfactual explanations 3 whether the approach can be used for tabular data, text or images, (iv) whether the approach is able to handle categorical data, (v) whether the approach guarantees validity, and finally, (vi) whether the approach is able to ensure actionability. The survey shows that most counterfactual explanation approaches are tailor-made for either tabular data, text, or images. The strategy used to generate the counterfactuals is categorized into four groups: (i) optimization-based, (ii) heuristic search-based, (iii) decision-tree-based, and (iv) instance-based. The majority of the approaches belong to the first two categories. In what follows we describe and briefly review the four groups of approaches. The optimization-based approaches usually first define a loss function of the form dist1(f(˜ x), y′) + λ·dist2(x,˜ x),(1) where f(·) is a prediction model, xis the original feature vector, ˜xis the counterfactual, λis some tuning parameter, y′= 1 (f(˜x)∈c) is an indicator for the desired decision (see footnote 2on page 2), and dist1(·,·) and dist2(·,·) are two appropriate distance functions (e.g., weighted L0,L1,L2norms or the median absolute deviation). The counterfactual ˜xis then found by e.g., stochastic gradient descent (Wachter et al,2017), integer programming (Ustun et al,2019), random walks (Laugel et al,2018), genetic algorithms (Dandl et al,2020;Rasouli and Chieh Yu,2022), or through a sequence of satisfiability (SAT) problems (Karimi et al,2020). The optimization-based approaches often restrict f(·) to be differentiable, meaning that they do not work for treebased classifiers like XGBoost (Chen and Guestrin,2016) or random forest (Breiman et al,1984). Further, these methods usually cannot guarantee that the counterfactual will be actionable, and most of them do not handle categorical features with more than two levels. Finally, most of the optimization-based methods produce counterfactuals that do not lie on the data manifold (Pawelczyk et al,2021). There are, however, some exceptions. The methods CEM-VAE (Dhurandhar et al,2018), REViSE (Joshi et al,2019), CLUE (Antor´an et al, 2021), C-CHVAE (Pawelczyk et al,2020), and CRUDS (Downs et al,2020) all use variational autoencoders (VAE) to find counterfactuals that are proximal and connected to input data. The approaches that use heuristic search strategies, of which two examples are VICE (Gomez et al,2020) and MOC (Dandl et al,2020), are quite similar to the optimization-based ones, but instead of using optimization to minimize a loss function, heuristic search strategies are used. Such strategies are usually much faster than optimization algorithms. The downside is that the solutions are not necessarily optimal. Moreover, the heuristic search-based methods have many of the same weaknesses as the optimization-based approaches listed above. The decision-tree based methods approximate the behavior of the prediction model with a decision tree and then exploit the logic revealed by the tree for building counterfactual explanations. One example is FT (Tolomei et al, 4Monte Carlo sampling of valid and realistic counterfactual explanations 2017). A weakness of this approach is that the logic learned by the surrogate decision tree need not be the same as that of the original prediction model. The last category, instance-based approaches, retrieves counterfactuals by selecting the most similar examples from a (simulated) data set. The approach proposed in this paper belongs to this category. According to Guidotti (2022) the instance-based methods are the best performers with respect to all the properties benchmarked in their survey. Four instance-based methods are discussed in Guidotti (2022). Three of these methods, CBCE (Keane and Smyth, 2020), NICE (Brughmans et al,2023), and FACE (Poyiadzi et al,2020), do not treat immutable features in an appropriate way. The fourth approach, NNCE (Wexler et al,2020), should not be used for privacy reasons, because it selects counterfactuals as actual rows in the training data set. As shown above, every existing method presents one or more limitations. Given that the instance-based approaches seem to be the most productive, we believe there is a demand for a method of this kind that successfully addresses and resolves the identified shortcomings. 1.2 Our contribution In this paper, we introduce MCCE (Monte Carlo sampling of valid and realistic Counterfactual Explanations for tabular data), a novel instance-based method that handles any type of prediction model and feature cardinality. MCCE consists of three main steps. In step 1, MCCE models the data distribution and the decision using an autoregressive generative model, e.g., MADE (Germain et al,2015), with conditionals estimated using decision trees. In step 2, MCCE generates a large set of samples that respect the learned data manifold using the model from step 1. Finally, in step 3, MCCE obtains counterfactual explanations by extracting the low-cost and valid samples. MCCE is different from previously proposed CE methods in that it: 1. Uses a simultaneous model for the features and the decision to ensure on-manifold and valid counterfactuals. It models the decision through an additional fixed binary feature y′= 1 (f(x)∈c), for a desired decision interval c(see footnote 2on page 2). When generating samples in step 2, MCCE conditions on y′= 1, which guides the samples to the correct decision. This increases the efficiency of the method substantially. To the best of our knowledge, this trick has not been exploited by CE methods before. 2. Generates examples that are guaranteed to obey the immutable features (e.g., age, sex), i.e., the samples are actionable. This guarantee stems from the way the model in step 1 is specified. 3. Uses decision trees to estimate the data distribution. This means that MCCE does not require the underlying prediction model to be a gradientbased classifier like CE methods that use VAEs. 4. Handles categorical features with an arbitrary number of levels. 5. Breaks up the task into independent steps that can easily be altered to target specific needs. Monte Carlo sampling of valid and realistic counterfactual explanations 5 The rest of the paper is organized as follows. In Section 2, we describe our method, MCCE, in detail. In Section 3, we compare the counterfactual explanations from MCCE with those of various competing methods and show how MCCE surpasses all other methods when it comes to speed and performance metrics. We also investigate how sensitive MCCE is to its main tuning parameter, study the characteristics of the data generated in the second step of the method, and examine the scalability of the different steps of the MCCE procedure. Section 4examines the privacy implications of MCCE before we conclude in Section 5. 2 Generating counterfactual explanations with MCCE We now present MCCE, our algorithm for generating counterfactuals. Suppose that we have a tabular data set X∈ X with dimension ntrain ×pand binary response Y∈ Y. Suppose further that we have a prediction model f(·), a set of predictions, f(xi), and a set of observations Hthat have obtained an undesirable decision, i.e., H:= {i:f(xi)/∈c}for decision interval c, that we wish to explain. Note that given x∈ X,f(x) is deterministic, i.e., always produces the same prediction. We divide the feature space, X, into mutable features Xmand immutable features Xfsuch that X=Xm× X f. Then, for each observation h∈Hwith factual xh, our aim is to create a set of J counterfactuals Eh={eh,1,...,eh,J }, such that: •Criterion 1:eh,j is on-manifold, i.e., p(Xm=em h,j |Xf=ef h,j)> ϵ, for some ϵ > 0; •Criterion 2:eh,j is actionable, i.e., does not violate any of the immutable features; •Criterion 3:eh,j is valid, i.e., f(eh,j)∈c, for the desired decision interval c; •Criterion 4:eh,j is low cost, i.e., close to the factual, xh. Given a training data set Xand these criteria, we propose to generate a set of counterfactuals with three steps: 1. Model the distribution of the mutable features given the immutable features and the decision, y′. In our implementation, we use an autoregressive generative model based on decision trees. 2. For each observation h∈H, generate a large data set of dimension K× (p+ 1) by sampling from this fitted distribution. Denote this data set Dh for each h∈H. 3. Create Ehfor each h∈Hby removing the rows in Dhthat do not obey criteria 3 and 4. Return Eh⊆Dhas the set of counterfactual explanations for this instance. Figure 1shows an outline of the three steps, and Sections 2.1-2.3 describe them in more detail. 6Monte Carlo sampling of valid and realistic counterfactual explanations Fig. 1 MCCE’s three steps illustrated on the introductory example with two immutable and two mutable features. Step 1: Fit a decision tree for each mutable feature iteratively on the previously fit features and the decision. Step 2: For each observation/prediction to explain, trace the trees down based on the immutable feature values, decision = 1, and previous sampled values, and then randomly sample a training observation in the leaf nodes and update Dh. Step 3: Among the sampled rows, find the instance closest to the original vector over the decision boundary. 2.1 Model Step 1 (model) To model the distribution of the qmutable features, Xm= {Xm 1, . . . , Xm q}, given the uimmutable features, Xf={Xf 1, . . . , Xf u}, and decision, Y′, we use an autoregressive generative model to decompose Xm| Xf, Y ′into products of conditional probability distributions as follows p(Xm|Xf, Y ′) = p(Xm 1|Xf, Y ′) q Y i=2 p(Xm i|Xf, Y ′, Xm 1, . . . , Xm i−1).(2) Then, we estimate each of the univariate conditional distributions in (2) separately. The main reason for rewriting p(Xm|Xf, Y ′) in this way is that Monte Carlo sampling of valid and realistic counterfactual explanations 7 it is typically much easier to correctly model distributions with only a single dependent variable, rather than many. The question is then how to model each of these conditional distributions. Options include the use of basic parametric probability distributions, copulas (Sklar,1959), non-parametric methods, and deep learning approaches such as MADE (Germain et al,2015). We choose to model the q−1 conditional distributions in (2) using q−1 decision trees. For example, to model Xm 2| Xf, Y ′, Xm 1, we fit a decision tree to Xm 2∼(Xf, Y ′, Xm 1) based on the data. We use decision trees because they are fast and robust to extrapolation, scale well, handle mixed (continuous and categorical) data automatically, handle arbitrary orders of interactions between the features, and do not require extra preprocessing or scaling of the data. Inspired by Drechsler and Reiter (2011) and Reiter (2005), we choose to use CART (categorical and regression trees) (Breiman et al,1984) to model these conditional distributions. CART is a tree algorithm that builds trees by recursively making binary splits on the feature space until a given stopping criterion is fulfilled. Depending on the class of the response, either the Gini impurity index or the squared-error loss can be used to measure the quality of the split. The tree can be controlled by setting other parameters like the maximum depth or number of features to check when splitting a node. We choose CART over other non-parametric models since CART is consistent (Scornet et al,2015;Chi et al,2022) and allows for high dimensionality. In addition, CART has previously been used for synthetic data generation with good results. For example, Reiter (2005) uses CART to impute sensitive data and then compares the generated observations with the true observations. Reiter show that CART generates data that are close to the corresponding population statistics and covers the corresponding data intervals. Further, Drechsler and Reiter (2011) compare different non-parametric approaches to generating data. They show that CART is one of two methods that gives the best results when it comes to recreating population means and variances. 2.2 Generation Step 2 (generation) Suppose that in step 1 we fit q−1 decision trees: T1, . . . , Tq−1where Tiis fit to Xm i∼(Xf, Y ′, Xm 1, . . . , Xm i−1). Step 2 consists of generating K×pdimensional data sets, Dh,h∈H, with samples from the fitted trees. The full generation procedure is outlined in Algorithm 1, but in brief, for each observation, h, we loop over each mutable feature, j, and sample a set of feature values that we append to Dh. To obtain the samples for feature j, we find the end node of the fitted tree Tjbased on the current row values of Dh, and sample from this node. While the feature values in Dhonly take the values of the features in the training data, our method combines the feature values in new ways such that the rows of Dhare rarely observed in the training data. See Section 4for a discussion on privacy concerns. Note that in practice, Algorithm 1is made more efficient by only following the trace down the branches of Tjfor each 8Monte Carlo sampling of valid and realistic counterfactual explanations Algorithm 1 Outline of step 2 in MCCE Input: u > 0 number of immutable features q > 0 number of mutable features (p=u+qis the total number of features) H > 0 number of observations to explain K > 0 number of samples from each tree T1, . . . , Tq−1fitted decision trees 1: for h∈Hdo 2: Let Dhbe a K×umatrix where each row is a copy of the vector of u immutable features xf h 3: Append a K×1 column vector to Dhwith y′= 1 repeated Ktimes 4: for 1≤j≤q−1do 5: Let dbe an empty vector of length K 6: for 1≤i≤Kdo 7: Find the end node of tree Tjbased on vector Dh[i, ]. 8: Let d[i] be a single (Monte Carlo) sample from the training observations belonging to this end node. 9: end for 10: Dh←[Dh|d]; ▷Append vector das a new column of Dh 11: end for 12: Remove the y′= 1 column vector from Dh 13: return Dh 14: end for Output: data sets Dhof dimension K×p unique row of Dhand then returning multiple samples from that end node (one for each duplicate row). This saves time in the first stages of the algorithm. Observations h∈Hwith the same immutable feature values can also share the same data set, Dh. 2.3 Post-processing Step 3 (post-processing) The last step of MCCE removes the rows of Dhthat do not satisfy the criteria listed at the start of Section 2. The rows that remain become the set of counterfactuals for observation h∈H. Criteria 1 and 2 are already fulfilled since the samples come from an approximation to the data distribution conditioned on the immutable features. In addition, most samples also fulfill criterion 3 since we condition on the decision when building the trees. Nevertheless, we remove the few rows of Dhwhere f(Dh[i, ]) /∈c. To fulfill criterion 4, we calculate sparsity and the Gower distance between each Dh[i, ] and factual xh. For factual x={x1, . . . , xp}and row Monte Carlo sampling of valid and realistic counterfactual explanations 9 d={d1, . . . , dp}, sparsity is equal to the number of features changed between xand dand Gower distance = 1 p p X j=1 δG(dj, xj)∈[0,1],(3) where δG(dj, xj) = (1 Rj|dj−xj|if xjis numerical, 1 dj=xjif xjis categorical, (4) and Rjnormalizes the jth feature so that it lies between 0 and 1. Although we could have used a distance function other than the Gower distance, it is our impression that the Gower distance is the most commonly used distance function in the counterfactual explanation literature, see e.g., Wachter et al (2017), Mothilal et al (2020) and Guidotti (2022). To generate one counterfactual for each h∈H, we find the minimum sparsity across all instances of Dhand remove the instances with a sparsity larger than this value. Then we find the instance in the remaining set with the smallest Gower distance. This becomes the counterfactual for observation h∈H. In this paper, we concentrate on generating a single counterfactual per instance to be explained. However, if we want to use our method to return several (say l) counterfactuals for each h∈H, we can set upper bounds for sparsity and the Gower distance and present the best linstances in Dhwith a sparsity and Gower distance less than these bounds. See Section 5for further discussion on this subject. 3 Experiments In this section, we perform four experiments. In the first experiment, we compare MCCE to a set of previously proposed methods for generating counterfactual explanations using four well-known real data sets. The second investigates how MCCE’s performance varies with the tuning parameter K, while the third studies the characteristics of the data generated in the second step of the algorithm. Finally, the fourth experiment considers the scalability of the different steps of the MCCE procedure. See also Appendices A.1 and A.2 for additional results for non-binarized data and a non-gradient based prediction model, and for a study of the importance of conditioning on the decision. All computations in this section are run on a desktop computer with a 16-core AMD Ryzen Threadripper 1950X, 3.4 GHz CPU, and an NVIDIA GeForce GTX 1080 Ti GPU, which runs Ubuntu 20.04 with Python 3.7. MCCE runs on CPU only, while most of the competing methods can benefit from being run on the GPU. Therefore, below, the competing methods are run using a GPU (GeForce GTX 1080 Ti) whereas MCCE is run on a CPU. 16 Monte Carlo sampling of valid and realistic counterfactual explanations Table 4 Experiment 2: Average and standard deviation (in parentheses) of performance metrics for counterfactuals generated with MCCE as we vary K. Data set: Adult, ntest = 1000, Kvaries K L0↓L1↓violation↓NCE ↑model(s)↓gener(s)↓post-proc(s)↓total(s)↓ 10 3.98 (1.1) 1.59 (0.94) 0 (0.0) 999 3.8 0.38 5.56 9.74 50 3.39 (0.94) 1.1 (0.77) 0 (0.0) 1000 3.8 0.58 5.61 9.99 100 3.18 (0.9) 0.93 (0.69) 0 (0.0) 1000 3.8 0.76 5.72 10.28 1000 2.7 (0.73) 0.56 (0.45) 0 (0.0) 1000 3.8 4.39 7.23 15.43 5000 2.51 (0.66) 0.45 (0.39) 0 (0.0) 1000 3.8 21.29 12.34 37.43 10,000 2.43 (0.64) 0.43 (0.37) 0 (0.0) 1000 3.8 42.61 17.71 64.11 25,000 2.39 (0.64) 0.41 (0.36) 0 (0.0) 1000 3.8 106.93 29.50 140.23 Monte Carlo sampling of valid and realistic counterfactual explanations 17 Figure 2shows histograms for four of the variables in this data set (white) with histograms of the real data superimposed (dark grey). Where the histograms overlap, the blend of white and dark grey gives a light grey color. As can be seen from the figure, the generated marginal distributions for both the continuous and categorical variables are very close to the data distributions. Appendix A.5 shows that the figures for the rest of the variables are very similar. In Figures 3-5we show the correlation matrices for the real and generated data, as well as the difference between the correlation matrices. It is evident that even the correlations are very well captured by the generative model. Another way of checking whether the generated data can easily be told apart from the original data is to compute the discriminator measure described by Borisov et al (2023). To do so, we fit a Random Forest model on a mix of generated training samples and samples from the real training data set, where the response indicates whether the observation comes from the generated or real data (when performing this experiment, the real response is discarded). Then, we use this model to predict the probability of belonging to the real data on a set of new samples and compute the model’s out-of-sample accuracy using a threshold of 0.5. If the discriminator cannot differentiate between the generated and real samples, the model’s accuracy will be around 50%. In our case, the mean accuracy over five trials with different random seeds for the Adult, GMC, FICO, and German Gredit data sets are 62.8%, 56.5%, 59.1%, and 55.7%, respectively. Although these numbers are slightly higher than 50%, the numbers for Adult and FICO are lower than the lowest accuracy reported for the same datasets in Borisov et al (2023)3. These experiments show that it is difficult to differentiate between counterfactuals generated by MCCE and the underlying data set. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Correlations of real data Fig. 3 Correlation matrix for the FICO data set. The areas of the circles are proportional to the absolute value of the corresponding correlation coefficients. 3In Borisov et al (2023) the FICO dataset is referred to as the HELOC dataset. 18 Monte Carlo sampling of valid and realistic counterfactual explanations -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Correlations of generated data Fig. 4 Correlation matrix for the data generated from the autoregressive generative model (K= 10,000). The areas of the circles are proportional to the absolute value of the corresponding correlation coefficients. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Difference between correlations for real and generated data Fig. 5 Difference between the correlation matrices in Figures 3and 4. The areas of the circles are proportional to the absolute value of the corresponding correlation coefficients. 3.4 Experiment 4: Is MCCE scalable? For explanability methods to be applicable in a wide range of practical situations, it is crucial that they scale well computationally. We now study the scalability of the different steps of the MCCE procedure. The first subsection provides actual measurements of the computation time for the three MCCE Monte Carlo sampling of valid and realistic counterfactual explanations 19 steps as we vary the different quantities, using a controlled simulation setup. The second subsection provides the corresponding theoretical computational complexity in big-Onotation. 3.4.1 Computational complexity in practice Let the p-dimensional feature vector Xbe simulated from a zero-mean Gaussian distribution with covariance Cov(Xj, Xk) = 0.5 and Var(Xj) = 1 for all j=k. In this simulated setup, we explain simple linear models of the form f(x) = Pp j=1 βjxj, where βj= 1 for all j. We generate counterfactual explanations for predictions of this model as we vary the following quantities: p,ntrain,ntest, and K. For each simulated model, we set the decision interval for an acceptable decision to c= (0,∞), which on average covers half of the observations regardless of the dimension p. We always only explain test observations that are originally outside of this interval. Moreover, we do not fix any of the features in the explanations, i.e., all features are mutable (q=p). As an initial broad overview, Table 5shows the computation times for the three steps of MCCE for all combinations of p= (5,30), ntrain = (1000,10,000), ntest = (1,50), and K= (10,000,100,000). The displayed computation times represent the mean of 10 repeated computations. Table 5 Overview of the computational times (in seconds) of the three steps of MCCE when varying ntest,ntrain ,p, and K. The displayed computation times represent the mean of 10 repeated computations. ntest ntrain p K model(s)↓gener(s)↓post-proc(s)↓total(s)↓ 1 1000 5 10,000 0.04 0.11 0.09 0.24 1 1000 5 100,000 0.04 0.75 0.12 0.91 1 1000 30 10,000 1.25 0.83 0.09 2.17 1 1000 30 100,000 1.25 5.10 0.15 6.50 1 10,000 5 10,000 0.48 0.07 0.09 0.64 1 10,000 5 100,000 0.48 0.46 0.12 1.06 1 10,000 30 10,000 12.05 0.95 0.10 13.10 1 10,000 30 100,000 12.05 5.39 0.15 17.59 50 1000 5 10,000 0.05 3.70 0.99 4.74 50 1000 5 100,000 0.04 38.26 3.36 41.66 50 1000 30 10,000 1.24 24.89 1.03 27.16 50 1000 30 100,000 1.24 264.63 6.21 272.08 50 10,000 5 10,000 0.48 2.07 0.70 3.25 50 10,000 5 100,000 0.49 25.16 3.01 28.66 50 10,000 30 10,000 12.06 24.87 1.04 37.97 50 10,000 30 100,000 12.00 268.76 6.03 286.79 Altering all quantities affects the total computation times. Even for as many as p= 30 features, performing the modeling with ntrain = 10,000 only takes 12 seconds. The dimension pdoes, however, play a crucial role for both 20 Monte Carlo sampling of valid and realistic counterfactual explanations the modeling time and the generation time. Further, the number of samples generated per test observation, K, and the number of test observations ntest, naturally affects the generation time significantly. The post-processing time is quite limited even for ntest = 50, p = 30 and K= 100,000. That is, for the model and quantity sizes used here, the generation time is the driving factor. In Figures 6-9we take a more detailed look at how each of the four quantities affects the runtime of the different steps when altering one at a time and keeping the others fixed. When the quantities are not altered, we fix them at the following values: ntest = 1, ntrain = 1000, p = 10, K = 10,000. modeling 5 20 40 60 80 100 0 3 6 9 dimension p Runtime (seconds) generation 5 20 40 60 80 100 0 1 2 3 dimension p Runtime (seconds) post−processing 5 20 40 60 80 100 0.0 0.1 0.2 0.3 dimension p Runtime (seconds) total 5 20 40 60 80 100 0 5 10 15 dimension p Runtime (seconds) Runtime as function of the data dimension, p Fig. 6 Runtime as function of dimension p. The other quantities are fixed to ntest = 1, ntrain = 1000, K = 10,000. From Figure 6, we see that increasing the dimension p, increases the modeling time faster-than-linearly, the generation time linearly, and the postprocessing time not at all. The higher rate of the modeling time indicates that for very high dimensions (>100), the modeling time may be the driving factor. Figure 7shows that the generation time scales linearly at a significant rate in terms of K. The picture is not as clear for the post-processing time, but this is less important, as the actual run times are rather small for this component. The modeling time is obviously unaffected by Kand is therefore not displayed. Figure 8shows that the generation time scales approximately linearly at quite a high rate. The effect of ntrain on the generation time is unclear from these simulations. For the smallest values of ntrain, this component is probably Monte Carlo sampling of valid and realistic counterfactual explanations 21 generation 5e+042e+05 4e+05 6e+05 8e+05 1e+06 4 8 12 16 K Runtime (seconds) post−processing 5e+042e+05 4e+05 6e+05 8e+05 1e+06 0.1 0.2 0.3 0.4 0.5 0.6 K Runtime (seconds) total Large K 5e+042e+05 4e+05 6e+05 8e+05 1e+06 5 10 15 K Runtime (seconds) 0.15 0.20 0.25 500 2000 4000 6000 8000 10000 K Runtime (seconds) 0.025 0.050 0.075 500 2000 4000 6000 8000 10000 K Runtime (seconds) Small K 500 2000 4000 6000 8000 10000 0.30 0.35 0.40 0.45 0.50 K Runtime (seconds) Runtime as function of the K Fig. 7 Runtime as function of the number of generated samples K. The other quantities are fixed to ntest = 1, p = 10, ntrain = 1000. The upper panels show K > 104, and lower panels show K≤104. The modeling time is obviously unaffected by Kand is therefore not displayed. influenced by inaccuracies in the computational timing. The post-processing time is not affected by ntrain and is therefore not displayed. Figure 9shows that both the generation and post-processing times naturally scale approximately linearly in terms of ntest. The modeling time is clearly not affected by ntest and is therefore not displayed. For this specific setup, the total computation time seems to increase by about 0.2 seconds per extra test observation. To showcase that these scalability results are relevant and valid in a broader context, we also run simulations with quantities mimicking the four real data sets from Section 3.1. Table A4 in Appendix A.6 displays the computation times for these simulations and shows a relatively good match between the computation times in the simulations and the real data experiments. This suggests that the scalability observed in the above simulation study is likely to extend roughly to real-world data scenarios. 3.4.2 Theoretical computational complexity The simulation study above shows that MCCE scales well in all four quantities. In this section, we discuss the theoretical computational complexity in terms of big-O notation. Recall that uand qare the number of immutable and mutable features, respectively. 22 Monte Carlo sampling of valid and realistic counterfactual explanations modeling 5e+04 2e+05 3e+05 4e+05 5e+05 0 50 100 150 n_train Runtime (seconds) generation 5e+04 2e+05 3e+05 4e+05 5e+05 0.2 0.3 0.4 0.5 0.6 n_train Runtime (seconds) total Large n_train 5e+04 2e+05 3e+05 4e+05 5e+05 0 50 100 150 n_train Runtime (seconds) 0.0 0.5 1.0 1.5 500 2000 4000 6000 8000 10000 n_train Runtime (seconds) 0.2 0.3 0.4 500 2000 4000 6000 8000 10000 n_train Runtime (seconds) Small n_train 500 2000 4000 6000 8000 10000 0.5 1.0 1.5 2.0 n_train Runtime (seconds) Runtime as function of n_train Fig. 8 Runtime as function of number of training samples ntrain. The other quantities are fixed to ntest = 1, p = 10, K = 10,000. The upper panels shows ntrain >104and the lower panels shows ntrain ≤104. The post-processing time is not affected by ntrain and is therefore not displayed. According to Hastie et al (2009, Ch 9.7), the training time for the i-th tree (with u+ifeatures) is O((u+i)ntrain log(ntrain)). Since we have q−1 such trees and Pq−1 i=1 (u+i) = (q−1)(u+ (q−2)/2), the computational complexity of step 1 is O(q(u+q)ntrain log(ntrain)). Generating a sample from one of the conditional distributions in step 1 requires traversing the decision tree from root to leaf, which requires going through roughly O(log(ntrain)) nodes (G´eron,2019). Hence, the computational complexity of step 2 is O(ntestKq log(ntrain)). In step 3, we compute validity, sparsity, and the Gower distance for every sample to find the best counterfactual. The computational complexity of the validity computation is non-trivial since it depends on the choice of prediction model. Assuming that the computation time of the prediction model increases no faster than linearly in u+q, the computational complexity of step 3 is O(ntest K(u+q)). Removing smaller order terms, the total computational complexity of MCCE becomes O(q(u+q)ntrain log(ntrain) + ntestK(qlog(ntrain) + u)), showing that MCCE scales very well in all quantities (i.e., quadratic in q, at rate ntrain log(ntrain) in ntrain, and linear in u, K, and ntest). Overall, these theoretical scalability terms are in accordance with our findings from the simulations: For the modeling time, the quadratic increase in q and log-linear increase in ntrain match the findings from our simulations. The Monte Carlo sampling of valid and realistic counterfactual explanations 23 generation 500 2000 4000 6000 8000 10000 0 500 1000 1500 n_test Runtime (seconds) post−processing 500 2000 4000 6000 8000 10000 50 100 150 n_test Runtime (seconds) total Large n_test 500 2000 4000 6000 8000 10000 0 500 1000 1500 n_test Runtime (seconds) 0 5 10 15 1 10 20 40 60 80 100 n_test Runtime (seconds) 0.0 0.5 1.0 1.5 1 10 20 40 60 80 100 n_test Runtime (seconds) Small n_test 1 10 20 40 60 80 100 0 5 10 15 n_test Runtime (seconds) Runtime as function of n_test Fig. 9 Runtime as function of number of test observations ntest. The other quantities are fixed to ntrain = 1000, p = 10, K = 10,000. The upper panels shows ntest >104, and upper panels shows ntest ≤104. The modeling time is clearly not affected by ntest and is therefore not displayed. linear increase in ntest, K, and qfor the generation time is also in accordance with our simulation results. The scalability in terms of ntrain was, as mentioned, fairly uncertain, partly due to the low computation time (<1 sec). Thus, it is not unreasonable that the theoretically found logarithmic rate is valid also for our implementation (with quite a small factor). For the postprocessing time, the computation times in our simulations are so small that the theoretically found linearity in q,Kand ntest seems reasonable. 4 Privacy Naive counterfactual explanations can leak sensitive or personal data in a variety of ways. First, if a counterfactual explanation is an exact copy of an individual in the training data, this individual may be recovered, leading to a breach of privacy. Table 6presents both the number and percentage of MCCE counterfactuals with exact copies in each of the datasets discussed in this paper. Since there are at most 1.1% of counterfactuals with copies, this type of privacy breach is rare for MCCE. However, in certain circumstances, it may be necessary for a strict assurance that the counterfactuals generated by MCCE do not match training observations. To achieve this, one can easily add a step in MCCE’s postprocessing step to remove all rows of Dhthat exactly match an observation 24 Monte Carlo sampling of valid and realistic counterfactual explanations Table 6 The number and proportion of counterfactuals being copies of original training observations for each dataset. Adult GMC FICO German Credit # counterfactuals 1000 1000 1000 200 # copies 11 4 0 2 Proportion of copies 1.1% 0.4% 0% 1.0% in the training data. Although this approach may slightly diminish the quality of the counterfactuals based on other performance measures, it provides a safeguard against privacy breaches of this type. Another privacy breach may occur due to so-called explanation linkage attacks (Goethals et al,2022) or membership inference attacks. A linkage attack is an attempt to re-identify individuals in a dataset by combining the generated data with background information. A membership inference attack attempts to re-identify individuals by exploiting the fact that their counterfactual explanation may reveal whether they were used in training the model (Pawelczyk et al,2023). Even though MCCE’s explanations are rarely observed in the underlying data, it is an open question as to whether MCCE counterfactuals are vulnerable to explanation linkage or membership inference attacks. One way to alleviate these potential privacy problems is to incorporate differential privacy (Dwork,2006). In the case of MCCE, incorporating differential privacy could entail introducing noise in its data generating step (see e.g., Mahiou et al (2022) or Nowok et al (2016)). Further work is required to investigate the efficiency and practical considerations of such an approach. 5 Conclusion To the best of our knowledge, MCCE is the first counterfactual method that models both the underlying data distribution and the decision. As shown in our experiments, this is a powerful setup because even without the post-processing step, MCCE generates samples that are almost always valid and actionable. In addition, modeling the data with an autoregressive generating model produces counterfactuals with lower costs than all the competing methods. Finally, as opposed to other on-manifold methods that typically use variational autoencoders and strict prediction model and data requirements, MCCE handles any type of prediction model and categorical features with more than two levels. 5.1 Future steps MCCE’s three steps are flexible and can be improved independently of each other. We will discuss possible changes to each step below. In MCCE’s first step, the decision trees can easily be replaced with another generation method like conditional generative adversarial nets (CGAN) (Mirza and Osindero, 2014) or conditional tabular GANs (CTGAN) (Xu et al,2019). Note, however, that both CGANs and GANs are restricted to gradient-based classifiers. Monte Carlo sampling of valid and realistic counterfactual explanations 25 If formal privacy guarantees are required, differential privacy can be built into MCCE’s data-generating process with little effort. In MCCE’s second step (the autoregressive generative model), the visit sequence has to be chosen. In the experiments conducted in this paper, the visit sequence corresponded to the order of variables in the original data sets. One might get even better correspondence between the generated and original data if the visit sequence for example is chosen based on the correlation between the variables, modeling the most correlated variables first. This would be an interesting and relevant topic for further work. Finally, MCCE’s third step can easily be adjusted to handle other performance metrics. Since the Gower distance penalizes changes to continuous features less than categorical features, one might choose to replace it with another type of distance function like the Heterogeneous or Interpolated Value Difference Metrics (Wilson and Martinez,1997). If, for some reason, the number of features being changed is irrelevant for a specific application, the initial filtering with the L0norm may also be dropped. We have developed MCCE for tabular data. One avenue for further work is to adapt MCCE to non-tabular data. As long as the prediction model is used for a binary (or multiclass) decision, a prerequisite for counterfactual explanations, it should be possible to adapt MCCE’s three steps to non-tabular data, such as time series, text, or images. One challenge might be that the number of samples needed to generate explanations with low enough cost becomes higher with more complex data. In this paper, the aim was to generate one counterfactual for each individual. However, as mentioned, MCCE can easily be extended to return the k∈Kbest counterfactuals if this is desired. To enhance the practical relevance of a set of multiple counterfactuals, it might be desirable for the generated counterfactuals to also exhibit diversity (Mothilal et al,2020). Finding the satisfactory balance between performance scores for single counterfactuals and their diversity remains a subject for future research. Acknowledgments. This work was supported by the Norwegian Research Council grant 237718 (BigInsight). This paper is supported by the European Union’s HORIZON Research and Innovation Programme under grant agreement No 101120657, project ENFIELD (European Lighthouse to Manifest Trustworthy and Green AI). Authors’ contributions. AR: Methodology, Software, Validation, Data Creation, Original Draft, Writing, Review & Editing, Visualization. MJ: Conceptualization, Methodology, Software, Validation, Original Draft, Writing, Review & Editing, Visualization. KA: Conceptualization, Methodology, Validation, Original Draft, Writing, Review & Editing, Visualization. AL: Conceptualization, Methodology, Writing, Original Draft, Review & Editing, Visualization. Availability of data and materials. The Adult, FICO and German Credit Data Sets can be downloaded from https://github.com/riccotti/Scamander/ 32 Monte Carlo sampling of valid and realistic counterfactual explanations MSinceMostRecentDelq Density 0 20 40 60 80 0.00 0.04 0.08 MaxDelq2PublicRecLast12M Density 02468 0.0 0.2 0.4 0.6 0.8 MaxDelqEver Density 2345678 0.0 0.2 0.4 0.6 0.8 NumTotalTrades Density 0 20 40 60 80 100 0.000 0.010 0.020 0.030 Fig. A3 Histograms for four of the variables in the generated data set (white) with the histograms for the real data superimposed (dark grey). Where the histograms overlap, the blend of white and dark grey gives a light grey color. NumTradesOpeninLast12M Density 0 5 10 15 0.0 0.1 0.2 0.3 0.4 0.5 PercentInstallTrades Density 0 20 40 60 80 100 0.000 0.010 0.020 MSinceMostRecentInqexcl7days Density −5 0 5 10 15 20 25 0.00 0.10 0.20 NumInqLast6M Density 0 10 20 30 40 50 60 70 0.00 0.05 0.10 0.15 Fig. A4 Histograms for four of the variables in the generated data set (white) with the histograms for the real data superimposed (dark grey). Where the histograms overlap, the blend of white and dark grey gives a light grey color. Monte Carlo sampling of valid and realistic counterfactual explanations 33 NumInqLast6Mexcl7days Density 0 10 20 30 40 50 60 70 0.00 0.05 0.10 0.15 0.20 NetFractionRevolvingBurden Density 0 50 100 150 200 250 0.000 0.005 0.010 0.015 NetFractionInstallBurden Density 0 100 200 300 400 500 0.000 0.004 0.008 NumRevolvingTradesWBalance Density 0 10 20 30 0.00 0.05 0.10 0.15 Fig. A5 Histograms for four of the variables in the generated data set (white) with the histograms for the real data superimposed (dark grey). Where the histograms overlap, the blend of white and dark grey gives a light grey color. NumInstallTradesWBalance Density −5 0 5 10 15 20 25 0.00 0.10 0.20 NumBank2NatlTradesWHighUtilization Density −5 0 5 10 15 0.00 0.10 0.20 PercentTradesWBalance Density 0 20 40 60 80 100 0.000 0.010 Fig. A6 Histograms for three of the variables in the generated data set (white) with the histograms for the real data superimposed (dark grey). Where the histograms overlap, the blend of of white and dark grey gives a light grey color. 34 Monte Carlo sampling of valid and realistic counterfactual explanations References Antor´an J, Bhatt U, Adel T, et al (2021) Getting a clue: A method for explaining uncertainty estimates. In: International Conference on Learning Representations Borisov V, Seffler K, Leemann T, et al (2023) Language models are realistic tabular data generators. In: Proceeedings of ICLR 2023 Breiman L, Friedman J, Olshen R, et al (1984) Classification and regression trees. Chapman and Hall Brughmans D, Leyman P, Martens D (2023) NICE: An algorithm for nearest instance counterfactual explanations. Data Mining and Knowledge Discovery pp 1–39 Chen T, Guestrin C (2016) XGBoost: A scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp 785–794 Chi CM, Vossler P, Fan Y, et al (2022) Asymptotic properties of highdimensional random forests. The Annals of Statistics 50(6):3415–3438 Dandl S, Molnar C, Binder M, et al (2020) Multi-objective counterfactual explanations. In: International Conference on Parallel Problem Solving from Nature, Springer, pp 448–469 Dhurandhar A, Chen PY, Luss R, et al (2018) Explanations based on the missing: towards contrastive explanations with pertinent negatives. In: Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp 590–601 Downs M, Chu JL, Yacoby Y, et al (2020) Cruds: Counterfactual recourse using disentangled subspaces. In: ICML Workshop on Human Interpretability in Machine Learning Drechsler J, Reiter JP (2011) An empirical evaluation of easily implemented, nonparametric methods for generating synthetic datasets. Computational Statistics & Data Analysis 55(12):3232–3243 Dwork C (2006) Differential privacy. In: Automata, Languages and Programming: 33rd International Colloquium, ICALP 2006, Venice, Italy, July 10-14, 2006, Proceedings, Part II 33, Springer, pp 1–12 Germain M, Gregor K, Murray I, et al (2015) Made: Masked autoencoder for distribution estimation. In: International Conference on Machine Learning, PMLR, pp 881–889 Monte Carlo sampling of valid and realistic counterfactual explanations 35 Goethals S, S¨orensen K, Martens D (2022) The privacy issue of counterfactual explanations: explanation linkage attacks. arXiv preprint arXiv:221012051 Gomez O, Holter S, Yuan J, et al (2020) Vice: Visual counterfactual explanations for machine learning models. In: Proceedings of the 25th International Conference on Intelligent User Interfaces. Association for Computing Machinery, New York, NY, USA, IUI ’20, pp 531–535 Guidotti R (2022) Counterfactual explanations and how to find them: literature review and benchmarking. Data Mining and Knowledge Discovery pp 1–55 G´eron A (2019) Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow, 2nd edition. O’Reilly Media, Inc Hastie T, Tibshirani R, Friedman JH, et al (2009) The Elements of Statistical Learning: Data Mining, Inference, and Prediction, vol 2. Springer Joshi S, Koyejo O, Vijitbenjaronk W, et al (2019) Towards realistic individual recourse and actionable explanations in black-box decision making systems. Safe Machine Learning workshop at ICLR Karimi AH, Barthe G, Balle B, et al (2020) Model-agnostic counterfactual explanations for consequential decisions. In: International Conference on Artificial Intelligence and Statistics, PMLR, pp 895–905 Karimi AH, Barthe G, Sch¨olkopf B, et al (2022) A survey of algorithmic recourse: contrastive explanations and consequential recommendations. ACM Computing Surveys 55(5):1–29 Keane MT, Smyth B (2020) Good counterfactuals and where to find them: A case-based technique for generating counterfactuals for explainable AI (XAI). In: Case-Based Reasoning Research and Development: 28th International Conference, ICCBR 2020, Salamanca, Spain, June 8–12, 2020, Proceedings 28, Springer, pp 163–178 Laugel T, Lesot MJ, Marsala C, et al (2018) Comparison-based inverse classification for interpretability in machine learning. In: International Conference on Information Processing and Management of Uncertainty in Knowledge-Based Systems, Springer, pp 100–111 Mahiou S, Xu K, Ganev G (2022) dpart: Differentially Private Autoregressive Tabular, a General Framework for Synthetic Data Generation. arXiv preprint arXiv:220705810 Mirza M, Osindero S (2014) Conditional generative adversarial nets. arXiv preprint arXiv:14111784 36 Monte Carlo sampling of valid and realistic counterfactual explanations Mothilal RK, Sharma A, Tan C (2020) Explaining machine learning classifiers through diverse counterfactual explanations. In: Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency, pp 607–617 Nowok B, Raab GM, Dibben C, et al (2016) synthpop: Bespoke creation of synthetic data in R. Journal of Statistical Software 74(11):1–26 Pawelczyk M, Broelemann K, Kasneci G (2020) Learning model-agnostic counterfactual explanations for tabular data. In: Proceedings of The Web Conference 2020, pp 3126–3132 Pawelczyk M, Bielawski S, Van den Heuvel J, et al (2021) Carla: A python library to benchmark algorithmic recourse and counterfactual explanation algorithms. arXiv preprint arXiv:210800783 Pawelczyk M, Lakkaraju H, Neel S (2023) On the privacy risks of algorithmic recourse. In: International Conference on Artificial Intelligence and Statistics, PMLR, pp 9680–9696 Poyiadzi R, Sokol K, Santos-Rodriguez R, et al (2020) Face: Feasible and actionable counterfactual explanations. In: Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society, pp 344–350 Rasouli P, Chieh Yu I (2022) CARE: Coherent actionable recourse based on sound counterfactual explanations. International Journal of Data Science and Analytics pp 1–26 Reiter JP (2005) Using CART to generate partially synthetic public use microdata. Journal of Official Statistics 21(3):441 Scornet E, Biau G, Vert JP (2015) Consistency of random forests. The Annals of Statistics 43(4):1716–1741 Sklar M (1959) Fonctions de repartition an dimensions et leurs marges. Publ inst statist univ Paris 8:229–231 Stepin I, Alonso JM, Catala A, et al (2021) A survey of contrastive and counterfactual explanation generation methods for explainable artificial intelligence. IEEE Access 9:11,974–12,001 Tolomei G, Silvestri F, Haines A, et al (2017) Interpretable predictions of treebased ensembles via actionable feature tweaking. In: Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. Association for Computing Machinery, New York, NY, USA, KDD ’17, pp 465–474 Monte Carlo sampling of valid and realistic counterfactual explanations 37 Ustun B, Spangher A, Liu Y (2019) Actionable recourse in linear classification. In: Proceedings of the Conference on Fairness, Accountability, and Transparency, pp 10–19 Verma S, Dickerson JP, Hines K (2021) Counterfactual explanations for machine learning: Challenges revisited. CoRR abs/2106.07756 Wachter S, Mittelstadt B, Russell C (2017) Counterfactual explanations without opening the black box: Automated decisions and the GDPR. Harv JL & Tech 31:841 Wexler J, Pushkarna M, Bolukbasi T, et al (2020) The what-if tool: Interactive probing of machine learning models. IEEE Transactions on Visualization and Computer Graphics 26(1):56–65. https://doi.org/10.1109/TVCG.2019. 2934619 Wilson DR, Martinez TR (1997) Improved heterogeneous distance functions. Journal of Artificial Intelligence Research 6:1–34 Xu L, Skoularidou M, Cuesta-Infante A, et al (2019) Modeling tabular data using conditional GAN. Advances in Neural Information Processing Systems 32