Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
55,274 characters · 14 sections · 38 citation commands
Multiply-Robust Causal Change Attribution
\usetikzlibrary{arrows.meta} \twocolumn[ \icmltitle{Multiply-Robust Causal Change Attribution}
\icmlaffiliation{aaa}{Amazon} \icmlaffiliation{mmm}{MIT Department of Economics.} \icmlaffiliation{xxx}{This work was completed while VQM was an intern at Amazon.}
\icmlcorrespondingauthor{Victor Quintas-Martinez}{[email removed]}
\icmlkeywords{Machine Learning, ICML} \vskip 0.2in ]
\printAffiliationsAndNotice
Analysts are often interested in identifying and quantifying the contribution of multiple possible causes of change to the performance metrics of large-scale systems, such as sales volume, throughput, or retention rate. As an example, a manufacturer compares data from 2023 and 2022 and realizes that net sales have increased. However, many factors have changed between these two time periods, including the characteristics of the product, competitors' prices, or market conditions. How much did each of these variables contribute to the increase in average sales? This is a change attribution question. When policymakers and business leaders are interested in using these insights for future changes, the change attribution problem necessarily becomes causal.
In the causal change attribution problem, we observe two samples of data, including an outcome and multiple explanatory variables. We are also given a causal Directed Acyclic Graph (DAG) that encodes the prior expert knowledge about the causal conditional (in-)dependence relationships among the variables in the data. The objective is to assign a score to each of the causal mechanisms in this DAG that quantifies its contribution to the change in the distribution of the outcome. Although in many instances we observe a change over time, we want to emphasize that our methods apply to more general settings involving a comparison of two samples --- for example, differences between groups (e.g., females and males) or geographical locations (e.g., East vs. West Coast of the US).
The contribution of each causal mechanism is generally very difficult to disentangle. We need to characterize what the distribution of the outcome variable would have been if we shifted only some causal mechanisms, but left others unchanged. This is a counterfactual distribution, which is not directly observable in the data. We overcome this fundamental challenge by employing a combination of regression and re-weighting methods. Moreover, our estimating equation is multipy robust, in the sense that it still recovers the parameters of interest even if some components of the model are misspecified. We show that our estimator is consistent and asymptotically normal under weak conditions on the ML algorithms used to learn the regression and the weights. The asymptotic variance is also consistently estimable, allowing us to compute valid standard errors, confidence intervals and $p$-values. To obtain these results we build on the existing double/debiased ML literature chernozhukov2018double, chernozhukov2021automatic, chernozhukov2023automatic.
There is a large body of work on the attribution problem efron2020prediction, yamamoto2012understanding, dalessandro2012causally, some of which has studied causal attribution liu2023need, lu2023evaluating, mougan2023explanation, zhao2023conditional, berman2018beyond, ji2016probabilistic, dawid2014fitting, shao2011data. An important branch of the literature has focused on developing suitable definitions of the contribution of each variable, e.g., with variations of Shapley values janzing2024quantifying, chen2023algorithms, budhathoki2022causal, jung2022measuring, sharma2022counterfactual or optimal transport methods kulinski2023towards, while remaining agnostic about estimation. Our contribution is related but distinct: the estimator we develop can be readily integrated into existing frameworks for causal change attribution, such as Shapley values, or it may be of interest in its own right. Shapley values based on our estimator will inherit its consistency and asymptotic normality properties; we also propose a convenient and computationally efficient bootstrap procedure to compute their standard errors.
Although other recent algorithms have been proposed to estimate Shapley Values, a majority of these algorithms are not multiply-robust, and generally they are based on regression only (learning the distribution of the outcome given the explanatory variables).\footnote{The only exception that we are aware of is jung2022measuring. However, their notion of $do$-Shapley values decomposes the effect of fixing a value for the explanatory variables, rather than a distribution for each causal mechanism, as we do in this paper. They consider only discrete-valued explanatory variables, and their approach does not seem easily generalizable beyond that.} As such, they will not perform well unless we have a high-quality estimator of the regression function. In high-dimensional or non-parametric settings, estimators for those objects will typically exhibit slow convergence rates, rendering traditional normal-based large sample inference (standard errors, confidence intervals or $p$-values) invalid. We provide the first formal analysis of the asymptotic properties of an estimator in the causal change attribution setting of budhathoki2021did, which allows us to perform valid inference in large samples (standard errors, confidence intervals and tests).
Causal change attribution is related to the classical problem of treatment effect estimation pearl2009causality,imbens2015causal,peters2017elements in the need to evaluate counterfactual distributions, but there are some fundamental differences. Within the treatment effects literature, the closest problem to ours is that of causal mediation. tchetgen2012semiparametric gave multiple-robustness results for mediation analysis with a single mediator (explanatory variable). The case with multiple mediators and an arbitrary DAG is substantially more challenging. Although causal mediation with multiple mediators has been studied before (e.g., daniel2015causal), we provide the first multiply-robust estimator in this setting. We also believe to be the first to explicitly build a connection between the causal change attribution problem and the causal mediation literature.
In this section we describe our multiply-robust methodology for causal change attribution. (ref) defines the causal model and the causal change attribution problem and its connections and contrasts with treatment effect estimation. (ref) gives an identification result in a simple example, which we then generalize in (ref). (ref) describe how to implement our estimator, and (ref) gives large-sample inference results. Finally, in (ref) we explain how our method can be incorporated within existing frameworks of causal change attribution, such as Shapley values. We provide the technical conditions and proofs for all lemmas and theorems in Appendix (ref).
We observe multiple i.i.d. measurements of the same variables $(T, \bm{X}, Y)$. Each observation consists of a sample indicator $T \in \{0, 1\}$, $K$ explanatory variables $\bm{X} := (X_1, \ldots, X_K)$ and an outcome of interest $Y \in \mathcal{Y} \subset \mathbb{R}$. The explanatory variables in $\bm{X}$ could be continuous, categorical, or even unstructured data types such as text or images.
We assume that the distribution of $(\bm{X}, Y) \mid T = t$ has a causal Markov factorization spirtes2000causation:
where $\mathrm{PA}_k$ are the parents (direct causes) of $X_k$ in the underlying causal DAG, other than $T$. Throughout, the superscript $(t)$ denotes conditioning on $T = t \in \{0, 1\}$. Each conditional distribution on the right-hand side of (ref) is called a causal mechanism peters2017elements. Without loss of generality, we assume that the explanatory variables are labeled so that $k < k'$ if $X_k \in \mathrm{PA}_{k'}$, i.e., the causal antecedents of $X_{k'}$ in the DAG have indices lower than $k'$. In practice, the researcher only needs to know one such causal ordering, rather than the full causal graph for $(\bm{X}, Y)$. Our causal model implies a distribution for $Y \mid T = t$ by marginalization, i.e., \[P_{Y}^{(t)} := \int P_{Y \mid \bm{X}}^{(t)} \prod_{k=1}^{K} \mathrm{d}P_{X_k \mid \mathrm{PA}_k}^{(t)}.\] The problem of change attribution tries to quantify how much of the differences between $P_{Y}^{(1)}$ and $P_{Y}^{(0)}$ are due to shifting each causal mechanism from its distribution at $T = 0$ to its distribution at $T = 1$. To clarify, intervening on a causal mechanism $P_{X_k \mid \mathrm{PA}_k}$ may impact the outcome both directly and indirectly through other explanatory variables which are causal descendants of $X_k$. We are interested on the total effect of changing $P_{X_k \mid \mathrm{PA}_k}$, without fixing the marginal distribution of its causal descendants. Finally, note that we allow the causal mechanisms to differ between samples, but the underlying DAG is assumed to be the same.
We want to draw a clear distinction between attribution and a related problem that has received a much larger attention in the literature, Average Treatment Effect (ATE) estimation under unconfoundedness. We depict the difference between the causal models underlying these two problems in (ref). In an ATE setting, $\bm{X}$ typically represents pre-treatment covariates, that affect both the outcome $Y$ and the propensity to receive a binary treatment $T$. To compute the causal effect of $T$ we need to control for those covariates, that is, compare units with similar values of $\bm{X}$. In contrast, causal change attribution acknowledges that both the distributions of $\bm{X}$ and of $Y$ given $\bm{X}$ vary depending on whether $T = 0$ or $T = 1$. The goal of causal change attribution is to quantify how much of the effect of $T$ on $Y$ goes through (is mediated by) each causal mechanism in this distribution. Under an unconfoundedness assumption about the assignment of $T$, the problem of causal attribution can be thought of as decomposing the ATE of $T$ on $Y$ into the Natural Direct Effect and the Natural Indirect Effect corresponding to each causal mechanism (pearl2009causality; daniel2015causal). We describe this in more detail, including the corresponding structural equation model and potential outcomes, in (ref).
We note, however, that our definition of attribution does not require $T$ to be a treatment that can be administered in the “interventional” sense. As discussed, $T$ could also be an indicator for group membership (e.g., female or male), time period (e.g., before and after a certain date), or geographical location (e.g., East vs. West Coast of the US). Under which conditions can this change attribution be interpreted as causal? We will assume that the researcher knows a causal ordering of the underlying true causal graph for $(\bm{X}, Y)$. This will allow us to quantify the contribution of each causal mechanism, rather than changes in the marginal distributions of explanatory variables. In practice, researchers can use a combination of causal discovery methods (e.g., the conditional independence test of zhang2011kernel) and domain knowledge to obtain the DAG. Finally, we also assume that there is no unobserved variable $U$ such that $T \rightarrow U$, $U \rightarrow \bm{X}$ and $U \rightarrow Y$. We discuss ways to relax this assumption in future research in (ref).
We begin with a simple yet illustrative example. Consider the simplest situation of $K = 1$. How much of the difference in means $\mathrm{E}[Y \mid T = 1] - \mathrm{E}[Y \mid T = 0]$ is due to changing $P_{X}$? How much of it is due to changing $P_{Y \mid X}$?
In order to answer these questions, we would like to know what the mean of $Y$ would be if we shifted $P_{X}$ to be as in sample 1, but left $P_{Y \mid X}$ unchanged as in sample 0, which we denote as: \[\theta^{\langle 1, 0 \rangle} = \int y \mathrm{d}P^{\langle 1,0 \rangle}_{Y}(y), \quad \text{for} \quad P^{\langle 1,0 \rangle}_{Y} = \int P_{Y\mid X}^{(0)} \mathrm{d} P_{X}^{(1)}. \] The fundamental challenge is that $P^{\langle 1,0 \rangle}_{Y}$ is a counterfactual distribution pearl2009causality, in the sense that we don't observe data sampled from it, and so we cannot estimate $\theta^{\langle1,0\rangle}$ directly as a sample average. It is still possible, however, to identify $\theta^{\langle1,0\rangle}$, as shown in the following lemma:
We show the intuition for (ref) graphically in (ref). Equation (ref) gives identification by regression. We learn the dependence between $Y$ and $X$ in sample 0 through a (non-parametric) regression, and then average that regression function over the $X$ in sample 1. This is essentially a non-parametric generalization of the Oaxaca-Blinder decomposition oaxaca1973male, blinder1973wage. It also appears in the causality literature as part of the mediation formula (see (ref) and pearl2009causality). Equation (ref) gives identification by re-weighting. We average $Y$ in sample 0, but we give more weight to observations whose $X$ is more likely to come from sample 1. Formally, the weights $\alpha(X)$ are Radon-Nykodim (RN) derivatives billingsley1995probability. When $X$ has a density or a probability mass function, these are simply the ratio of densities or probability mass functions between the two samples, respectively. The re-weighting idea has been used in the literature on covariate shift problems shimodaira2000improving, bickel2009discriminative, but we are not aware of causal change attribution methods that leverage this insight.
The following result combines regression and re-weighting methods to obtain a more robust identifying equation. It is essentially a version without pre-treatment covariates of the efficient influence function given in tchetgen2012semiparametric for causal mediation analysis with a single mediator. We restate it here in our notation, because it will make the intuition of our novel result for a general causal graph ((ref)) clearer.
Equation (ref) is doubly robust in the following sense: it still identifies the parameter of interest even if one of the regression function or the weights is misspecified. The term $\mathrm{E}_{(0)}[a(X) e(X, Y)]$ is a “debiasing” term, as in the double/debiased machine learning literature chernozhukov2018double, consisting of an average of the non-parametric regression error $e(Y, X)$ weighted by $a(X)$. If the regression function is correctly specified, this average will be zero regardless of the weights (by the Law of Iterated Expectations). On the other hand, if the regression function is not correctly specified but the weights are, the second term will account and correct for the misspecification of $g(X)$. We refer the reader to the proof in (ref) for details.
We are now ready to introduce the main result. Let $\bm{c} := \langle c_1, \ldots, c_K, c_{K+1} \rangle \in \{0,1\}^{K+1}$ denote a change vector, where $c_k = 1$ if we shift the $k$-th causal mechanism to be as in sample 1, and $c_k = 0$ otherwise. The last entry $c_{K+1}$ indicates whether we want to shift the conditional distribution of the outcome, $P_{Y \mid \bm{X}}$. We will denote by $P^{\bm{c}}_{Y}$ the distribution: \[P^{\bm{c}}_{Y} := \int P_{Y \mid \bm X}^{(c_{K+1})}\prod_{k=1}^K \mathrm{d}P^{(c_k)}_{X_k \mid \mathrm{PA}_k}.\] For example, in (ref) we considered $P^{\langle 1, 0 \rangle}_{Y}$, where we shifted $P_X$ to be as in sample 1 but kept $P_{Y \mid X}$ to be as in sample 0. Of all the possible $P^{\bm{c}}_Y$ for $\bm{c} \in \{0,1\}^{K+1}$, only two are observed directly in the data: $P^{(0)}_{Y} := P^{\langle 0, \ldots, 0, 0\rangle}_{Y}$ and $P^{(1)}_{Y} := P^{\langle 1, \ldots, 1, 1\rangle}_{Y}$; the rest are counterfactual distributions.
We will quantify differences in the distribution of the outcome through the change in some functional (summary measure) of the form $\theta(P_Y) = \int h(y) \mathrm{d}P_{Y}(y)$ for some $h \colon \mathbb{R} \to \mathbb{R}$. Important examples of such functionals are the mean $h(y) = y$, the second moment $h(y) = y^2$ (which, combined with the mean, allows to obtain the variance), and the CDF at a point $u$: $h_u(y) = \mathds{1}\{y \leq u\}$ (which can be inverted to obtain quantiles and Wasserstein distances between two distributions). Our main results concern identification, estimation and inference on $\theta^{\bm{c}} : = \theta(P_Y^{\bm c})$ for a counterfactual distribution described by change vector $\bm{c} \in \{0,1\}^{K+1}$.
To state our main results, we adopt the following notation. We denote with a bar $\Bar{\bm{X}}_{k} := (X_1, \ldots, X_k)$ for any $k \leq K$. In a slight abuse of notation, we write $\Bar{\bm{X}}_{K+1} = (\bm{X}, Y)$ and define $g_{K+1}(\Bar{\bm{X}}_{K+1}) := h(Y)$ as a convention to make mathematical expressions more compact.
The intuition for this result is similar to (ref). The parameter of interest can be estimated by regression, but now we need to use sequentially nested regressions $\gamma_k(\Bar{\bm{X}}_k)$ (i.e. regressions where the outcome itself is a regression function). This fixes the desired distribution for each causal mechanism. Since we are estimating $K$ regression functions, we require $K$ debiasing terms, which appear additively. Each of these debiasing terms includes a “weight” with true value $\alpha_k(\Bar{\bm{X}}_k)$, which is a product of the RN derivatives for the distributions that are shifted with respect to $c_{k+1}$. Including these $K$ debiasing terms makes equation (ref) multiply robust: for each causal mechanism, only the corresponding regression function or the corresponding weights need to be correctly specified, but not necessarily both. We also note that one of our conditions for multiple robustness is that $g_k(\Bar{\bm{X}}_k) = \gamma_k(\Bar{\bm{X}}_k)$ where $\gamma_k(\Bar{\bm{X}}_k) := \mathrm{E}_{(c_{k+1})} [g_{k+1}(\Bar{\bm{X}}_{k+1}) \mid \Bar{\bm{X}}_k]$. This is weaker than setting $\gamma_k(\Bar{\bm{X}}_k) = \mathrm{E}_{(c_{k+1})} [\gamma_{k+1}(\Bar{\bm{X}}_{k+1}) \mid \Bar{\bm{X}}_k]$, because we do not require that $g_{k+1}$ is correctly specified, as long as $a_{k+1}(\Bar{\bm{X}}_{k+1}) = \alpha_{k+1}(\Bar{\bm{X}}_{k+1})$.
In this section, we describe how to implement the multiply-robust estimating equation (ref) in practice.
\paragraph{Step 1. Estimate the Regressions} The nested regression functions $\gamma_k(\Bar{\bm{X}}_k)$ can be estimated by any parametric or non-parametric ML regression algorithm, including LASSO, ridge, random forests, neural networks, boosting, or any ensemble method combining multiple of these. We work recursively, starting with an estimator $\hat \gamma_K(\Bar{\bm{X}}_K)$ for the regression of $h(Y)$ on $\Bar{\bm{X}}_K := \bm{X}$, using the the data from sample $c_{K+1} \in \{0,1\}$ (i.e., the distribution we want to fix for $P_{Y \mid \bm X}$). Next, we obtain $\hat \gamma_{K-1}(\Bar{\bm{X}}_{K-1})$ by regressing $\hat \gamma_K(\Bar{\bm{X}}_K)$ on $\Bar{\bm{X}}_{K-1}$, using the the data from sample $c_{K} \in \{0,1\}$ (i.e., the distribution we want to fix for $P_{X_K \mid \mathrm{PA}_K}$). We proceed this way up until we obtain $\hat \gamma_1(X_1)$.
\paragraph{Step 2. Estimate the Weights} We propose two different approaches to estimate the RN derivatives that appear in the weights $\alpha_k(\Bar{\bm{X}}_k)$: classification and automatic estimation. Both methods target the weights directly, rather than estimating a density or probability mass function in each sample and then taking the ratio. For conciseness, we only describe estimation of $\mu(\Bar{\bm{X}}_j) := \mathrm{d}P^{(1)}_{\Bar{\bm{X}}_j}/\mathrm{d}P^{(0)}_{\Bar{\bm{X}}_j}(\Bar{\bm{X}}_j)$. The reciprocal RN derivative can be estimated analogously by exchanging the indices 0 and 1. Conditional RN derivatives for $X_j \mid \Bar{\bm{X}}_{j-1}$ can be obtained by dividing the RN derivative for $\Bar{\bm{X}}_j$ by the RN derivative for $\Bar{\bm{X}}_{j-1}$.
\paragraph{2a. Classification} By Bayes' rule, we can express:
where $\beta(\Bar{\bm{X}}_j) := \Pr(T = 1 \mid \Bar{\bm{X}}_j)$ and $p := \Pr(T = 1)$ (the posterior and prior probabilities of $\Bar{\bm{X}}_j$ coming from sample 1, respectively). For a given $p$, this defines a one-to-one mapping between $\mu(\Bar{\bm{X}}_j)$ and $\beta(\Bar{\bm{X}}_j)$.
This suggests the following methodology to estimate $\mu(\Bar{\bm{X}}_j)$. First, we train a classification ML algorithm (logistic regression with LASSO or ridge penalties, neural networks, random forests, etc.) to predict $T$ based on the concatenated data for $\Bar{\bm{X}}_j$. Second, we replace $\beta(\Bar{\bm{X}}_j)$ in Bayes' rule with the predicted posterior probabilities, and we replace $(1-p)/p$ by its empirical analog, $n_0/n_1$, where $n_t$ is the number of observations in the $T = t$ sample. The posterior probabilities can also be calibrated by cross-validation, as discussed, e.g., in niculescu2005predicting.
Although this method of estimating RN derivatives is not new (c.f., for instance, sugiyama2012density, or arbour2021permutation), the application to building a multiply-robust moment function is novel.
\paragraph{2b. Automatic Estimation} We specialize general results derived by chernozhukov2021automatic to the causal change attribution problem. In particular, we can characterize the RN derivative $\mu(\Bar{\bm{X}}_j)$ as the solution of: \[\mu(\Bar{\bm{X}}_j) = \arg\min_m \mathrm{E}_{(0)}[m(\Bar{\bm{X}}_j)^2] - 2\mathrm{E}_{(1)}[m(\Bar{\bm{X}}_j)]. \] The literature refers to this method as “automatic” estimation of the weights because the loss function above depends only in $m(\Bar{\bm{X}}_j)$, without using the explicit form of $\mu(\Bar{\bm{X}}_j)$ as an RN derivative. Thus, we could minimize a sample version of the criterion above over some class of functions (e.g. linear with LASSO or ridge penalties, neural networks, random forests, etc.).
\paragraph{Step 3. Estimate the Target Parameter} Finally, we build an estimator of the counterfactual parameter $\theta^{\bm{c}}$ by replacing $g(\Bar{\bm{X}}_k)$ with $\hat\gamma_k(\Bar{\bm{X}}_k)$ and $a(\Bar{\bm{X}}_k)$ with $\hat\alpha_k(\Bar{\bm{X}}_k)$ in a sample analog of equation (ref):
where $\hat{e}_k(\Bar{\bm{X}}_{i,k+1}) : = \hat\gamma_{k+1}(\Bar{\bm{X}}_{k+1}) - \hat\gamma_{k}(\Bar{\bm{X}}_k)$ for $k = 1, \ldots, K$, and we use $\widehat{\mathrm{E}}_{(t)}[f(W)] = n_t^{-1} \sum_{i : T_i = t} f(W_i)$ to denote the sample average on sample $t \in \{0, 1\}$.
In this section, we show consistency and asymptotic normality of our estimator $\hat\theta^{\bm c}$, and explain how to perform large-sample inference on the true parameter $\theta^{\bm c}$. First, we notice that the estimator in (ref) is the sum of two averages over different samples: \[\hat\theta^{\bm c} = \widehat{\mathrm{E}}_{(0)}[\hat\psi_{(0)}(\bar{\bm{X}}_{K+1})] + \widehat{\mathrm{E}}_{(1)}[\hat\psi_{(1)}(\bar{\bm{X}}_{K+1})],\] where $\hat\psi_{(t)}(\bar{\bm{X}}_{K+1})$ contains the terms that depend on the data from sample $t \in \{0,1\}$ (the exact expression is given in the proof of (ref)). This characterization will be useful in the results below in terms of writing the asymptotic variance of $\hat\theta^{\bm c}$. We define an “oracle” version of the same object, $\psi_{(t)}(\bar{\bm{X}}_{K+1})$, by replacing $\hat\gamma_k(\bar{\bm{X}}_k), \hat\alpha_k(\bar{\bm{X}}_k)$ with their respective true values $\gamma_k(\bar{\bm{X}}_k), \hat\alpha_k(\bar{\bm{X}}_k)$ for $k = 1, \ldots, K$.
The following theorem gives an asymptotic normality result for $\hat\theta^{\bm c}$, which allow us to perform valid inference on the true parameter $\theta^{\bm c}$ in large samples (e.g., provide standard errors or confidence intervals):
We now discuss two ways to use the $\theta^{\bm c}$ for the problem of causal attribution: Shapley values and along a causal path.
\paragraph{Shapley Values} Shapley values compute the effect of shifting a causal mechanism averaging over each possible combination of which other causal mechanisms to change. Formally, let $\bm e_k$ be a $(K+1)$-vector with a 1 in the $k$-th entry and 0 everywhere else. The Shapley value associated with $P_{X_k \mid \mathrm{PA}_k}$ is: \[\mathrm{SHAP}_k = \sum_{\bm{c} : c_k = 0} \frac{1}{(K+1) {\binom{K}{\sum_j c_j}}} (\theta^{\bm{c} + \bm{e}_k} - \theta^{\bm{c}}).\]
\paragraph{Along a Causal Path} Define, for each $k = 1, \ldots, K+1$, a vector $\bm{b}_k \in \{0,1\}^{K+1}$ with entries $1, \ldots, k$ equal to 1, and the rest equal to 0. In this method, we define the contribution of the $k$-th causal mechanism to be: \[\mathrm{PATH}_k = \theta^{\bm b_k} - \theta^{\bm b_{k-1}}.\] In other words, we define the contribution of the $k$-th causal mechanism as the effect of shifting from $P_{X_k \mid \mathrm{PA}_k}$ fixing its causal antecedents to be as in sample 1, but its causal descendants to be as in sample 0.
\paragraph{How to choose between these two approaches?} As discussed in the literature budhathoki2021did, the effect of changing the one causal mechanism may depend on which other mechanisms have already been changed; Shapley values take that into account, whereas the method along a causal path does not. Whether this is a relevant concern depends on the application: in some cases, the impact of a given explanatory variable may be approximately the same regardless of the distribution of other explanatory variables. On the other hand, causal attribution along a causal path requires estimating $\theta^{\bm c}$ for only $K$ change vectors $\bm{c}$ (and, by the simplification in (ref), each requires training only one regressor and one classifier), whereas computing Shapley values requires estimating $\theta^{\bm c}$ for $2^{K+1}$ change vectors $\bm{c}$ (each of which could require training up to $2K$ machine learners). Some computationally efficient approximations to $\mathrm{SHAP}_k$ exist (e.g, kolpaczki2023approximating).
Let $\widehat{\mathrm{SHAP}}_k$ and $\widehat{\mathrm{PATH}}_k$ be estimators of $\mathrm{SHAP}_k$ and $\mathrm{PATH}_k$, respectively, that replace the true $\theta^{\bm c}$ with the multiply-robust estimator $\hat{\theta}^{\bm c}$ in (ref). The following corollary gives a useful result for large-sample inference on the causal attribution parameters:
This is a consequence of $\widehat{\mathrm{SHAP}}_k$ and $\widehat{\mathrm{PATH}}_k$ being finite linear combinations of $\hat{\theta}^{\bm c}$. The asymptotic variance could be computed explicitly, although this might be cumbersome, especially for Shapley values, since it requires the $2^{K+1}\times 2^{K+1}$ covariance matrix of $\hat{\theta}^{\bm c}$ for $\bm c \in \{0,1\}^{K+1}$. As an alternative, the multiplier bootstrap of belloni2017program can be used. We describe the procedure in (ref). An advantage of the multiplier approach is that each bootstrap replication does not require re-training the ML estimators of the regression function and weights.
In this section, we discuss two sets of Monte-Carlo simulation results and a real-data application of our method to understanding the gender wage gap. The code for the simulations and the empirical application is available at \url{https://github.com/victor5as/mr_causal_attribution}.
\paragraph{Design 1} We now present simulation evidence on the performance of our method. Our first simulation design is based on the causal model of (ref), with DAG $X_1 \rightarrow X_2$, $X_1 \rightarrow Y$ and $X_2 \rightarrow Y$. We give details about the data-generating process in (ref).
(ref) presents the results of our simulation (over 1,000 draws) for different choices of learners for the regression and the weights. (ref) shows the Mean Absolute Error (MAE) for a scenario where both the regression and the weights are correctly specified parametric learners: OLS with quadratic terms for the regressions and logistic classifier with quadratic terms to build the weights.\footnote{This is the correct specification for the weights, since the log-likelihood ratio between two Gaussian distributions with possibly different variances is a quadratic polynomial.} Comparing our multiply-robust (MR) estimator to using regression or re-weighting only, we can see that its performance is statistically indistinguishable from the best of the two methods (regression).
Tables (ref) and (ref) show what happens when either the weights or the regressions are misspecified, while the specification of the other component of the model is still correct. As a particular form of misspecification, we omit the quadratic terms from either the OLS regression or the logistic classifier. In either case, we observe that MR performs as well as the correctly-specified method, while the misspecified method exhibits much higher errors. This illustrates the multiple robustness property of our estimator. When both the weights and the regression are misspecified parametric models, as in (ref), our theoretical guarantees no longer apply. However, we still show these results for completeness. Although it seems to be the case that the error in Shapley values achieved by the MR methods is about as low as the best of the two other methods, it is significantly larger than the correctly-specified benchmark of (ref).
Next, we turn our attention to flexible, non-parametric learners for the unknown regression functions and weights, to mimic a more realistic situation in which we do not want to make strong parametric assumptions about these. First, we consider neural networks ((ref)). We use MLPRegressor and MLP\-Classifier from the Python library sk\-learn. We use the default settings for all hyperparameters, except that we allow for early stopping and use 3 hidden layers of 100 neurons each. Second, we consider gradient boosting ((ref)), using Gradient\-Boosting\-Regressor and Gradient\-Boosting\-Classifier from sk\-learn with the default settings. In both cases, we calibrate the predicted probabilities after classification with \texttt{CalibratedClassifierCV} also from \texttt{sk\-learn}. The results clearly illustrate (ref): the bias of the MR non-parametric estimator will, in general, vanish faster than the convergence rates of regression-only or re-weighting-only methods. As such, we see that in general MR has equal or lower error than the best of the other two methods. Moreover, the MAE is often statistically indistinguishable from that of the correctly-specified parametric benchmark (ref).
\paragraph{Design 2} In a second set of experiments, we consider the effect of increasing hte number of explanatory variables $K$. We consider a line DAG, $X_1 \rightarrow X_2 \rightarrow \cdots \rightarrow X_K \rightarrow Y$. We give more details about the data-generating process in (ref). We use LASSO to learn the regression functions, and Logistic Regression with $\ell_1$ penalty (Logit-LASSO) for the weights. In both cases, the penalty is selected by cross-validation.
(ref) shows the average (across 500 simulations) of the worst-case (across parameters) absolute error, that is, the average of $\max_{k = 1, \ldots, K+1} |\mathrm{PATH}_k - \widehat{\mathrm{PATH}}_k|$. The pattern is almost identical if we look at the Mean Absolute Error instead. The MR estimator performs significantly better than the other two methods, even though the regression and the weights are correctly specified (i.e., they are truly linear/logistic and sparse).
Following budhathoki2021did, we ask how much of the gender wage gap can be attributed to differences in education or occupation. We use data from the Current Population Survey (CPS) 2015. After applying the same sample restrictions as chernozhukov2018sorted, the resulting sample contains 18,137 male and 14,382 female employed individuals. On this data, the (unconditional) gender wage gap in this sample is of \$7.91/hour (standard error 0.01). In other words, female workers earn 24% less on average than male workers; the difference is statistically significant.
We assume the same causal graphical model as budhathoki2021did, which is captured by the DAG in (ref). We randomly split the data into a training set used to estimate the regression function and weights (40%), a calibration set to calibrate the predicted probabilities from our classification algorithm (10%), and an evaluation set used to obtain the main estimates of the $\hat\theta^{\bm c}$ parameters (50%). To estimate the regression and the weights, we use Hist\-Gradient\-Boosting\-Regressor and Hist\-Gradient\-Boosting\-Classifier from sklearn in Python. We calibrate the probabilities by isotonic regression.
We plot our estimated Shapley values are in (ref). The contribution of education is estimated to be positive and significant. In contrast, the distribution of occupation given education has a significant negative effect, slightly larger in magnitude than the contribution of education, but of opposite sign. Most of the gender wage gap is therefore explained by differences in the distribution of wages given education and occupation. This is consistent with previous findings in the literature. For example, chernozhukov2018distribution, documents the same fact for the UK.
One could think of $P(\text{wage} \mid \text{occup, educ})$ as “residual” variation, i.e., the remaining gender gap which is not explained by education or occupation. It may capture the effect of additional variables that are not included in our analysis. For example, there may be gender differences in labor market experience due to women temporarily exiting the labor force to take care of young children. Unfortunately, experience is not available in census data, and so we cannot quantify its impact, but excluding it from the analysis does not bias our results, because it is plausibly not a causal antecedent of education and occupational choice. It may also be due to other factors, including gender discrimination in promotion or gender differences in wage bargaining.
To gain more insight into the results of our attribution exercise, it is natural to ask two questions: (1) how different is the distribution of education and occupation between the two groups, and (2) how do these differences translate into differences in average wage. We present some summary statistics related to these questions in the appendix. (ref) compares the distribution of educational attainment for males and females. It appears that females tend to be slightly more educated, with higher rates of college and advanced degree graduation. This could explain the positive Shapley value associated with education; for comparison, college graduates earn \$12.60/hour more than high school graduates on average. (ref) shows the distribution of occupation conditional on having a college degree. Females are more predominant in administrative, education or healthcare professions, whereas males are more likely to work in management and sales. For comparison, managers earn \$16.66/hour more than educators and \$6.29/hour more than healthcare practitioners on average. At the same time, there are wide differences across genders that cannot be explained by education and occupation alone: female college graduate managers earn \$13.77/hour less on average than their male counterparts, which is consistent with the finding that $P(\text{wage} \mid \text{occup, educ})$ has the largest Shapley value in our causal change attribution exercise.
In this paper, we develop a new estimator for causal change attribution measures, which combines regression and re-weighting methods in a multiply-robust estimating equation. We provide the first results on consistency and asymptotic normality of ML-based estimators in this setting, and discuss how to perform inference. Moreover, our method can be used to estimate Shapley values, which will inherit the consistency and large-sample distribution properties.
Finally, we suggest one direction for future research. The causal interpretation of the change attribution parameters relies on an assumption of no unobserved confounding. The sensitivity bounds in Chernozhukov2022long could be adapted as a way to test the robustness of a causal change attribution study to unobserved confounding.
The authors would like to thank Victor Chernozhukov, Dominik Janzing, Eric Tchetgen-Tchetgen, and seminar audiences at Amazon and MIT for providing helpful feedback. We also would like to thank Patrick Blöbaum for help with code integration in DoWhy.
This paper presents work whose goal is to advance the fields of Causal Inference and Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
The “causal” interpretation of our change attribution measures is grounded on certain assumptions about the Directed Acyclic Graph (DAG) that codifies causal dependence among the variables in the data. We are precise about these assumptions ((ref)), and discuss some ways to relax them in future extensions ((ref)). Our theoretical results also rely on regularity conditions, which we state formally in the Appendix and explain in more intuitive terms in the main body of the paper. Practitioners need to be mindful that, if these assumptions do not hold, our proposed algorithm is not guaranteed to lead to accurate results.