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.
86,101 characters · 20 sections · 28 citation commands
Effect or Treatment Heterogeneity? Policy Evaluation with Aggregated and Disaggregated Treatments
\doublespacing
The analysis of causal effects is at the heart of empirical research in economics, political science, the biomedical sciences, and beyond. To evaluate and design policies, interventions, or programs for units with different background characteristics, it is crucial to develop a thorough understanding of the heterogeneity present in causal relationships. There is now a large literature that develops and applies identification and estimation strategies for causal or treatment parameters that explicitly take into account such heterogeneity, see \citeA{Athey2017} or \citeA{Abadie2018EconometricEvaluation} for recent overviews.
Most attention is on effect heterogeneity of binary treatments, while less is given to treatment heterogeneity. However, many binary treatments in applications can be conceived as heterogeneous in the sense that they summarize (many) underlying effective treatments that impact the outcome of interest. In such cases it is not clear whether effect heterogeneity as defined in the canonical binary treatment setting reflects heterogeneous effects or heterogeneity in the effective treatments. This paper proposes new estimands to disentangle these sources of heterogeneity in a general setting where the analyzed binary indicator does not coincide with the effective treatments. The distinction between sources of heterogeneity is crucial for evaluating and improving assignment mechanisms. Consider the following two scenarios: \footnote{See also Supplementary Appendix (ref) for a numerical example.}
Scenario 1 (binarized treatments): Multiple or continuous treatments are ex-post subsumed into a binary indicator (e.g. different smoking intensities become “smoking yes/no”). Such aggregations are often motivated by simplicity or data availability, but can have unintended consequences: First, discovered effect heterogeneity can be a spurious byproduct of aggregation and thus falsely be attributed to unit background characteristics. Second, actual effect heterogeneity could be masked as a consequence of the aggregation.
Scenario 2 (multiple treatment versions): A binary treatment takes different versions after assignment, e.g. access to a training program with multiple specializations. Here, effect heterogeneity could result from different version targeting and not from different effectiveness of the versions themselves. This distinction is crucial for policy makers to assess the quality of the version assignment mechanism.
In this paper we propose a novel method for decomposing canonical effect heterogeneity into new estimands that are representative of (i) heterogeneous effects and (ii) heterogeneity from different underlying treatment compositions. These decomposition parameters serve as summary measures to evaluate the consequences of (dis)aggregating treatment variables for the causal analysis. Furthermore they provide a simple framework for comparing the quality of treatment version assignments and their heterogeneity across units or groups.
We develop a simple but flexible nonparametric method for estimation and statistical inference for the decomposition parameters. Our framework allows for the use of machine learning techniques such as random forests, deep neural networks, or high-dimensional sparse regression models in the estimation of the nuisance parameters. We provide high-level conditions regarding the required rates for machine learners, their interaction with the nonparametric decomposition step, and the number of effective treatments $J$. We also provide sufficient conditions for explicit example estimators.
The decomposition can be used to conduct hypothesis tests that consider all effective treatments simultaneously. This allows to test selection and effect heterogeneity without the need for multiple testing procedures. It compares favorably to conventional multi-valued treatment effect analysis under many effective treatments $J\rightarrow\infty$, expanding sets of nuisance parameters, and extreme propensity scores. In particular, regular inference is still achievable even if propensity scores are arbitrarily close to zero (limited overlap). This result is obtained by leveraging local superefficiency properties of probability estimators. The large sample theory extends to other parameters that combine unbiased signals with machine learning inputs and estimated weights. Monte Carlo simulations suggest that coverage rates are close to nominal in finite samples.
We provide two applications of our decomposition method, one for each leading scenario: First, we show that parts of the finding that the detrimental effect of smoking on birth weight is largest for white mothers can be explained by white mothers smoking more heavily conditional on being smokers. Similarly, different effects for different age groups are partly due to teenage mothers smoking less intensively compared to older mothers. Second, we investigate the lower effectiveness of access to the Job Corps training program for women compared to men. We find evidence that the well-documented gender gap is largely explained by the vocational training curriculum, which focuses more on lower paying service jobs for women and more on higher paying craft jobs for men. Imposing the same mix of vocational training as part of our decomposition removes 73% of the total gender differences in the effect on earnings.
The paper is structured as follows: Section (ref) discusses the related literature. Section (ref) outlines the decomposition of the causal effect parameters and discusses their identification. Section (ref) contains the estimation and inference method. Section (ref) introduces the technical assumptions and discusses the large sample properties. Section (ref) provides the Monte Carlo study. Section (ref) contains the application. Section (ref) concludes. We also provide an \href{https://github.com/MCKnaus/causalDML/blob/master/R/HK_decomposition.R}{implementation in R} and replication notebooks.\footnote{For Section (ref) see \href{https://mcknaus.github.io/assets/code/Replication_NB_smoking.nb.html}{mcknaus.github.io/assets/code/Replication_NB_smoking.nb.html} and for Section (ref) \href{https://mcknaus.github.io/assets/code/Replication_NB_JC.nb.html}{mcknaus.github.io/assets/code/Replication_NB_JC.nb.html} on GitHub.}
The proposed decomposition complements the literature that considers (dis)aggregated binary treatments. \citeA{Lechner2002ProgramPolicies} discusses how to aggregate average effects of multiple treatments into composite treatment effects. \citeA{Hotz2005PredictingLocations} and \citeA{Hotz2006EvaluatingProgram} investigate the consequences of summarizing different training components in one binary indicator and emphasize the potential lack of external validity under latent treatment heterogeneity. \citeA{McCall2016Government-SponsoredAdults} discuss the challenges to determine the optimal degree of coarsening of multi-valued treatments in applications. Similarly, a recent stream of papers formalizes structural causal models and interpretations of compound treatments Cole2009TheInference,VanderWeele2009ConcerningInference,Hernan2011CompoundInference,Petersen2011CompoundGraphs. \citeA{VanderWeele2013CausalTreatment} note that non-homogeneous treatments violate the second component of the “Stable Unit Treatment Value Assumption” \cite<SUTVA,>{Rubin1980RandomizationComment}: no-multiple-versions-of-treatment, which requires a homogeneous treatment or at least the treatment variation irrelevance assumption of \citeA{VanderWeele2009ConcerningInference}. \citeA{VanderWeele2013CausalTreatment} formalize a setting where this assumption is violated and provide several new identification results and estimands. Aggregating heterogeneous treatments has also been discussed in the context of instrumental variables Angrist1995Two-stageIntensity,Marshall2016CoarseningEstimates,Andresen2021Instrument-basedRestriction,Harris2022InterpretingEducation, regression discontinuity designs Cattaneo2016InterpretingCutoffs, and models with spillovers and interactions Manski2013IdentificationInteractions,Vazquez-Bare2022IdentificationExperiments. These papers mostly discuss the consequences of (dis)aggregation of treatments on unconditional estimands and their connection to (weighted) causal effects. Our paper focuses on the consequences of (dis)aggregation on effect heterogeneity.
The focus on effect heterogeneity is motivated by the surging literature that develops \cite<e.g.>{Athey2016,Athey2017a,Kunzel2017,Knaus2021} and applies \cite<e.g.>{Davis2020RethinkingJobs,Knaus2022HeterogeneousApproach,Buhl-Wiggers2022SomeInterventionb} flexible machine learning methods to the estimation of heterogeneous causal effects. We build on the double/debiased machine learning framework by \citeA{Chernozhukov2018}. They use Neyman-orthogonal score functions and sample splitting in conjunction with machine learning methods for estimation of low-dimensional parameters that depend on nuisance quantities.
Regarding heterogeneity analysis, there is now a series of papers that obtain (functional) parameters by localizing these score functions using (nonparametric) regression or machine learning methods Lee2017,Zimmert2019NonparametricConfounding,Colangelo2020DoubleTreatments,Kennedy2020OptimalEffects,Semenova2021DebiasedFunctions,Fan2022EstimationData,Knaus2022DoubleUnconfoundedness,Heiler2022HeterogeneousPolarization. Our theoretical contribution builds on the structural function approach by \citeA{Semenova2021DebiasedFunctions} with least squares series estimation Newey1997ConvergenceEstimators,Belloni2015SomeResults,Cattaneo2020LargeEstimators. We extend some of the inferential results by \citeA{Semenova2021DebiasedFunctions} to settings where pseudo-outcomes are constructed as a weighted average of Neyman-orthogonal scores with (estimated) weights and potentially many treatments.
The paper is also related to the literature regarding inference on effect parameters under extreme propensity scores or “limited overlap” Khan2010IrregularEstimation,Rothe2017RobustOverlap,Ma2020RobustWeighting,Hong2020InferenceOverlap,Heiler2021ValidScores. Limited overlap occurs by construction when allowing for “many treatments” $J\rightarrow\infty$. In this case, the set of nuisance parameters is expanding and classic multi-valued treatment effect parameters \cite<e.g.>{Cattaneo2010EfficientIgnorability} are irregularly identified which complicates inference. The decomposition method, however, always yields three aggregate (functional) parameters independently of $J$. As a consequence, regular estimation and inference regarding heterogeneity is still feasible as long as $J$ does not grow too fast. In finite samples, determining what constitutes a many treatments setup is difficult as $J$ is always a finite number and a small lower bound for propensities are hard to distinguish from a zero lower bound Rothe2017RobustOverlap. Thus, a method that is robust to a potentially large number of treatments provides safeguard for empirical practice.
Assume we observe independent data $(Y_i,D_i,T_i,X_i)$ for $i=1,\dots,n$. $Y_i$ denotes the outcome of interest, $D_i \in \{0,1\}$ is the analyzed binary indicator, $T_i \in \mathcal{T} = \{0,1,\dots,J\}$ indicates the effective treatment\footnote{Note that \citeA{Manski2013IdentificationInteractions} also uses the term “effective treatments” in the context of interference. Like in our setting, it describes the treatments that create variation in potential outcomes. In the following the term “treatment” refers to effective treatment if not stated differently.}, and $X_i$ contains confounding variables. We consider settings that are characterized by two features: (i) Not $D_i$, but the effective treatment $T_i$ has a direct influence on the outcome creating potential outcomes $Y_i(t)$ for each $t \in \mathcal{T}$. Thus, we assume SUTVA with respect to the effective treatment such that $Y_i = \sum_t \mathbbm{1}(T_i =t) Y_i(t)$. (ii) Conditional on $T_i$, the binary indicator $D_i$ is deterministic, i.e. it perfectly separates the support $\mathcal{T}$. We use directed acyclic graphs \cite<DAGs, see e.g.>{Pearl1995CausalResearch} to outline our main scenarios:
Figure (ref) outlines the causal structure of Scenario 1 where the binary indicator variable $D_i$ is the result of an ex-post aggregation and not directly related to the outcome in a structural sense. In practice, this aggregation is often conducted after the outcome realizes, which makes it unlikely for $D_i$ to affect $Y_i$ directly. This is indicated by a missing arrow from $D_i$ to $Y_i$.
The DAG in Figure (ref) depicts the causal structure of Scenario 2 where a randomized binary treatment $D_i$ precedes the confounded allocation of treatment versions $T_i$. Here, $D_i$ is not an ex-post variable with regards to $Y_i$. $Y_i$ and $D_i$ are associated as the latter determines which treatment versions are available, but has no direct effect beyond that. Its effect is completely mediated through the treatment versions $T_i$.
We denote $D_{t,i} = \mathbbm{1}(T_i =t)$ to indicate that unit $i$ is observed in treatment $t$ and define $e_t(x) = P(D_{t,i} = 1 | X_i = x)$ as corresponding propensity score. Without loss of generality, we assume throughout that $T_i=0$ denotes a homogeneous control condition. Thus, the binary indicator is defined as $D_i=\sum_{t\neq 0}D_{t,i}$ and $D_{0,i} = 1-D_i$ in what follows.
We are interested in cases with causal structures as described in Section (ref) but analysis limited to binary $D_i$. Here, typical quantities of interest are conditional average treatment effects ($CATE$) or aggregations thereof like the average treatment effect ($ATE$). Canonical strong ignorability assumptions for $D_i$ are then imposed to exploit quantity $\tau(x) = E[Y_i | D_i = 1, X_i = x] - E[Y_i | D_i = 0, X_i = x]$ for identification of the $CATE$. However, when $D_i = 1$, the potential outcome is not uniquely defined unless $J=1$. Therefore, the question is what does this $\tau(x)$ actually identify? Given the setting outlined in Section (ref), we can backwards engineer the actually identified estimand in terms of potential outcomes of the effective treatment:
Equation (ref) shows that the estimand consists of three components: First, a weighted average of $CATEs$ of the effective treatments, $\tau_t(x) = E[Y_i(t) - Y_i(0) | X_i = x]$, with weights depending on the conditional probability of the respective effective treatment. Second, a weighted average of effective treatment specific selection effects. Third, a selection effect into the control group. The selection effects are positive if units with characteristics $x$ that are actually observed in treatment $t$ show higher potential outcomes than the general population described by $x$, or negative if vice versa. The second and third term is relevant if there is selection into the effective treatments even after conditioning on observed confounders. This can e.g. occur in the case of a randomized binary treatment in Scenario 2 where the selected heterogeneity variables $X_i$ might not include all confounders for the treatment versions.
The decomposition in (ref) highlights that the interpretation of the underlying estimand becomes more nuanced in the presence of heterogeneous treatments. What is supposed to be an easily interpretable $CATE$ depends now on the potentially unknown distribution of effective treatments and selection into those treatments. Thus, without further assumptions, heterogeneous effects attributed to the binary indicator can be driven by different $CATEs$, different compositions of the effective treatments, different selection effects of the effective treatments, or combinations thereof.
The development of an identifiable decomposition for a parameter such as (ref) requires conditional independence or related assumptions. For example, the leading scenarios in Figure (ref) and Figure (ref) imply the same conditional independence relationship between effective treatment and potential outcomes despite not being Markov equivalent:
Condition (ref) implies the more conventional “weak unconfoundedness” assumption for multi-valued treatments \cite<see e.g.>{Cattaneo2010EfficientIgnorability,Yang2016PropensityTreatments}. The latter is sufficient for the decomposition proposed in Section 3.3. Therefore we maintain it throughout the paper together with a common support assumption:
Assumption (ref) is a standard assumption in the multiple treatments setting Imbens2000TheFunctions,Lechner2001. It imposes that (a) the set of conditioning variables contains all confounders and (b) there are comparable units across all treatments. Under Assumption (ref) $E[Y_i(t) | D_{t,i}=1,X_i=x] = E[Y_i(t) | X_i=x]~\forall~t \in \mathcal{T}$ and the selection effects in (ref) disappear. The underlying estimand then simplifies to
We call this estimand the natural conditional average treatment effect $nATE(x)$ because it is the result of the actual or “natural” effective treatment composition. It is important to note that, even under Assumption (ref), the differences between units characterized by $x$ and $x'$ can result from different treatment shares, different treatment $CATEs$, or both. We thus could detect seemingly heterogeneous effects, even if the treatment $CATEs$ are constant within treatments but not homogeneous between treatments, i.e. $\tau_t(x) = \tau_t~\forall~ t \in \mathcal{T}, x \in \mathcal{X}$ but $\tau_t \neq const. ~\forall~ t \in \mathcal{T}$, as long as the probabilities to be observed in the different effective treatments are heterogeneous. This fundamentally affects the interpretation of heterogeneous effects even if the underlying effective treatments are not observable. If they are observable, however, we can further decompose heterogeneous effects of the binary indicator in what follows.
In this section we demonstrate how to disentangle actual effect heterogeneity and heterogeneity driven by selection into effective treatments. We propose to decompose the $nATE(x)$ in two parts:
where $\pi_t=E[D_{t,i}]$ are the unconditional treatment probabilities.\footnote{ In principle, analogous decompositions could be constructed with alternative weights for the effective treatments, e.g. equal weighting $1/J$. However, the unconditional effective treatment probabilities ensure that $nATE(x) = rATE(x)$ in the case of completely randomized effective treatments. This shows resemblance to the comparison between the canonical $ATE$ and $ATET$: $ATE$ is an average effect under hypothetical random assignment, $ATET$ under actual treatment assignment. Both coincide under a completely randomized binary treatment.} The first component on the right hand side fixes the composition of the effective treatments at the population value. It resembles a situation where effective treatments are randomly allocated using the population level selection probabilities. Thus, we refer to it as the random conditional average treatment effect $rATE(x)$. All heterogeneity in $rATE(x)$ is driven by “real” effect heterogeneity within treatments, $\tau_t(x) \neq \tau_t(x')$ for some $x,x'\in\mathcal{X}$, as the underlying treatment composition is held fixed. In other words, differences in $rATE(x)$ describe effect heterogeneity compositionis paribus. Thus, we can exploit potential heterogeneity in $rATE(x)$ to test for classic (or “within”) effect heterogeneity.
The second component of the decomposition $\Delta(x)$ is the part of $nATE(x)$ stemming from the interaction of non-constant effective treatment probabilities and different effective treatments having different effects (“between” treatment effect heterogeneity). The decomposition is redundant, i.e. $\Delta(x) = 0~\forall~x\in \mathcal{X}$, under (i) effective treatment composition homogeneity $\frac{e_t(x)}{\sum_{t \neq 0}e_t(x)} - \frac{\pi_t}{\sum_{t \neq 0}\pi_t} = 0~\forall~t \in \mathcal{T}$ and $x \in \mathcal{X}$, (ii) treatment variation irrelevance $E[Y_i(t) | X_i=x] = E[Y_i(t')| X_i=x]~\forall~x\in \mathcal{X},~ t, t' \in \mathcal{T}$ VanderWeele2009ConcerningInference , or (iii) if positive and negative components net out to zero. Hence, $\Delta(x) \neq 0$ is a necessary condition for unequal treatment probabilities and between treatment effect heterogeneity and thus a violation of SUTVA. Furthermore, heterogeneity in $\Delta(x)$ is a necessary condition for heterogeneous assignment probabilities, within treatment effect heterogeneity, or both. Thus, the decomposition addresses a variety of relevant policy questions. The focus on such necessary conditions offers statistical advantages over testing related conditions in the standard multi-valued treatment effect setup when there are many effective treatments, see Section (ref).
Under Assumption (ref), the conditional average potential outcome of treatment $t$ is identified as $\mu_t(X_i) \equiv E[Y_i(t) | X_i] = E[Y_i(t) | D_{t,i} = 1, X_i] = E[Y_i | D_{t,i} = 1, X_i] $ and accordingly the decomposition terms are identified as:
Aggregations or projections of the three estimands are thus also identified. In particular, let $Z_i = f(X_i)$ denote a (low dimensional) function (e.g. subset) of confounders supported on $\mathcal{Z}$ and define
Focusing on specific subgroups defined by $Z_i = z$ provides concise, predictive summaries of heterogeneity or allocation differences without compromising on the dimensionality of confounders $X_i$ and is standard in the literature on effect heterogeneity Chernozhukov2017GenericExperiments,Semenova2021DebiasedFunctions. The unconditional decomposition terms $nATE = E[nATE(X_i)]$, $rATE = E[rATE(X_i)]$, and $\Delta = E[\Delta(X_i)]$ are special cases thereof.\footnote{The $rATE$ is a special case of composite treatment effects Lechner2002ProgramPolicies. If $J\rightarrow \infty$, it can approximate integrated dose-responses of continuous treatments \cite<e.g.>{Kennedy2017Non-parametricEffects}. $\Delta$ is similar to the population average prescriptive effect in the context of policy learning Imai2021ExperimentalRules.} Thus, we focus on the former throughout the paper. Estimation and inference methods are presented in Section (ref).
The interpretation of $\Delta(x)$ depends on the scenario: In Scenario 1, $\Delta(x)$ and its aggregates have descriptive interpretation. They describe how much of $nATE(x)$ is driven by an underlying effective treatment mix that deviates from the population mix. A non-constant $\Delta(x)$ indicates that the binarization has consequences for detected heterogeneous effects. Thus, it helps to understand heterogeneity resulting from the binarization. In Scenario 2, $\Delta(x)$ and its aggregates provide information for assignment evaluation. Positive (negative) values indicate that assignment of treatment versions is better (worse) than random assuming that individuals act equivalently under the hypothetical random assignment compared to the observational assignment Heckman2020Epilogue:Revisited. A non-constant $\Delta(x)$ indicates that the selection quality of versions varies across different groups. Thus, the estimand provides an evaluation of the actual assignment mechanism.
In this section, we outline a flexible estimation approach for the (conditional) decomposition terms and propose a method for conducting valid statistical inference. The method accommodates the use of modern machine learning and other non- or semiparametric methods in the estimation of the required nuisance parameters.
We propose to approximate the conditional expectations of the decomposition terms $g(z)$ by a linear combination of transformations $b(z)$ of heterogeneity variables $z$, i.e.
where $\beta_0$ is the parameter vector of the best linear predictor given as solution to equation $E[b(Z_i)(g(Z_i) - b(Z_i)'\beta_0)] = 0$. $r_g(z)$ is the approximation error and $b(z)$ can be basis transformations of the regressors of interest such as polynomials, splines, wavelets, or other functions. The number of components in $b(\cdot)$ is allowed to grow with the sample size which allows us to be agnostic about the shape of the true $g$-function.
Let in the following $\eta = \eta(x) = (\mu_0(x),\dots,\mu_J(x),e_0(x),\dots,e_J(x))'$ denote the vector of nuisance quantities and write $\eta = \eta_i = \eta(X_i)$ with subscript and argument suppressed whenever it does not cause confusion. Also define the unconditional selection probability vector $\pi = (\pi_0,\dots,\pi_J)$.
We follow the general idea of \citeA{Semenova2021DebiasedFunctions} to construct “Neyman-orthogonal” scores $\psi_i(\eta,\pi)$ such that $g(z) = E[\psi_i(\eta,\pi)|Z_i=z]$. These scores are defined by having an (approximate) zero Gateaux derivative with respect to the underlying nuisance parameters at the true parameter vector Chernozhukov2018. The robust scores for the three decomposition parameters considered here are weighted combinations of the well-known Neyman-orthogonal scores for average potential outcomes Robins1995SemiparametricData, also known as augmented inverse probability weighting (AIPW) scores:
where $\psi_i^{[t]}({\eta})$ is the score of the treatment $t$ specific average potential outcome and $\Psi_i({\eta})$ is the score for the group described by the binary indicator. Table (ref) shows how to combine these scores to form unbiased signals of the decomposition parameters. These combinations retain Neyman-orthogonality with respect to $\eta$, see Appendix (ref), but inference has to be adjusted for uncertainty in the estimation of $\pi$, see Section (ref).
Consider now the projection of the score functions onto the space spanned by the $k$-dimensional transformation of $Z_i$, $b(Z_i)$. This yields the estimator
where the score of a decomposition term with estimated nuisance quantities $\psi_i(\hat{\eta},\hat{\pi})$ serves as pseudo-outcome in the corresponding least squares regression on $b(Z_i)$. For $\hat{\pi}$ we use simple sample averages, i.e. $\hat{\pi}_t = n^{-1} \sum_{i=1}^n D_{t,i}$. Estimation of $\hat{\eta}$ can be done via modern machine learning methods or other non- and semiparametric estimation methods with good approximation qualities for the functions at hand. For details regarding the technical assumptions, consider Section (ref). We require that all components in $\hat{\eta}$ are obtained via $K$-fold cross-fitting:
Cross-fitting controls the potential bias arising from overfitting using flexible machine learning methods without the need to evaluate the complexity of the function class that contains true and estimated nuisance quantities. If finite parametric models such as linear or logistic are assumed for the nuisance quantities, the proposed methodology can be applied without cross-fitting.
Under suitable assumptions, the predictions using estimator (ref) are consistent for $g(z)$. Moreover, it is possible to conduct asymptotically valid inference, i.e. for any $z_0 = z_{0,n}$ we can construct $(1-\alpha)\%$ confidence intervals for the decomposition parameter as
where $q_{1-\alpha/2}$ denotes the $(1-\alpha/2)$-quantile of the standard normal distribution and $\hat{\Omega}$ is a consistent sample estimator of the asymptotic variance $\Omega$ (see Section (ref) and Appendix (ref)). The estimator explicitly takes into account the additional uncertainty from estimating the unconditional treatment probabilities in the decomposition terms. The interval in (ref) is also valid for the best linear predictor $b(z_0)'\beta_0$ under misspecification if the approximation error is not too large. It provides asymptotically accurate confidence intervals around the true $g$-function if the approximation error vanishes at a suitable rate as the number of basis functions or transformations increases. For the technical details consider Section (ref).
In this section, we present and discuss the large sample properties of the proposed decomposition estimator. First, we introduce the relevant definitions. We then discuss the assumptions required for (i) all decomposition parameters, (ii) $nATE$, and (iii) $rATE/\Delta$ specifically and their connections to the literature. We contrast (ii) and (iii) as the $nATE$ tends to require less restrictive conditions compared to $rATE/\Delta$. We then present the main Theorem and outline potential extensions.
In the following, quantities like $\psi_i()$ or $r_g()$ are used in their generic sense, i.e. for a given choice of decomposition parameter $nATE$, $rATE$ or $\Delta$. $a \lesssim b$ means $a/b = O(1)$ and $a\lesssim_P b$ means $a/b = O_p(1)$. For a general matrix $M$ denote its largest (smallest) eigenvalue by $\lambda_{max}(M)$ ($\lambda_{min}(M)$). Let $\eta \in T$ where $T$ is a convex subset of some normed vector space. Denote the realization set of the estimated nuisance quantities by $\mathcal{H}_n = \mathcal{E}_n \times \mathcal{M}_n \subset T$ with $ \mathcal{E}_n = E_{0,n}\times E_{1,n}\times\dots\times E_{J,n}$ and $ \mathcal{M}_n = M_{0,n}\times M_{1,n}\times\dots\times M_{J,n}$, where $E_{t,n}$ and $M_{t,n}$ are the realization sets that contain estimates $\hat{e}_t(X_i)$ and $\hat{\mu}_t(X_i)$ with probability $1-u_n$. Define their $L_q$ error rates
and the slowest $L_q$ rates over all treatments $s_{n,q} = \sup_{t \neq 0} s_{t,n,q}$ and $m_{n,q} = \sup_{t \neq 0} m_{t,n,q}$. Note that $\psi_i^{[t,0]}(\eta) = \psi_i^{[t]}(\eta) -\psi_i^{[0]}(\eta)$. By definition $\varepsilon_i = \psi_i(\eta,\pi) - E[\psi_i(\eta,\pi)|Z_i]$ where $\psi_i(\eta,\pi)$ corresponds to the score function of the decomposition parameter of choice from Table (ref). Denote $g(z)= E[\psi_i(\eta,\pi)|Z_i=z]$ where $g \in \mathcal{G}$ with $\mathcal{G} = \mathcal{G}_n$ being a function class potentially depending on $n$. Thus, $g(z) = b(z)'\beta_0 + r_g(z)$ where $\beta_0$ is the parameter of the best linear predictor defined as the root of equation $E[b_i(g(Z_i) - b_i'\beta_0)] = 0$ with $b_i = b(Z_i)$ being the $k$-dimensional basis functions. Also define the potential outcome mean error $\varepsilon_i(t) = Y_i(t) - E[Y_i(t)|X_i]$ and its conditional variance $\sigma_t^2(X_i) = E[\varepsilon_i(t)^2|X_i]$. For the $nATE$ machine learning bias components, we define
and equivalently for $rATE/\Delta$ with score functions according to Table (ref). Remainder terms $R_n$ are defined in Appendix (ref). Let $\gamma_t = E[b_i\psi_i^{[t,0]}(\eta)] = E[b_i\tau_t(X_i)]$, $\gamma = (\gamma_1 \dots \gamma_J)$, and define $a_i = (a_i^{[1]} \dots a_i^{[J]})'$ with $a_i^{[t]} = (1-\pi_0)^{-2}(D_{t,i}(1-\pi_0) + D_{0,i}\pi_t - \pi_t)$. Now let $Q = E[b_ib_i']$ and define
We now present the assumptions required for all decomposition parameters. They are meant to hold uniformly over $n$ if not stated otherwise: {Decomposition Assumptions:}
}
A.1 rules out multicollinearity of the basis functions used for the nonparametric heterogeneity analysis in the last stage. A.2 is a mild heterogeneity restriction on the potential outcomes. A.3 is crucial: It is concerned with the degree of overlap for a general number of treatments $J = J_n$. In particular, it assumes that there is strong overlap for the control group and the aggregate treatment, i.e. control and aggregate treatment propensities are uniformly bounded away from zero. However, the propensities for treatments $t=1,\dots,J$ are allowed to be arbitrarily close to zero as long as they vanish at most at a rate proportional to their respective unconditional treatment selection probability $\pi_t$. This allows for limited overlap at each treatment level which is necessary when their number is allowed to increase with the sample size, i.e. $J \rightarrow \infty$. We suggest to assess Assumption A.3 empirically by analyzing the (estimated) distribution of $e_t(x)/\pi_t$ for all $t$: If these re-scaled scores have sufficient density bounded away from zero by the same standard used to assess conventional propensity score distributions Heiler2021ValidScores, then the assumption is likely to hold, see Appendix (ref) and (ref) for examples based on the empirical applications from Section (ref). Assuming a homogeneous $1/J$-rate for all $\pi_t$ is without loss of generality: If the product converges to zero for some $t$, it vanishes from relevant first-order approximations and estimation properties are eventually determined by the treatments that obey Assumption A.3. Moreover, if A.3 only applies to a smaller finite subset of treatments, it effectively corresponds to strong overlap for these particular $t$ and thus estimators behave analogously to standard AIPW for a control potential outcome. Note that the growth of $J$ is restricted to rate $o(n)$ such that consistent estimation of unconditional multi-valued treatment effects is still possible, albeit at a slower rate compared to the strong overlap case similar to \citeA{Hong2020InferenceOverlap}.
A.4 says that the worst relative prediction for the cross-fitted propensities is bounded on the realization set with high probability. This is a non-standard assumption, in particular when $J\rightarrow \infty$. It is likely to hold for frequency based methods, i.e. estimators that use some form of (weighted) average within the cells defined by $D_{t,i}$ for $t=0,\dots,J$ to construct propensities including advanced machine learning methods. A sufficient, but by no means necessary, condition is uniform consistency of $\hat{e}(x)$ over $\mathcal{X}$ at rate $o(e_t(x)^{-1})$. This can be shown to hold e.g. for single-index models, see Theorem 2 and Theorem 3 in \citeA{Ma2022TestingOverlap}, and nonparametric kernel regression Heiler2022NonparametricFrequencies under weak conditions. A key point is that these estimators inherit a local superefficiency property from $\hat{\pi}_t$, i.e. faster convergence rate $|\hat{\pi}_t - \pi_t| \lesssim_P (nJ)^{-1}$ in regimes with many treatments/vanishing unconditional selection probabilities. A.4 then requires the estimators to have a consistency rate increased by a factor of $\sqrt{J}$ compared to the finite $J$ case. For parametric estimators this holds as long as $J = o(n)$ while for kernel regression, for example, under the usual smoothness assumptions with $\mathcal{X}$ of dimension $d$ and a bandwidth $h$, it requires that $\sqrt{\log(n)nh^d/J} = o(1)$. We provide some more intuition about Assumption A.4 and the links between large $J$ and superefficient nuisance parameter estimation in Section (ref).
We now present the assumptions required for $nATE$ followed by $rATE/\Delta$:
{$nATE$ Assumptions:}
For some $m > 2$ , we have that:
}
B.1 imposes some regularity on the tails of the conditional potential outcomes. B.2 defines the $L_2$ and uniform approximation rates using the basis functions for function class $\mathcal{G}$. If the basis is sufficiently rich to span $\mathcal{G}$, we say it is correctly specified and $c_k \rightarrow 0$ as $k \rightarrow \infty$. However, our results allow for the case of misspecification, i.e. $c_k \not\rightarrow 0$. This is a standard characterization in the literature on nonparametric series methods, see e.g. \citeA{Belloni2015SomeResults} for more details and examples.
B.3 is crucial: It requires high-quality approximation capabilities of the first-stage machine learning methods for the nuisance quantities. In the case of a finite-dimensional, bounded basis $\sup_{z\in\mathcal{Z}}||b(z)||_{\infty} < C$ and finite $J$, the conditions can be simplified to $\sqrt{nk}s_{0,n,2}(s_{n,2} + m_{n,2}) = o(1)$. This means that the products of the nuisance quantities for the conditional control propensity and treatment propensities/potential outcome means have to converge at least at rate $o((nk)^{-1/2})$ identical to conditional ATE estimation in \citeA{Semenova2021DebiasedFunctions}. With fixed basis, as in unconditional binary ATE estimation, it reduces to the well-known requirement that nuisance functions attain rate $o(n^{-1/4})$ Chernozhukov2018. For many treatments $J\rightarrow\infty$, flexible $k$, and/or machine learning estimators, the convergence requirements can be more demanding. We discuss these cases and corresponding rate requirements in Section (ref).
B.4 controls the approximation error from linearization taking into account the unknown design matrix $Q$ of increasing dimension. The condition is equivalent to linearization in conventional least squares series estimation Belloni2015SomeResults.\footnote{This suggests that, for specific series methods such as splines Huang2003LocalRegression and local partitioning estimators Cattaneo2020LargeEstimators, the rate can be improved to $\sqrt{\xi_k^2 \log k/n}(1+ \sqrt{\log k}l_kc_k)$, see also \citeA{Belloni2015SomeResults}, Section 4 and \citeA{Cattaneo2020LargeEstimators}, Remark SA-4 of their supplemental appendix.} Note that this rate does not depend on $J$ as, for linearization, the treatment dimension enters only through estimation of the expanding set of nuisance parameters. Once the difference between true and estimated nuisance parameters is controlled for via B.3, there is no difference to the standard series estimation/binary ATE case with no or known nuisances.
B.5 controls the rate of the basis function relative to approximation error such that the Lindeberg condition holds. Note that this rate is required to be faster by a factor of $J$ relative to conventional series estimation. This is due to the fact, that the tails of the summands that determine the first-order asymptotics are selected from a combination of $J$ different potential outcome errors $\varepsilon_i(t)$ for $t=1,\dots,J$. Thus, the conditions for the many treatments case are somewhat stronger then the ones expected for series estimation or (conditional) ATE estimation under a moment assumption such as B.1 and finite $J$.
{$rATE/\Delta$ Assumptions}
For some $m > 2$ , we have that:
}
We discuss and contrast Assumptions C.1--C.6 with B.1--B.5: C.1 and C.2 are equivalent to B.1 and B.2 with potentially different $m$, $c_k$, and $l_k$. C.3 controls for the estimation of nuisance parameters. In the case of a bounded basis, the condition for $B_n$ reduces to $\sqrt{nk}Js_{n,2}m_{n,2} = o(1)$ which is expected to be equivalent to the $nATE$ machine learning bias rate when $J$ is finite. However, for large $J$, estimating potential outcome means at $m_{n,2}$ is generally slower than estimating control propensities at $s_{0,n,2}$. In the parametric case, for example, we have that $m_{n,2} = \sqrt{J/n} = \sqrt{J}s_{0,n,2}$, see Section (ref). An equivalent argument holds for $\Lambda_n$ leading to an additional $\sqrt{J}$ factor compared to the $nATE$ case. Thus, the product rates have to be faster by a factor of $\sqrt{J}$ in this case, i.e. $rATE/\Delta$ generally require somewhat higher quality first-stage learners in comparison to the $nATE$.
C.4 provides the error from the linearization. Note that there is an additional term $R_{n,\pi}$ due to the moment functions not being Neyman-orthogonal with respect to the unconditional weights $\pi$. It puts an additional restriction on the growth of the number of treatments. The first condition reduces to $\sqrt{kJ^3/n} = o(1)$ in case of parametric nuisance quantities. $R_{n,Q}$ corresponds to the $nATE$ case plus an additional term of order $\sqrt{\xi_k\log k/n}k^{1/2}J^{1/4}$. This is a result of the interaction between estimation error from estimating the unconditional weights with design matrix $Q$.\footnote{Again, for specific series such as splines or local partitioning, we conjecture that a faster rate of $\sqrt{\xi_k\log k/n}\sqrt{\log k}J^{1/4}$ is attainable.} The additional factors are only of order $J^{1/2}$ and $J^{1/4}$ compared to the $nATE$. This is due to the superefficiency of the unconditional probability estimator whenever $J\rightarrow\infty$.
C.5 is similar to B.5 but requires more stringent conditions on $J$ and $\xi_k$ compared to the $nATE$. The Lindeberg condition for $nATE/\Delta$ is driven by the tails of a weighted combination of moment functions from many treatment groups which can have high variance when $J$ is large. It is more restrictive compared to the $nATE$, as, for the latter, the weight for each $t$-specific moment function $\psi_i^{[t,0]}(\eta)$ is the actual treatment propensity $e_t(X_i)$, see Table (ref). Thus, inverse propensity scores disappear leading to a lower variance for the $nATE$ explaining the additional $J$-dependent factors between B.5 and C.5.
The first condition in C.6 rules out the degenerate case where $nATE = rATE$. Naturally, if this applies, the weaker Assumptions B.1--B.5 can be used instead. The second condition excludes the hypothetical case where the sum of noise plus approximation error is perfectly negatively correlated with the ($\gamma$-weighted) error from estimating the unconditional weights $\pi$. Both restrictions are expected to always hold in practice and can also be assessed by looking at the empirical analogues of $\Omega$, $\Omega_1$, and $\Omega_2$.
For the estimation of the asymptotic variance, we also assume that A.V holds. The corresponding high-level conditions and discussion can be found in Appendix (ref).
} A.V can require somewhat stronger moment and growth conditions for basis and/or number of treatments. For example, for the $nATE$, they reduce to the same rates required by \citeA{Semenova2021DebiasedFunctions}, Theorem 3.3, condition (ii) with factor $n^{1/m}$ replaced by $(nJ)^{1/m}$. Under finite $J$, they are again equivalent. We obtain the following Theorem:
Theorem (ref) demonstrates the asymptotic validity of the confidence intervals proposed in (ref). The result accommodates the case of misspecification often present in applied research. It is most useful under the additional undersmoothing condition that makes any misspecification bias vanish sufficiently fast.\footnote{ In particular, when $\mathcal{G}$ is in a $s$-dimensional ball on $\mathcal{Z}$ of finite diameter (a H\"older class of smoothness order $s$) then the condition simplifies to $n^{1/2}k^{-(\frac{1}{2} + \frac{s}{d})}\log(k) \rightarrow 0$, see also \citeA{Belloni2015SomeResults}, Comment 4.3 for additional details. Note that such undersmoothing does in general not admit IMSE optimal $k$ choices. Alternatively, bias-correction methods could be employed Cattaneo2020LargeEstimators.}
Theorem (ref) extends to alternative combinations of Neyman-orthogonal scores other than $rATE$ or $\Delta$. In particular, the results for the $rATE$ can directly be applied to any convex combination of conditional average treatment effects as long as the weights are either (i) deterministic sequences (relative to $J$) or (ii) can be estimated at the same rate as $\pi_t$. This can be useful when comparing heterogeneity of a given selection mechanism to alternative, hypothetical (estimated or true) allocation policies different from random selection as considered in this paper even when there are many different treatments.
In this section we provide some basic examples and intuition about the properties of probability and nuisance function estimation when $J$ is large and how this relates to the machine learning bias Assumptions B.3/C.3. We first discuss the necessity of Assumption A.3 and the consequences for the unconditional probability estimates. We then show how these properties translate into different convergence rates for propensity scores and potential outcome means under simplified parametric assumptions. We then discuss the explicit requirements for the machine learning bias Assumptions B.3/C.3 for the flexible high-dimensional nuisance parameter case using Lasso methodology under approximate sparsity and many treatments.
Consider the second part of Assumption A.3: If $J$ is large, then $J\pi_t \lesssim 1$ is a necessary requirement. Because if the product diverges, $J \pi_t > 1$ causes a contradiction with the constraint $\sum_{t\neq 0}\pi_t = 1 -\pi_0$. In principle, one could allow for some $t\neq 0$ such that $J\pi_t = o(1)$. However, restricting A.3 to hold for all $t\neq 0$ is without loss of generality as otherwise it would be asymptotically equivalent to a regime where only a smaller subset of treatments $J' < J$ obey A.3. Thus, all further assumptions and rates would be equivalent with $J$ being replaced by the new number of asymptotically relevant treatments $J'$.
(Local) superefficiency of frequency estimators $\hat{\pi}_t = \frac{1}{n}\sum_{i=1}^nD_{t,i}$ and other nuisance parameters is not a new discovery and has been exploited and discussed in different places in the literature, see e.g. \citeA{Stoye2009MoreParameters}. For general $P(D_{t,i}=1) = \pi_t$, note that $V[D_{t,i}] = \pi_t(1-\pi_t)$. Hence, $ \sqrt{n}(\hat{\pi}_t - \pi_t) \overset{d}{\rightarrow} \mathcal{N}(0,\pi_t(1-\pi_t))$ which implies that $|\hat{\pi}_t - \pi_t| \lesssim_P (n/\pi_t)^{-1/2} \lesssim (nJ)^{-1/2}$ due to Assumption A.3. Thus under the many treatments regime $J \rightarrow \infty$, the frequency estimator is superefficient, i.e. converges at a quicker rate than $n^{-1/2}$.
Superefficiency of $\hat{\pi}_t$ spills over to frequency-based/parametric estimators of propensity scores. For example, consider the case where $\mathcal{X}$ is discrete (finite-dimensional) with $f_x := P(X_i=x) > 0$ for all $x \in \mathcal{X}$. Consider a simple frequency-based estimator for the treatment propensity $t$ as
where the additional indicator in the denominator assures existence. By standard arguments, it follows that, for each $x \in \mathcal{X}$, $\sqrt{n}(\hat{e}_t(x) - e_t(x)) \overset{d}{\rightarrow} \mathcal{N}(0,e_t(x)(1-e_t(x))/f_x)$ which implies that $|\hat{e}_t(x) - e_t(x)| \lesssim_P (n/e_t(x))^{-1/2} \lesssim (nJ)^{-1/2}$ due to Assumption A.3. Thus, the frequency-based finite-dimensional/parametric propensity score has the same superefficiency property as the unconditional frequency estimator.
Parametric estimators of potential outcome means, however, are not superefficient. On the contrary, convergence rates are generally slower under the many treatments regime. For example, consider a frequency-based estimator similar to the one for the propensity score:
where $\mathcal{X}$ is again assumed to be discrete (finite-dimensional) with $f_x := P(X_i=x) > 0$ for all $x \in \mathcal{X}$. Without loss of generality, assume that $E[(Y_i(t) - \mu_t(X_i))^2|X_i=x] = \sigma^2$ (Assumption B.1/C.1 would suffice as well). Again, by standard arguments, we have that, for all $x\in\mathcal{X}$, $\sqrt{n}(\hat{\mu}_t(x) -\mu_t(x)) \overset{d}{\rightarrow} \mathcal{N}(0,\sigma^2/e_t(x)f(x))$ which implies that $|\hat{\mu}_t(x) -\mu_t(x)| \lesssim_P (ne_t(x))^{-1/2} \lesssim (n/J)^{-1/2}$ due to Assumption A.3. Thus, the estimator converges at a slower than parametric rate. In fact, using the rates from Section (ref), yields, for any $t\neq 0$, $s_{n,t,2} = Jm_{t,n,2} = \sqrt{J}s_{0,n,2}$. Thus, corresponding components in the machine learning bias assumptions B.3/C.3 will be of equal rate in the parametric case.
Here we provide some intuition regarding the use of nuisance function estimation using (frequency-based) Lasso in high-dimensional approximately sparse models. In particular, we say potential outcome means are generated by $\mu_t(x) = x'\theta_0 + r_{\mu_t}(x)$ where $x \in \mathcal{X} \subseteq \mathbb{R}^{p_{\mu_t}}$ (and equivalently for $e_t(x)$ with logistic link). $p_{\mu_t}$ denotes the number of available regressors that is allowed to be high-dimensional and grow with $n$. Assume that the typical regularity conditions for Lasso hold as in \citeA{Semenova2021DebiasedFunctions}, Lemma B.1. Denote $s_{e_t}^*$ and $s_{\mu_t}^*$ as the corresponding sparsity indices that obey these assumptions. For simplicity, let the number of available regressors and sparsity indices coincide for propensity score and potential outcome estimation, i.e. $s_{e_t}^* = s_{\mu_t}^* \equiv s^*$ and $p_{e_t} = p_{\mu_t} \equiv p$. For the machine learning bias for the $nATE$ in Assumption B.3 we then conjecture that
based on the same argument as for the parametric frequency-based estimation above.
Thus $B_n^{[nATE]} = o(1)$ requires that sparsity indices have to obey
For $rATE/\Delta$, it is similarly required that
Thus, sparsity conditions for $rATE/\Delta$ are stronger than for $nATE$ in the many treatments regime. This reflects the different variability due to different weighting between the decomposition parameters as the $nATE$ weights minimize variance. Comparing rate requirements for $nATE$ to the ones \citeA{Semenova2021DebiasedFunctions}, we find that $s^*$ here must be slower by a factor of $\sqrt{J}$ compared to their Lemma B.1. This is the price for the expanding set of nuisance parameters when estimating treatment propensities and potential outcome means for each treatment level separately instead of imposing the binary treatment structure to begin with. Moreover, the nonparametric heterogeneity step adds an additional $\sqrt{k}$ to the sparsity requirements compared to standard double machine learning estimation of the binary ATE in \citeA{Chernozhukov2018}. An analogous derivation can be conducted for $\Lambda_n$ as well. Note that the given sparsity assumption here is for each treatment probability separately. In practice, we might want to impose some (group-based) sparsity across treatments to improve estimation when many treatments are available. In this case, rates can be improved depending on the total complexity of the propensity scores Farrell2015. We leave extensions along these lines for future work.
In this section we analyze the finite sample performance of the analytical confidence bounds proposed in Section (ref). In particular, we evaluate the empirical coverage rates of the corresponding confidence intervals in setups with heterogeneous effective treatment probabilities for all the decomposition parameters. We consider the case of three effective treatment levels and a best linear predictor for the heterogeneity analysis using different sample sizes and total number of confounding variables including high-dimensional designs. In the heterogeneity step, we regress the estimated pseudo outcomes on a single confounder and evaluate the coverage rates for the parameters of the linear predictor. All nuisance parameters are estimated via 2-fold cross-fitting using $\ell_1$-regularized linear regression for the outcome models as well as $\ell_1$-regularized multinomial logistic regression for the propensity scores.\footnote{We have also conducted similar simulations for correctly specified parametric models. Coverage rates are similar or slightly better in small $n$/large $p$ setups. Results are available upon request.} Tuning parameter selection is done via 5-fold cross-validation. The true models satisfy the necessary sparsity assumptions required for high-quality approximation of the machine learning methods Belloni2013LeastModels,Farrell2015,Belloni2016Post-selectionControls. We consider two parameterizations: Design A has linear log-odds and potential outcomes in the heterogeneity dimension while Design B also includes nonlinear components (sign, trigonometric, polynomial, and rectified linear functions). For more details on the designs please consider Appendix (ref).
Table (ref) contains the coverage rates of the confidence intervals based on (ref) using double machine learning at a significance level of $5\%$ for both designs with sample sizes $n=1000, 5000$ and number of confounders $p=10,100,1000,5000$. For $rATE(z)$ and $nATE(z)$ results are always very close to the nominal coverage rate in both designs. For $\Delta(z)$, there is some undercoverage for the intercept $\alpha$ in design A which increases in the number of parameters and decreases with the sample size as expected. The slope parameter coverage for $\beta$ for $\Delta(z)$ is very close to nominal for any sample size, confounding dimension, or design. Overall the inference based on the asymptotic approximation in (ref) seems to be mostly reliable in finite samples.
The detrimental effect of smoking on birth weight and its economic costs are well documented \cite<see e.g.>[and references therein]{Almond2005TheWeight,Abrevaya2006EstimatingApproach}. Beyond the standard average effects it is important to understand the heterogeneous effects to e.g. identify for which subgroups interventions to reduce smoking during pregnancy would be most beneficial. \citeA{Abrevaya2006EstimatingApproach} documents that the negative effect of smoking is less pronounced for black compared to white mothers in a standard subgroup analysis. A variety of papers analyze heterogeneous effects of smoking as a function of mother's age Abrevaya2015EstimatingEffects,Lee2017,Zimmert2019NonparametricConfounding,Fan2022EstimationData. They all document increasingly negative effects with higher age. The aforementioned studies consider “smoking yes/no” as the binary treatment. \citeA{Cattaneo2010EfficientIgnorability} notes that smoking is not a homogeneous treatment, but that the negative effects become more extreme for higher intensities of smoking. Thus, the binary indicator “smoking” represents only an aggregation of smoking intensities which directly affect birth weight. This corresponds to Scenario 1. We investigate whether the heterogeneous effects documented in the literature can be at least partly explained by different smoking intensities of different groups.
We analyze the dataset of \citeA{Almond2005TheWeight} used by \citeA{Cattaneo2010EfficientIgnorability} with five intensities of smoked cigarettes per day as the effective treatment $T_i \in \mathcal{T} = \{0,1-5,6-10,11-15,16-20,>20\}$, the binary indicator defined as $D_i = \mathbbm{1}(T_i > 0)$, the outcome $Y_i$ being birth weight in gram, and the confounders $X_i$ including age, education, ethnicity, and marital status of mother and father as well as health indicators and pregnancy history of the mother.\footnote{We thank Matias Cattaneo for sharing the full data. A random subsample is available on his \href{https://github.com/mdcattaneo/replication-C_2010_JOE}{GitHub repository}.} The dataset comprises 511,940 observations after removing the 0.1% of the observations with missing values in relevant variables and 52 confounders. The nuisance parameters are estimated with 2-fold cross-fitting using an ensemble learner of the unconditional mean, random forests, lasso and ridge regression with 2-fold cross-validated weights. For the propensity scores, we use logistic lasso and ridge.
Smoking behavior differs along the heterogeneity variables ethnicity and age showing that white and older smoking mothers smoke more heavily.\footnote{Appendix (ref) and in particular Figure (ref) provides the smoking distributions by heterogeneity variables.} Combined with the result of \citeA{Cattaneo2010EfficientIgnorability} that different smoking intensities have different effects, this suggests that at least part of the heterogeneity could be explained by different smoking intensities.
Figure (ref) contains the result of the decomposition for the heterogeneity variable “ethnicity”. The upper panel shows the decomposition for each subgroup. It is obtained by running an OLS regression of the estimand specific pseudo-outcome on a set of four dummy variables indicating ethnicity of the mother without a constant. The $nATE$ in the left part corresponds to standard subgroup analysis. Like previous studies, we find that smoking reduces the birth weight of newborns more for white women than for Blacks, Hispanics and others. Given that smoking is a binarized treatment, it is not clear how much is really effect heterogeneity and how much is driven by the fact that subgroups differ in their smoking intensity. The decomposition term $rATE$ fixes the intensity of smoking for all subgroups at the population level. It provides the subgroup specific effect of smoking if all groups had the same smoking intensity. Under this harmonized smoking intensity the negative average effect of smoking is smaller for white women and larger for the others. $\Delta$ in the right graph quantifies the difference between $nATE$ and $rATE$. It shows relatively small differences suggesting that different smoking intensities are not the main driver of the differences between white mothers and the other groups. However, they are also not negligible as the lower panel of Figure (ref) shows. It quantifies the heterogeneous effects by subtracting the effects for white mothers from the other three groups. We observe that a significant portion of the difference between black/hispanic mothers and white mothers is driven by different smoking intensities. For black vs. white mothers the difference in the $nATE$ is 69 gram of which 12% are due to different smoking intensities ($\Delta = 8$). For hispanic vs. white mothers it explains around 17% ($\Delta = 14$).
Figure (ref) depicts the heterogeneity analysis along age. We use B-splines as basis functions of age. We select the nodes and order via leave-one-out cross-validation for each parameter and apply the most flexible/low-bias model for all parameters to ensure that the $rATE$ and $\Delta$ curves add up to the $nATE$ curve. The left panel of Figure (ref) replicates the well-established findings of previous papers that the $nATE$ is much smaller for younger mothers than for older mothers. In the extreme case where different smoking intensities would fully explain the heterogeneous $nATE$, we would see a flat $rATE$ curve in the middle graph. However, we only observe that the effect of teenage mothers would be more negative if we harmonize smoking intensity over all age groups. Overall, only a relatively small part of the heterogeneous effects of the binarized smoking indicator can be attributed to different smoking intensities and the larger part seems to be driven by different age groups actually being affected differently.
We illustrate Scenario 2 with an evaluation of the Job Corps (JC) program. JC operates since 1964 and is the largest training program for disadvantaged youth aged 16-24 in the US \cite<see>[for a detailed description]{Schochet2001NationalOutcomes,Schochet2008DoesStudy}. The roughly 50,000 participants per year receive an intensive treatment as a combination of different components like academic education, vocational training, and job placement assistance. Participants plan their educational and vocational curricula together with counselors. This means that although the variable “access to JC” is a binary indicator, different versions of JC participation are conceivable. Heterogeneous effects might thus be driven by different effectiveness of JC for different groups, by different tailoring of the curriculum, or a combination thereof.
We investigate this based on data from an experiment in 1994-1996 Schochet2019ReplicationStudy.\footnote{The data is available as public use file via \href{https://doi.org/10.3886/E113269V1}{https://doi.org/10.3886/E113269V1}.} This experiment is basis of a variety of studies looking at different aspects of JC. Many of them report gender differences in the effectiveness of the programs with women benefiting less than men from access to JC \cite<e.g.>{Schochet2001NationalOutcomes,Schochet2008DoesStudy,Flores2012EstimatingCorps,Eren2014WhoProgram,Strittmatter2019HeterogeneousApproach}. One potential explanation for this finding is that men and women focus on average on different vocational training within JC. In particular men receive more often training for higher paying craft jobs, while women focus more often on training for the service sector Quadagno1995TheCorps,Inanc2017GenderGap.\footnote{Appendix (ref) and in particular Figure (ref) provides the distribution of trainings by gender.} We apply our decomposition method to investigate this potential explanation of the gender gap in program effectiveness.
We analyze the intention to treat effect (ITT) of the binary variable indicating random access to JC ($D_i$) on weekly earnings four years after random assignment ($Y_i$). We consider 11 versions of the effective treatment ($T_i$): (i) No JC if eligible individuals did not participate (non-compliers), (ii) JC without vocational training if eligible individuals entered JC but did not receive vocational training, (iii-ix) training for jobs in the clerical, health, auto mechanics, welding, electrical/electronics, construction, or food sector, (x) other vocational training, (xi) training for multiple sectors.
Nuisance parameters are estimated with the same ensemble as in Section (ref) using 5-fold cross-fitting. We control for 55 covariates that include pre-treatment information about labor market history, socio-economic characteristics, education, health, crime, and JC related variables. These control variables overlap mostly with those of \citeA{Flores2012EstimatingCorps} who also employ an unconfoundedness strategy.\footnote{Considering second-order interactions results in a total of 1428 variables after screening for nearly empty cells (less than 1% observations) and nearly perfectly correlated variables (correlation higher than 0.99).} In total we work with a sample of 9,708 observations.
The unconditional $nATE$, corresponding to the ITT of eligibility for JC on monthly earnings, is estimated at \$14.2 (S.E. 3.8), which is an increase of 7% in line with previous studies. The unconditional $rATE$ is larger (\$17.4, S.E. 4.1) suggesting that hypothetical random allocation of the curricula would yield higher average outcomes compared to the actual assignment. However, the unconditional difference $\Delta$ is insignificant ($\$-3.1$, S.E. 1.8). This suggests that, on average, the selection of versions is not statistically distinguishable from random allocation.
Figure (ref) depicts the decomposition of the gender specific effects. We observe that the effect for women with the actual composition of vocational training ($nATE$) is not significant at $\alpha=0.05$, but under the hypothetical treatment composition of the population would show a clear positive effect ($rATE = \$16.0$). The gender gap in effectiveness basically disappears when both groups receive the same hypothetical mix of vocational training. The right part of Figure (ref) suggests that 73% of the gender gap in the effectiveness of JC is due to different training curricula. This means that the worse than average performance of the assignment mechanism seen in the unconditional parameters is mostly driven by women. While the assignment to vocational training for men is as well targeted as random assignment, for women it is even worse. This indicates that there is room for improvement to target vocational training in general and for women in particular. Our results suggest that removing the worse than random targeting of vocational training for women could decrease the gender gap in the effectiveness of access to JC.
The method proposed in this paper provides a practical way of decomposing effect heterogeneity obtained from analyzing a binary treatment indicator that does not coincide with the effective multi-valued treatment. The approach likely extends to other causal parameters and identification strategies such as continuous effective treatments, selection on unobservables/instrumental variables, or mediation analysis. It would also be interesting to see whether the ideas could be further developed to find the most relevant dimensions of effective treatments for cases with multiple treatment versions instead of requiring the researcher to manually specify them.
The conceptual and empirical results highlight that potential treatment heterogeneity underlying the analyzed binary indicator should be taken more seriously and explicitly discussed in applications, especially when interpreting heterogeneous effects. The decomposition provides one principled way to do this. It requires to observe the effective treatment. Thus, data collection can anticipate the goal of better understanding treatment heterogeneity by recording effective treatment information beyond a binary indicator. Furthermore, the decomposition shows that reducing the analysis to such binary indicators, while facilitating the analysis, comes at the cost of a more intricate interpretation of empirical results.
{\setstretch{1.45} }