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.
65,581 characters · 21 sections · 48 citation commands
Causal Machine Learning for Moderation Effects
\begingroup \let\relax \endgroup
\thispagestyle{empty} \setcounter{page}{1}
Detecting and interpreting heterogeneity in treatment effects is crucial for understanding the impact of interventions and (policy) decisions. Researchers have recently developed many methods to estimate heterogeneous treatment effects. However, there are still limitations in interpreting these effects. For example, suppose that the effect of a particular training program for the unemployed is larger for women than for men. By comparing the average treatment effects for these two groups, without taking into account the different distribution of other covariates of men and women, such as education or labor market experience, we might implicitly compare a group with more extended labor market experience (men) with a group with shorter labor market experience (women). Therefore, obtaining a balanced distribution of relevant characteristics across different groups may be crucial to ensure proper comparisons and draw meaningful conclusions. In the specific example, we might want to ensure that both groups have the same average years of labor market experience. This approach allows us to isolate the differences in causal effects between groups in a way that is not confounded by certain other covariates. However, the difference may still be due to differences in unobservable characteristics that vary between the two groups of interest. Determining whether a particular variable causes differences in treatment effects is often of interest. In the literature, these variables are called causal moderator variables. In addition to balancing distributions of all covariates that confound the moderation effect across groups, additional assumptions must hold to interpret the differences in treatment effects causally, i.e., so that the group variable can be considered an unconfounded moderator.
\defcitealias{Chernozhukov:2018}{Chernozhukov, Chetverikov, Demirer, Duflo, Hansen, Newey, & Robins, 2018}
This paper discusses how to estimate differences in treatment effects between groups in an unconfoundedness setting. First, a new parameter called balanced group average treatment effect (BGATE) is introduced. This parameter is a group average treatment effect (GATE) with a specific distribution of a-priori-determined covariates. It is beneficial for comparing the treatment effects of two groups with each other. This results in the difference of two BGATEs called $\Delta$BGATE. We demonstrate how it relates to a difference of two group average treatment effects ($\Delta$GATE), discuss its identification and propose different estimators for discrete moderators (for simplicity, we refer to the variable for which we find heterogeneous treatment effects as moderator) and discrete treatments. For the estimator based on double/debiased machine learning (subsequently abbreviated as DML) \citepalias{Chernozhukov:2018}, we show that it is asymptotically normal centred at the true value, which allows valid inference. Additionally, we show how the $\Delta$BGATEs can be estimated using automatic debiased machine learning (Auto-DML) and a reweighting method. DML relies on using a double-robust score function that depends on the propensity score, which can be problematic if the propensity score is extreme Lechner:2024,Busso:2014,Frolich:2004. Auto-DML uses Riesz representers instead of the propensity score; hence, it does not suffer from the problem of extreme propensity scores Chernozhukov:2022a,Chernozhukov:2018b. The reweighting method allows to estimate the $\Delta$BGATE with whatever method the user prefers for the GATE estimation, as the $\Delta$BGATE equals the $\Delta$GATE on the reweighted data. A small-scale simulation study and an empirical example using administrative labor market data from Switzerland demonstrate the practicability of the different estimators.
The paper proceeds as follows: In Section (ref), we define the parameter of interest for the case of a binary treatment and a binary moderator and discuss its interpretation and identification under unconfoundedness. Section (ref) proposes different estimation strategies and shows the asymptotic properties of the DML estimator. Section (ref) depicts the design and results of a small-scale Monte Carlo simulation. In Section (ref), we demonstrate how the new parameter can be used in practice and Section (ref) concludes. Appendix (ref) shows how the $\Delta$BGATE can be interpreted in a causal way. The identification and asymptotic proofs are depicted in Appendix (ref). The different algorithms for the estimation are explained in Appendix (ref). Appendix (ref) shows more details about the simulation study, and Appendix (ref) provides more details about the empirical example.
Several authors have recently contributed to the topic of effect heterogeneity \citep*[e.g.,][]{Tian:2014, Wager:2018, Athey:2019, Kuenzel:2019, Nie:2021, Semenova:2021, Knaus:2022b, DiFrancesco:2022, Foster:2023, Kennedy:2023}. The proposed methods make it possible to estimate heterogeneities between fine-grained subgroups more accurately. For a recent review of the various methods and their performance, see Knaus:2021. As a result, many applied papers now use such methods to detect heterogeneities \citep*[e.g.,][]{Davis:2017, Knaus:2022, Cockx:2023}. However, decision-makers are often more interested in heterogeneities for a small subset of covariates than at the finest possible granularity. Therefore, several papers show how to detect (low-dimensional) heterogeneities at the group level, called “Group Average Treatment Effects” (GATEs) (a GATE is a conditional average treatment effect (CATE) with a small number of conditional variables, often only a single one). Several approaches have been developed to estimate GATEs. Abrevaya:2015 show how to identify GATEs nonparametrically under unconfoundedness and estimate them using inverse probability weighting estimators. Athey:2019 and Lechner:2018 modify a random forest algorithm to adjust for confounding to estimate heterogeneous treatment effects. Semenova:2021 use the DML framework to find heterogeneity based on linear models. Zimmert:2019 and Fan:2022 develop a two-step estimator that allows estimating GATEs nonparametrically. The first step is estimated using machine learning methods, and in the second step, they apply a nonparametric local constant regression.
We add to this literature by proposing a new parameter of interest: the difference between two GATEs with (partially) balanced characteristics. By taking the difference, we can disentangle the difference in treatment effects from the difference in the distribution of covariates between the two groups. Hence, this paper is related to the Kitagawa-Oaxaca-Blinder decomposition (KOB), often used to decompose differences in outcomes between two groups into one part that is due to the difference in the distribution of other covariates and another part that is due to the variable of interest Kitagawa:1955, Oaxaca:1973, Blinder:1973. Vafa:2024 show that such decompositions suffer from omitted variable bias if the model for estimating them is not complex enough. They suggest a methodology for wage gap decompositions using foundation models, such as large language models. Yu:2023 introduce a new decomposition method that allows to identify how a treatment variable contributes to a group disparity in an outcome variable. They decompose it into different prevalence of treatment, different treatment effects across groups, and different patterns of selection into treatment, and they show how the decomposition can be estimated non-parametrically. Similarly, Chernozhukov:2013 show how to model counterfactual distributions based on regression methods. They consider two scenarios: making changes to either the distribution of covariates associated with the outcome or the conditional distribution of the outcome given the covariates. The $\Delta$BGATE fits into the second scenario. Additionally, Chernozhukov:2023 introduce automatic debiased machine learning for covariate shifts. This approach is related as the $\Delta$BGATE can be interpreted as a $\Delta$GATE on a population with a shifted covariate distribution. While Chernozhukov:2023 consider the case of two distinct datasets that do not overlap and are statistically independent, we have one dataset and shift the distribution of covariates within this dataset.
Chernozhukov:2018b consider the problem of interpreting heterogeneous treatment effects from another angle. They suggest sorting the estimated partial effects by percentiles and comparing the covariate means of observations falling in the different percentiles to see which individuals with which characteristics are most positively and negatively affected. In contrast, we are in this paper not interested in finding the characteristics of the individuals with the most positive and negative effects but in comparing the treatment effects of the two groups while balancing the distribution of some covariates. If the goal is to to learn the characteristics of (relative) winners and losers from a treatment, their approach suits well. However, if the goal is to understand how much a variable drives the differences in treatment effects between two groups, our approach is more appropriate.
\defcitealias{Chernozhukov:2018}{Chernozhukov et al. (2018)}
As has already become apparent from the discussion of the first part of the relevant literature, the DML literature is closely related to our work. \citetalias{Chernozhukov:2018} developed this framework, which allows using machine learning methods for causal analysis. Machine learning algorithms may introduce two biases: a regularization and an overfitting bias. The main idea of DML is that by using Neyman-orthogonal score functions, we can overcome the regularization bias, and by using cross-fitting, an efficient form of sample splitting, we can overcome the overfitting bias. Many papers have adapted the general DML framework for different settings, for example, for continuous treatments \citep*[e.g.,][]{Kennedy:2017, Semenova:2021}, for mediation analysis \citep*{Farbmacher:2022}, for panel data \citep*[e.g.][]{Clarke:2023} or for difference-in-difference estimation \citep*[e.g.,][]{Zimmert:2018, SantAnna:2020}. We contribute to this literature by using this highly flexible framework to estimate the new parameter of interest.
Last, we add to the literature on moderation effects. Discussions of moderation effects are primarily found in the psychology and political science literature \citep*[e.g.,][]{Gogineni:1995, Frazier:2004, Fairchild:2010, Marsh:2013, Bansak:2021a, Blackwell:2021}. Common approaches to analyzing moderation effects specify interactions in a regression or subgroup analysis. Without additional assumptions, these parameters cannot be interpreted causally. Bansak:2021 studies causal moderation effects in an experimental setting by showing what identifying assumptions are needed. Because of the randomization of treatment, it is possible to estimate the causal effect of the moderator on the outcome separately for the treated and control units and then subtract both estimates to obtain the moderator effect. Such differences cannot be interpreted causally if the moderator variable influences some covariates, which is often the case. In addition, Bansak:2022 have show how to identify and estimate subgroup differences in a regression discontinuity design.
The causal moderation framework used in this paper is based on the potential outcome framework of Rubin:1974. A causal effect is defined as the difference between two potential outcomes, whereas for a unit, we only observe one of these potential outcomes. Therefore, finding a credible counterfactual is problematic. We observe $N$ i.i.d. observations of the independent random variables $H_i = (D_i, Y_i, Z_i, X_i)$ according to an unknown probability distribution $\mathds{P}$. Here, the focus is on a treatment $D_i$ and a moderator $Z_i$ which, for simplicity, are assumed to be binary (realizations of the treatment variable are $d \in \{0,1\}$, and of the moderator variable $z \in \{0,1\}$). The formal theory presented in Appendix (ref) is based on fixed numbers of discrete values for $Z_i$ and $D_i$. As usual, the potential outcomes are indexed by the treatment variable: $(Y_i^{0}, Y_i^{1})$. Finally, a set of $k \in \{1, \dots, p\}$ covariates $X_{i,k}$ might simultaneously affect treatment allocation and potential outcomes, where $X_i = (X_{i,1}, \dots, X_{i,p})$. Additionally, potential covariates and potential moderators are defined as: $(X_i^{0}, X_i^{1}, Z_i^{0}, Z_i^{1})$.
Since only realizations of one of the potential outcomes are observed, we can never consistently estimate realizations of the individual treatment effect (ITE), $\xi_i = Y_i^1- Y_i^0$. However, under suitable assumptions, the identification of, for example, the average treatment effect (ATE) $\theta = \mathrm{E}[Y_i^1- Y_i^0]$ is possible Imbens:2009. It is often interesting to additionally investigate different aspects of the heterogeneity of the $\xi_i$ which can be captured by so-called conditional average treatment effects (CATE). A CATE measures the average treatment effect conditional on a (sub-) set of covariates $X_i$. The individualized average treatment effect (IATE) and the group average treatment effect (GATE) are specific CATEs. The IATE measures the treatment effect at the most granular aggregation level. Namely, it compares the average effect of the treatment for all individuals with a specific value of all relevant covariates used. Formally, the IATE is defined as follows:
The GATE measures the treatment effect at the group level, i.e. at a more aggregated level than the IATE, but still at a finer level than the ATE. Formally, the GATE is defined as follows:
As long as the interest lies only in describing effect heterogeneity, IATEs and GATEs are sufficient. However, if the interest lies in the difference in treatment effects between the two groups, the difference between the two GATEs ($\Delta$GATE), i.e.
may be difficult to interpret because the two groups may differ in the distribution of other covariates $X_i$.
Thus, a new parameter is introduced, the balanced group average treatment effect (BGATE). This parameter is called BGATE, as it is designed to balance the distribution of other variables within the groups (defined by the different values of $Z_i$) we want to compare to each other. The variables used to balance the GATEs are denoted as $W_i$. $W_i$ is part of $X_i$. If $W_i$ is empty, or $W_i$ is independent of $Z_i$, the BGATE reduces to the GATE. Thus, the new parameter of interest, denoted by $\theta^B(z)$, is defined as
and its difference, $\theta^{\Delta B}$, as
A $\Delta$BGATE ($\theta^{\Delta B}$) represents the difference between two groups, adjusting the distribution of some other covariates ($W_i$) in both groups to the overall population distribution. The $\Delta$BGATE usually shows associational moderation effects. Under certain assumptions, we define a causal balanced group average treatment effect ($\Delta$CBGATE) that can be interpreted causally. The discussion of this parameter is referred to Appendix (ref).
To clarify the distinction between a $\Delta$GATE and a $\Delta$BGATE, the $\Delta$GATE is decomposed into two components: the $\Delta$BGATE, representing the direct effect of the moderator variable, and the compositional effect, arising from differences in the distributions of $W_i$ across the groups. The compositional components capture the differences in the distribution of $W_i$ in both $Z_i$ subsamples weighted by their relative importance:
The derivation of this decomposition is shown in Appendix (ref). This decomposition is similar to the KOB-decomposition Kitagawa:1955, Oaxaca:1973, Blinder:1973 used to decompose differences in outcomes between two groups into one part that is due to the difference in the distribution of other covariates and one part that is due to the variable of interest. It is often used in the context of analyzing gender wage differences Blau:2017.
The $\Delta$BGATE is similar to the direct effect in the KOB-decomposition. However, in a KOB-decomposition, the distribution of the balancing variables $W_i$ is not adjusted to the overall population distribution but to the distribution of a specific group. If the distribution of $W_i$ does not differ across groups, the direct effect in the KOB-decomposition and the $\Delta$BGATE are equal. In that case, the $\Delta$BGATE also equals the $\Delta$GATE. Adjusting the distribution of the covariates to the overall population distribution is intuitive if we are interested in the population and not in a specific part of the population.
To identify the GATE, BGATE, $\Delta$GATE or $\Delta$BGATE in an unconfoundedness setting, usual identifying assumptions are needed Imbens:2004:
For the proof of Lemma (ref) see Appendix (ref).
\defcitealias{Chernozhukov:2018}{(e.g. Chernozhukov et al., 2018, 2022)} Since we are interested in differences in treatment effects, the estimation strategy focuses on the $\Delta$BGATE. The estimator can easily be adapted to estimate the BGATE and can be estimated using different estimators, from causal machine learning \citetalias[e.g.][]{Chernozhukov:2018} to other semiparametric estimators Ai:2007, Ai:2012. In this paper, we explain three methods based on causal machine learning and prove analytically that the DML estimator is $\sqrt{N}$-consistent and asymptotically normal.
\defcitealias{Chernozhukov:2018}{Chernozhukov et al. (2018)}
The identification results suggest a three-step estimation strategy. To obtain a flexible estimator that allows for a potentially high-dimensional vector of covariates, the first suggested estimator relies on the methodology of DML as proposed by \citetalias{Chernozhukov:2018}. In the first step, the usual double robust score function is estimated. In the second step, the score function is regressed on the two indicator variables defined by the different values of the moderator variable $Z_i$ and the covariates we want to balance $W_i$. Last, we take the difference between the two groups defined by the moderator and average over the variables we balance. This approach is close to the approach of Kennedy:2023 for estimating CATEs with the DR-learner. However, instead of estimating a CATE, we average over $W_i$.
The estimated doubly robust score function is given by the following expression:
with
$\hat \delta(h)$, $\hat \mu_d(z, x)$, $\hat \lambda_z(w)$ and $\hat \pi_d(z,x)$ denote the estimated values of $\delta(h)$, $ \mu_d(z,x)$, $\lambda_z(w)$ and $ \pi_d(z,x)$, respectively. Furthermore, notice that $\hat g_z(w) = \hat \mathrm{E}[\hat \delta(H_i)| Z_i = z, W_i = w]$ is the regression of $\hat \delta(H_i)$ on $Z_i$ and $W_i$ and that $\tilde g_z(w) = \hat \mathrm{E}[\delta(H_i)| Z_i = z, W_i = w]$ is the corresponding oracle regression of $\delta(H_i)$ on $Z_i$ and $W_i$. Hence, $\hat \mathrm{E} [\dots | \dots]$ denotes a generic regression estimator, which can be linear or non-linear, depending on the presumed data-generating process. Last, the estimated nuisance parameters are $\hat\eta = (\hat\mu_d(z,x)$, $\hat \pi_d(z,x)$, $\hat \lambda_z(w),$ $\hat g_z(w))$.
As explained above, the score function $\hat \delta(h)$ has to be estimated. The product of the nuisance function errors for $\hat \delta(h)$ must converge faster than or equal to $\sqrt N$, and cross-fitting with K-folds ($K > 1$) must be used. In the second estimation step, the product of the nuisance function errors has to converge with $\sqrt N$ and cross-fitting with J-folds ($J > 1$) has to be used. If certain conditions are met, the estimator is $\sqrt N$-consistent and asymptotically normal (see Subsection (ref)). Because $E[\phi^{\Delta B}(H_i; \theta^{\Delta B}, \eta)] = 0$, the variance of $\hat \theta^{\Delta B}$ is given by
and is estimated by $\widehat{\operatorname{Var}}(\hat \theta^{\Delta B}) = \frac{1}{N}\sum_{j = 1}^J \sum_{i \in S_j}[ \hat \phi^{\Delta B}(H_i; \theta^{\Delta B}, \hat\eta)]^2$ with $S_j$ being a random fold in the second estimation step. Algorithm (ref) in Appendix (ref) depicts the implementation of this DML estimator.
We investigate the asymptotic properties of the estimator and impose the following assumptions:
\defcitealias{Chernozhukov:2018}{Chernozhukov et al., 2018} Assumptions (ref) to (ref) made are standard in the DML literature \citepalias{Chernozhukov:2018}. The only difference is that these assumptions are applied for the first and the second estimation step. Assumption (ref) is needed to ensure that the product of the nuisance function errors converges faster than or equal to $\sqrt{N}$. $L_2$ convergence rates of various machine learning methods adhere to these properties when sparsity conditions are met. For example, Belloni:2013 shows that under approximate sparsity, meaning that the sorted absolute values of the coefficients decay quickly, the error of the Lasso estimator is of order $O\left(\sqrt{\frac{s \log(\max(p,N))}{N}}\right)$ with $p$ being the regressors and $s$ the unknown number of true coefficients in the oracle model. Similar rates also depending on the sparsity level, have been shown for shallow regression trees or honest and arbitrarily deep regression forests Wager:2016, Syrgkanis:2020, boosting in sparse linear models Luo:2016 or a class of deep neural nets Farrell:2021. For detailed information on how sparsity conditions depend on the parameters of the predictors, we refer to the cited references. Assumption (ref) is needed because, in the second estimation step, we regress an already estimated quantity $\hat \delta$ on $Z_i$ and $W_i$. Stability can be perceived as a type of stochastic equicontinuity for a nonparametric regression. Kennedy:2023 proves that linear smoothers, such as linear regression, random forests or nearest neighbour matching, are stable.
Given these assumptions, we can derive the main theoretical result:
Theorem (ref) states that the estimator is $\sqrt N$-consistent and asymptotically normal. See Appendix (ref) for the proof.
The previous section shows that DML relies on a double-robust score function that depends on the propensity score. As extreme propensity scores can be problematic and distort the results Lechner:2024, Busso:2014, Frolich:2004, we propose a second estimation strategy that does not rely on the propensity score. We use the Auto-DML framework Chernozhukov:2022a, Chernozhukov:2022b that relies on using Riesz representers instead of propensity scores. Hence, we proceed in three steps, similarly to the DML estimation strategy outlined above. The usual double robust score function using the Riesz representer is estimated in the first step. In the second step, the score function is regressed on the two indicator variables defined by the different values of the moderator variable $Z_i$ and the covariates we want to balance $W_i$. Last, we take the difference between the two groups defined by the moderator and average over the distribution of $W_i$ in the population.
The estimated double robust score function is given by the following expression:
with
Again, the score function $\hat \delta(h)$ has to be estimated, the product of the nuisance function errors for $\hat \delta(h)$ must converge faster than or equal to $\sqrt N$, and cross-fitting with K-folds ($K > 1$) must be used. Similarly, in the second estimation step, the product of the nuisance function errors has to converge with $\sqrt N$ and cross-fitting with J-folds ($J > 1$) has to be used.
The variance of $\hat \theta^{\Delta B}_{Riesz}$ is given by
and is estimated by $\widehat{\operatorname{Var}}(\hat \theta^{\Delta B}_{Riesz}) = \frac{1}{N}\sum_{j = 1}^J \sum_{i \in S_j}[ \hat \phi^{\Delta B}_{Riesz}(H_i; \theta^{\Delta B}_{Riesz}, \hat\eta)]^2$. Algorithm (ref) in Appendix (ref) depicts the implementation of the Auto-DML estimator. The proof of the Auto-DML estimator for the BGATE is beyond the scope of this paper.
Last, we suggest an estimation strategy independent of the specific method a researcher wants to use to estimate GATEs. The idea is to change the data so that the distribution of the balancing variables is the same across the groups defined by the moderator variable. After having balanced the data, the BGATE can be estimated using any method for estimating GATEs as the GATE on the balanced data equals the BGATE.
An example of a reweighting strategy is as follows: In the first step, each observation is duplicated such that we have an identical observation for each group defined by the moderator variable. In the second step, a nearest-neighbour matching procedure is performed. For each observation in the expanded dataset, the covariate values for the balancing variables $W_i$ are compared with those of all other observations that share the same treatment assignment. Using the Mahalanobis distance metric, each observation is matched to the "nearest" observation from the original dataset with similar covariate values. This ensures that the matched pairs are comparable in terms of observed covariates. When $W_i$ is a vector of several variables, we adjust for their covariance structure using the inverse covariance matrix (i.e., using the so-called Mahalanobis distance for matching). The covariate values of the matched observations are then assigned to the duplicated observation in the reference dataset. This substitution allows each observation to represent a counterfactual scenario where it maintains similar covariate characteristics but with the opposite treatment assignment. The exact procedure for reweighting the data is depicted in Algorithm (ref) in Appendix (ref).
After reweighting the data, any estimator for $\Delta$GATE can be employed to estimate the $\Delta$BGATE (e.g., the DML-based version outlined as Algorithm (ref) in Online Appendix (ref)). However, the variance estimator must be adjusted to account for the different weights each observation receives. The derivation of the variance is shown in Online Appendix (ref). While the formal proof of the asymptotic properties of the reweighting strategy is beyond the scope of this paper, it performs well in the simulations.
We start with simulating a p-dimensional covariate matrix $X_{i,p}$ with p=10. The first two covariates are drawn from a uniform distribution $X_{i,0}, X_{i,1} \sim \mathcal{U}[0,1]$ and the remaining covariates from a normal distribution $X_{i,2} \dots, X_{i,p-1} \sim \mathcal{N}\left(0.5, \sqrt{1/12}\right)$. All covariates have a mean of 0.5 and a standard deviation of $\sqrt{1/12}$. The moderator variable $Z_i$ is drawn from a Bernoulli distribution with probability $P(Z_i = 1|X_{i,0}, X_{i,1}) = (0.1 + 0.8 \beta(X_{i,0} \times X_{i,1};2,4))$. $\beta(X_{i,0} \times X_{i,1};2,4)$ denotes the cdf of a beta distribution with the shape parameters $a = 2$ and $b = 4$. The propensity score is created similarly as in Kuenzel:2019 and Wager:2018. The treatment variable $D_i$ is drawn from a Bernoulli distribution with probability $P(D_i = 1|X_{i,0}, X_{i,1}, X_{i,2}, X_{i,5}, Z_i) = \left(0.2 + 0.6 \beta( \frac{X_{i,0} + X_{i,1} + X_{i,2} + X_{i,5} + Z_i}{5}; 2, 4)\right)$.
The response functions under treatment and non-treatment, and the two states of the moderator variable are specified. The highly non-linear non-treatment response function is specified similarly as in Nie:2021 and is given by
The response functions under treatment depend on $Z_i$ and are defined as:
They are chosen such that the $\Delta$BGATE is different from the $\Delta$GATE. Last, we simulate the potential outcomes as $Y_i^{d}(z) = \mu_{d}(z,X_i) + e_{i,d,z}$ $\forall z \in \{0,1\}$ with noise $e_{i,d,z} \sim \mathcal{N}(0,1)$. Summing up, the data consists of an observable quadruple $(y_{i,r}, d_{i,r}, z_{i,r}, x_{i,r} )$.
The parameters of interest are two $\Delta$BGATEs and a $\Delta$GATE, namely:
$X_{i,0}$ is unbalanced across the two moderator groups, whereas $X_{i,2}$ is balanced. Hence, $\theta^{\Delta B}_{X_0}$ is different from $\theta^{\Delta G}$, whereas $\theta^{\Delta B}_{X_2} $ is equal to $\theta^{\Delta G}$. We generate $R = 2,000$ samples of size $N = 1,000$, $R = 1,000$ samples of size $N = 2,500$, $R = 500$ samples of size $N = 5,000$, and $R = 250$ samples of size $N = 10,000$. Since the variance is doubled when the sample size is halved, we make the number of replications proportional to the sample size (to reduce the computational costs of the simulations). The true values are estimated on a sample with $N = 1,000,000$. The parameters of interest are estimated by the three different estimation strategies outlined above, i.e. DML, Auto-DML, and the reweighting strategy.
The algorithm used for the DML estimation strategy is depicted in Algorithm (ref) in Appendix (ref). We use two folds in both estimation steps ($K = 2$, $J = 2$). All nuisance functions in the DML estimation strategy are estimated using random forests (number of trees: 1000). As shown by Bach:2024, tuning the learners used to estimate the nuisance parameters is crucial. Table (ref) in Appendix (ref) shows the hyperparameters used.
The algorithm used for the Auto-DML estimation strategy is depicted in Algorithm (ref) in Appendix (ref), and we use again two folds in the first and two folds in the second estimation step ($K = 2$, $J = 2$). The algorithm is implemented by using a neural net as proposed in Chernozhukov:2022c. More details about implementing the specific neural net can be found in Appendix (ref).
Last, the effects are estimated by first reweighting the data so that the distribution of the balancing variables is the same across the groups defined by the moderator variable. The algorithm used for reweighting is depicted in Algorithm (ref) in Appendix (ref). After having balanced the data, we use the DML estimation strategy to estimate the $\Delta$GATE with the adjusted variance estimator as explained in Subsection (ref) (Algorithm (ref) in Appendix (ref)).
Table (ref) shows the simulation results for the different sample sizes, parameters of interest, and estimators. Results with additional performance measures can be found in Appendix (ref).
Comparing the different estimators shows that DML has the smallest root mean squared error (rmse) in almost all cases, followed by Auto-DML. As expected, the reweighting strategy often leads to a slightly higher standard deviation (std), which results in a higher rmse. The coverage for Auto-DML is sometimes worse, as the bias of the effect or the standard error is too large. It is possible that the performance of the Auto-DML can be improved by tuning the neural networks for this particular estimation task. Such an additional investigation is, however, left to future research. Additionally, the sample size for using a neural net should be large enough, which might not be the case for a sample size of 1,250. This might be a reason why the coverage for Auto-DML is poor for the estimation of the $\Delta$GATE (78.3%). Concluding, the DML estimator performs best in the simulations. However, this finding may not generalize to cases where DML is known to have performance issues, such as when propensity scores become extreme. In such cases the suggested Auto-DML estimator or the reweighting method might be more appropriate.
Comparing the results across different sample sizes, the std and rmse approximately halve by increasing the sample size by four, which suggests a $\sqrt{N}$ convergence of the estimator already for the comparatively small sample size used in the simulations. This pattern can be observed for all estimators.
Next, we compare the different parameters of interest. Figure (ref) depicts the distributions of the errors for $\hat \theta^{\Delta B}_{X_0}$ and $ \hat \theta^{\Delta G}$ if the effect of interest is $ \theta^{\Delta B}_{X_0}$. Estimating $\hat \theta^{\Delta G}$ leads to a different result, as the variable $X_{i,0}$ is not balanced across the two groups of $Z_i$.
In contrast, $X_{i,2}$ is already balanced across the two groups of $Z_i$. Figure (ref) depicts the distributions of the errors for $\hat \theta^{\Delta B}_{X_2}$ and $\hat \theta^{\Delta G}$ with the effect of interest being $\theta^{\Delta B}_{X_2}$. As expected, estimating $\hat \theta^{\Delta G}$ leads to the same result as estimating $\hat \theta^{\Delta B}_{X_2}$ because $X_{i,2}$ is already balanced. Hence, if the covariate(s) is (are) not balanced, it is important to differentiate between the two effects.
An empirical example is taken from the literature that evaluates the effects of active labour market policies (ALMP) on unemployed individuals. Card:2018 summarize the estimates of more than 200 studies on the effects of ALMPs in a meta analysis. Generally, they find higher positive effects for females, long-term unemployed, and low-income individuals. Several recent studies use causal machine learning to analyze the heterogeneous impacts of such ALMPs. For example, Cockx:2023 examine Belgian unemployed individuals, finding positive medium-term program effects, especially for recent immigrants with low local language proficiency. Burlat:2024 reports that positive impacts of technical training in Eastern France vary by education. Knaus:2022 use Swiss administrative data from 2003 to evaluate ALMPs with various causal machine learning methods. They find effect heterogeneity in the first six months after the start of the job search programs. The heterogeneity relates to labor market characteristics and nationality. Individuals with disadvantaged labor market characteristics, like low-income or low-labor market attachment, benefit more from the programs. Similarly, foreigners benefit more. These heterogeneous effects can be explained by the indirect costs of the programs (due to not searching intensively for a job during a program and therefore needing more time to find a job), which are lower for more disadvantaged and foreign individuals.
\defcitealias{Knaus:2020}{Lechner, Knaus, Huber, Frölich, Behncke, Mellace & Strittmatter (2020)}
To explore the method in an ALMP setting, we use the publicly available dataset from \citetalias{Knaus:2020}. Using this adminiatrative data, we study the effect of job search programs on the employment status of Swiss unemployed in 2003. Knaus:2022 provide an extensive description of the data. The dataset provides information on individuals' participation in a job search program, where $(d = 1)$ indicates participation and $(d = 0)$ indicates non-participation. Two outcome variables are considered: a short-term outcome, $(Y_i^{short})$, which represents the number of months an individual was employed during the first six months following the start of the program, and a medium-term outcome, $(Y_i^{long})$, which represents the number of months employed during the last six months available in the dataset (months 28 to 33 after the program's start). The effect is claimed to be identified in an unconfoundedness setting. Using the same specification as Knaus:2022, we include several covariates on the individuals' socioeconomic background and the labor market history $(X_i)$.
A concern regarding the treatment definition arises because caseworkers can assign individuals to the program anytime during their unemployment spell. Therefore, individuals with more favourable labor market characteristics might be overrepresented in the control group, as they already found a job by the time a caseworker would have assigned them to the program. To overcome this issue, we follow Knaus:2022 and predict (pseudo) treatment start dates for each individual in the control group. To ensure consistency in treatment definitions between participants and non-participants, we restrict the analysis to individuals who remain unemployed at their (pseudo) treatment start dates. The final sample consists of 84,582 unemployed individuals. We use this dataset to illustrate the proposed method and check whether the heterogeneities found by Knaus:2022 are due to the variables identified by the authors or whether other underlying variables possibly confound them.
For the sake of brevity, we only consider effect heterogeneity concerning nationality and employability. The employability variable reflects the caseworkers' assessment of whether the unemployed individual is easy or difficult to reemploy. The first part of Table (ref) shows selected covariate means for Swiss and non-Swiss individuals among the treated and non-treated individuals. The second part of the table shows the same for easy and hard to employ individuals among the treated and non-treated individuals. This helps to better understand which variables might account for the variation in treatment effects of the two groups. More descriptives can be found in Appendix (ref).
There are some differences in covariates related to the previous labor market history, namely in past income, previous job and qualifications and the number of unemployment spells in the last two years. Furthermore, some sociodemographic characteristics, such as being married, also differ. Finally, as expected, there are more foreign than Swiss individuals with a mother tongue other than German, French, Italian or Raeto-Romansh. Similarly, hard to employ individuals are more likely to not have a Swiss mother tongue. Other variables, such as age or gender, are already well balanced, so balancing these covariates should not change the effect much. Nevertheless, we include age and gender to avoid that they will be no longer balanced after balancing some (previously unbalanced) other covariates (since they might correlate with the newly balanced covariates).
For the empirical application we focus on the DML estimator, as it performs best in the simulation study and we have asymptotic properties for it. The first two columns of Table (ref) show the different effects considered in the analysis. These effects include the $\Delta$GATE, a $\Delta$BGATE with already balanced covariates, such as age and gender, a $\Delta$BGATE with additionally adding marital status, an extended $\Delta$BGATE balancing additionally unbalanced covariates like past annual income, previous job, and qualification variables. Then, we further add mother tongue to the variables to be balanced. Finally, the analysis considers a $\Delta$BGATE that balances all covariates included in the study.
Columns three to six of Table (ref) depict the results for the different effects in the first six months, a period that is commonly called “lock-in” period. As a reference point, the average treatment effect is $\hat \theta_{lock-in} = -0.785$ ($0.021$). $\hat \theta^{\Delta G}$ shows that the difference in treatment effects between Swiss and non-Swiss individuals is statistically and economically significant. Hence, it seems that the program works better for foreigners. As pointed out above, the interpretation is not straightforward because the two groups have unbalanced covariates. After explicitely balancing already balanced sociodemographic characteristics like age, gender, and also marital status the coefficient remains relatively stable. When balancing the labor market history, including covariates like past income and previous job details, there is a notable reduction in the coefficient. This directly translates into an increase in the compositional effect (column 6) that shows the part of $\hat \theta^{\Delta G}$ that comes from differences in the distribution of covariates. As anticipated, additionally, balancing mother tongue results in an even lower difference between Swiss and non-Swiss individuals of only $0.088$, which is not statistically significant at the 5% significance level. Balancing all covariates included in the analysis does not reduce the difference further. Hence, it becomes evident that mother tongue disparities and differences in the labor market history significantly contribute to the variance in treatment effects between Swiss and non-Swiss individuals.
Similarly, we find that the program has a significantly different effect on easily employable individuals compared to hard to employ individuals (columns 7 to 10). It works better for hard to employ individuals and the $\hat \theta^{\Delta G}$ is $-0.371$ ($0.055$). Again, balancing sociodemographic characteristics does only slightly change the result. Balancing the labor market history reduces the difference in treatment effects to $-0.184$ ($0.067$), which is halving the difference. Interestingly, the difference does not reduce to zero after balancing the labor market history. Hence, this employability variable created by the caseworkers does not seem to be based only on the labor market history. Balancing mother tongue does not change the result further, however, balancing for all remaining covariates included in the analysis does reduce the difference to basically zero.
The average medium-term effect is $\hat \theta_{medium} = 0.012$ ($0.033$), which is very small and not significantly different from zero. We also do not find any relevant differences between Swiss and non-Swiss individuals in the medium-term. However, there is a difference between easily and hard to employ individuals. It seems that even in the medium-term, the program works better for hard to employ individuals. This difference becomes even more pronounced when balancing other covariates. This directly relates to the findings of the meta-analysis of Card:2018, where the job search program appears more effective for such individuals.
Concluding, these results highlight the importance of carefully interpreting group average treatment effects, as the difference in treatment effects between some groups diminish or vanish after balancing the distribution of other covariates.
This paper presents a novel approach for analyzing and interpreting treatment heterogeneity in an unconfoundedness setting. We introduce a parameter called $\Delta$BGATE for measuring the difference in treatment effects between different groups while accounting for variations in covariates. This paper proposes an estimator based on DML for discrete treatments and moderators, demonstrating its consistency and asymptotic normality under standard conditions. Additionally, we outline two alternative estimation strategies. The first one is the so-called Auto-DML, which is a DML-type estimator with the additional advantage that the already well-documented non-robustness of DML estimators to extreme estimated propensity scores can be avoided. The second alternative estimator is a reweighting method that allows the use of any estimator that is consistent for the GATE applied to the reweighted data. A simulation study shows the practicability of these estimation strategies. An empirical example illustrates the proposed estimand and underlines that seeming causal heterogeneity may be caused by an underlying different distribution of other covariates. The proposed new parameter allows a more informative interpretation of heterogeneity and, thus, a better understanding of the differential impact of decisions. Future research could extend the estimation approach to continuous treatments and moderators. This paper shows identification in an unconfoundedness setting. It would be interesting to extend it to an instrumental variable setting for the treatment, the moderator, or both. More extensive research on how to tune the RieszNet and providing a thorough analysis of the asymptotic properties of Auto-DML and the proposed reweighting estimator are fruitful areas for further research. More extensive simulation studies would lead to a more comprehensive picture of the finite sample properties of the proposed estimators.