Contamination Bias in Linear Regressions

Contamination Bias in Linear Regressions

\pagenumbering{Alph}

\pagenumbering{arabic} \setstretch{1.25}

Introduction

Consider a linear regression of an outcome $Y_i$ on a vector of treatments $X_i$ and a vector of flexible controls $W_i$. The treatments are assumed to be as good as randomly assigned conditional on the controls. For example, $X_i$ may indicate the assignment of individuals $i$ to different interventions in a stratified \ac{RCT}, with the randomization protocol varying across some experimental strata indicators in $W_i$. Or, in an education \ac{VAM}, $X_i$ might indicate the matching of students $i$ to different teachers or schools with $W_i$ including measures of student demographics and lagged achievement which yield a credible selection-on-observables assumption. The regression might also be the first stage of an \ac{IV} regression leveraging the assignment of multiple decision-makers (e.g.\ bail judges) indicated in $X_i$, which is as-good-as-random conditional on some controls $W_i$. These sorts of regressions are widely used across many fields in economics.\footnote{Prominent \acp{RCT} where randomization probabilities vary across strata include Project STAR krueger1999experimental and the RAND Health Insurance Experiment manning1987health. Prominent \ac{VAM} examples include studies of teachers NBERw14607,chetty2014measuring, schools angrist2017leveraging,ahpw20,mountjoy2020returns, and healthcare institutions abaluck2020mortality,geruso2020all. Prominent “judge \ac{IV}” examples include kling2006incarceration,maestas2013does,dobbie2015debt.}

This paper shows that such multiple-treatment regressions generally fail to estimate convex weighted averages of heterogeneous causal effects, and discusses solutions to this problem. The problem may be surprising given an influential result in angrist98, showing that regressions on a single binary treatment $D_i$ and flexible controls $W_i$ estimate a convex average of treatment effects whenever $D_i$ is conditionally as good as randomly assigned. We show that this result does not generalize to multiple treatments: regression estimates of each treatment's effect are generally contaminated by a non-convex average of the effects of other treatments. Thus, the regression coefficient for a given treatment arm incorporates the effects of all arms.

We first derive a general characterization of such contamination bias in multiple-treatment regressions.\footnote{Our use of the term “contamination” follows sun2021estimating, and differs from its use in some analyses of clinical trials keoghbrown07 to describe settings where members of one treatment group receive the treatment of another group---what economists typically call “non-compliance”. Our “bias” terminology refers to an implication of our result: if a given treatment has constant effects, but the other treatment effects are heterogeneous, the regression estimand is inconsistent for the given treatment effect.} We show the core problem by focusing on the special case of a set of mutually exclusive treatment indicators, though our characterization applies even when the treatments are not restricted to be binary or mutually exclusive. To separate the problem from the typical challenge of \ac{OVB}, we assume a best-case scenario where the covariate parametrization is flexible enough to include the treatment propensity scores (e.g., with a linear covariate adjustment, we assume that the propensity scores are linear in the covariates). This condition holds trivially if the only covariates are strata indicators. Under these conditions, we show that the regression coefficient on each treatment identifies a convex weighted average of its causal effects plus a contamination bias term given by a linear combination of the causal effects of other treatments, with weights that sum to zero. Thus, each treatment effect estimate will generally incorporate the effects of other treatments, unless the effects are uncorrelated with the contamination weights. Since these weights sum to zero some are necessarily negative---further complicating the interpretation of the coefficients.

Contamination bias arises because regression adjustment for the confounders in $W_i$ is generally insufficient for making the other treatments ignorable when estimating a given treatment's effect, even when this adjustment is flexible enough to avoid \ac{OVB}. To see this intuition clearly, suppose the only controls are strata indicators. \ac{OVB} is avoided when the treatments are as good as randomly assigned within strata. But because the treatments enter the regression linearly, the angrist98 result implies that the causal interpretation of a given treatment's coefficient is only guaranteed when its assignment depends linearly on both the strata indicators and the other treatment indicators. With mutually exclusive treatments, this condition fails because the dependence is inherently nonlinear---the probability of assignment to a given treatment is zero if an individual is assigned to one of the other treatments, regardless of their stratum, but strata indicators affect the treatment probability otherwise. Such dependence generates contamination bias.

Contamination bias also arises under an alternative “model-based” identifying assumption that---rather than making assumptions on the treatment's “design” (i.e.\ propensity scores)---posits that the covariate specification spans the conditional mean of the potential outcome under no treatment, $Y_i(0)$. In a linear model with unit and time fixed effects, this reduces to the parallel trends restriction often used in \ac{DiD} and event study regressions. It is common for $X_i$ to include multiple indicators in such settings---for example, the leads and lags relative to a treatment adoption date used to support the parallel trends assumption or estimate treatment effect dynamics.\footnote{Alternatively $X_i$ may indicate multiple contemporaneous treatments, as in certain “mover” regressions.} We show that replacing the restriction on propensity scores in our characterization with an assumption on $Y_i(0)$ generates an additional issue: the own-treatment weights are negative whenever the implicit propensity score model used by the regression to partial out the covariates and the other treatments fits probabilities greater than one. This result shows that the negative weighting and contamination bias issues documented previously in the context of two-way fixed effects regressions goodman2021difference, sun2021estimating, de2020aer, de2020two, callaway2021difference, borusyak2021revisiting,wooldridge2021mundlak,hull2018movers are more general---and conceptually distinct---problems.\footnote{Our analysis also relates to issues with interpreting multiple-treatment \ac{IV} estimates behaghel2013robustness, kirkeboen2016field, kline2016headstart,hull2018isolateing,leeMultivalued,bhuller_sigstad.} Negative weighting arises because regressions leveraging model-based restrictions on $Y_i(0)$ may fit treatment probabilities exceeding one. Contamination bias arises because additive covariate adjustments don't account for the non-linear dependence of a given treatment on the other treatments and covariates. This generates a different form of propensity score misspecification: a non-zero fitted probability of a given treatment, even when one of the other treatments is known to be non-zero.\footnote{While our results are framed in the context of a causal model, we show how analogous results apply to descriptive regressions which seek to estimate averages of conditional group contrasts without assuming a causal framework---as in studies of outcome disparities across multiple racial or ethnic groups, studies of regional variation in healthcare utilization or outcomes, or studies of industry wage gaps.}

We then discuss three solutions to the contamination bias problem, and their trade-offs. These solutions apply when the propensity scores are non-degenerate, such as in an \ac{RCT} or other “design-based” regression specification.\footnote{Solving the contamination bias problem under model-based identification approaches requires either targeting subpopulations of the treated or applying substantive restrictions on the conditional means of potential outcomes under treatment. We do not explore this case as it has already been studied extensively in the \ac{DiD} context de2020two, sun2021estimating, callaway2021difference,borusyak2021revisiting,wooldridge2021mundlak.} First, a conceptually principled solution is to adapt approaches to estimating the \ac{ATE} of a conditionally ignorable binary treatment to the multiple treatment case cattaneo10, chernozhukov2018double,ChNeSi21,ReZu20, GrPi22. For example, one could run a regression that includes interactions between the treatments and demeaned controls, or combine such regression with inverse propensity score weighting for doubly-robust estimation. Such \ac{ATE} estimators work well under strong overlap of the covariate distribution for units in each treatment arm. But they may be imprecise under limited overlap or be outright infeasible with overlap failures---common scenarios in observational studies crump2009dealing.

This practical consideration motivates an alternative approach: estimating a weighted average of treatment effects, as regression does in the binary treatment case, while avoiding the contamination bias with multiple-treatments. We derive the weights that are easiest to estimate, in the sense of minimizing a semiparametric efficiency bound under homoskedasticity. This \ac{EW} scheme is always convex; it corresponds to weighting schemes previously proposed in crump2006dealing, liBalancing, and lili19. The weights also coincide with the implicit linear regression weights when the treatment is binary (i.e.\ the angrist98 case). In the multiple treatment case, the \ac{EW} scheme that allows the weights to be treatment specific can be implemented by a simple second solution: a linear regression which restricts estimation to the individuals who are either in the control group or the treatment group of interest. Since the weights are treatment-specific, these one-treatment-at-a-time regressions preclude direct comparisons across treatment arms. The third solution is to impose common weights across treatments in the \ac{EW} scheme; these weights can be implemented using a weighted regression approach. We show how researchers can gauge the extent of contamination bias in practice and implement these tools with a new R and Stata package, multe.\footnote{The package is available at CRAN (R) and \url{https://github.com/gphk-metrics/stata-multe} (Stata).}

We study the empirical relevance of contamination bias in nine applications: six \acp{RCT} with stratified randomization and three observational studies of racial disparities. We find economically and statistically significant bias in two of the three observational studies with no evidence for bias in any of the experimental studies. In a detailed analysis of one experiment---the Project STAR trial---we show that the lack of contamination bias is driven by small variation in the contamination weights, rather than limited effect heterogeneity. This analysis highlights the importance of conducting contamination bias diagnostics, particularly in observational studies where covariates are expected to generate high variability in propensity scores, and thus likely in contamination weights.

We structure the rest of the paper as follows. (ref) illustrates contamination bias in a simple stylized setting. (ref) characterizes the general problem, and discusses connections to previous analyses. (ref) provides three solutions, and gives guidance for measuring and avoiding contamination bias in practice. (ref) illustrates these tools in nine applications. (ref) concludes. (ref) collects all proofs and extensions. (ref) discusses the connection between our contamination bias characterization and that in the \ac{DiD} literature. Details on the applications and additional exhibits are given in (ref).

Motivating Example

We build intuition for the contamination bias problem in two simple examples. We first review how regressions on a single randomized binary treatment and binary controls identify a convex average of heterogeneous treatment effects. We then show how this result fails to generalize when we introduce an additional treatment arm. We base these examples on a stylized version of the Project STAR experiment, which we return to as an application in (ref). The simple structure of these examples helps isolate the core mechanisms of contamination bias. Later sections consider non-experimental settings with richer control specifications, both theoretically and empirically.

Convex Weights with One Randomized Treatment

Consider the regression of an outcome $Y_i$ on a single treatment indicator $D_{i}\in\{0,1\}$, a single binary control $W_i\in\{0,1\}$, and an intercept:

By definition, $U_i$ is a mean-zero regression residual that is uncorrelated with $D_i$ and $W_i$. For example, analysing the Project STAR trial, krueger1999experimental primarily studied the effect of small class size $D_i$ on the test scores $Y_i$ of kindergartners indexed by $i$. Project STAR randomized students to classes within schools, with the fraction of students assigned to small classes varying by school due to the varying number of total students in each school. To account for this, krueger1999experimental included school fixed effects as controls. Such specifications are often found in stratified \acp{RCT} with varying treatment assignment rates across a set of pre-treatment strata. If we imagine two such strata, demarcated by a binary indicator $W_{i}$, then (ref) corresponds to a stylized two-school version of a Project STAR regression.

We wish to interpret the coefficient $\beta$ in terms of the causal effects of $D_i$ on $Y_i$. For this we use potential outcome notation, letting $Y_i(d)$ denote the test score of student $i$ when $D_i=d$. Individual $i$'s treatment effect is then given by $\tau_{1i}=Y_i(1)-Y_i(0)$, and we can write realized achievement as $Y_i=Y_i(0)+\tau_{1i}D_i$. Since treatment assignment is random within schools, $D_i$ is conditionally independent of potential outcomes given $W_i$: $\left(Y_i(0), Y_i(1)\right)\perp D_i\mid W_i$.

angrist98 showed that regression coefficients like $\beta$ identify a convexly-weighted average of within-strata \acp{ATE}. In our Project STAR example, this result shows that:

gives a convex weighting scheme, and $\tau_{1}(w)=E[Y_i(1)-Y_i(0)\mid W_i=w]$ is the \ac{ATE} in school $w\in\{0,1\}$. Thus, in our example the coefficient $\beta$ identifies a weighted average of school-specific small classroom effects $\tau_{1}(w)$ across the two schools.

(ref) can be derived by applying the \ac{FWL} Theorem. The multivariate regression coefficient $\beta$ can be written as a univariate regression coefficient from regressing $Y_i$ onto the population residual $\tilde{D}_i$ from regressing $D_i$ onto $W_i$ and a constant:

where we substitute the potential outcome model for $Y_i$ in the second equality. Since $W_i$ is binary, the propensity score $E[D_i\mid W_i]$ is linear and the residual $\tilde{D_i}$ is mean independent of $W_i$ (not just uncorrelated with it): $E[\tilde{D_i}\mid W_i]=0$. Therefore,

The first equality in (ref) follows from the law of iterated expectations, the second equality follows by the conditional random assignment of $D_i$ and the third equality uses $E[\tilde{D_i}\mid W_i]=0$. Hence, the first summand in (ref) is zero. Analogous arguments show that

where $\operatorname{var}(D_i\mid W_i)=E[\tilde{D}_i^2\mid W_i]$ gives the conditional variance of the small-class treatment within schools. Since $E[\operatorname{var}(D_i\mid W_i)]=E[E[\tilde{D}_i^2\mid W_i]]=E[\tilde{D}_i^2]$, it follows that we can write the second summand in (ref) as

proving the representation of $\beta$ in (ref).

The key fact underlying this derivation is that the residual $\tilde{D}_i$ from the auxiliary regression of the treatment $D_i$ on the other regressors $W_i$ is mean-independent of $W_{i}$. By the \ac{FWL} theorem, treatment coefficients like $\beta$ can always be represented as in (ref) even without this property. We next show, however, that the remaining steps in the derivation of (ref) fail when an additional treatment arm is included. This failure can be attributed to the fact that the auxiliary \ac{FWL} regression delivers a treatment residual that is uncorrelated with---but not mean-independent of---the other regressors. The lack of mean independence leads to an additional term in the expression for the regression coefficient.

Contamination Bias with Two Randomized Treatments

In reality, Project STAR randomized students to three mutually exclusive conditions within schools: a control group with a regular class ($D_i=0$), a treatment that reduced class size ($D_i=1$), and a treatment that introduced full-time teaching aides ($D_i=2$). We incorporate this extension of our stylized example by considering a regression of student achievement $Y_i$ on a vector of two treatment indicators, $X_{i} = (X_{i1}, X_{i2})'$, where $X_{ik}=\1{D_i=k}$ indicates assignment to treatment $k=1,2$. We continue to include a constant and the school indicator $W_i$ as controls, yielding the regression

The observed outcome is now given by $Y_i=Y_i(0)+\tau_{i1}X_{i1}+\tau_{i2}X_{i2}$, with $\tau_{i1}=Y_i(1)-Y_i(0)$ and $\tau_{i2}=Y_i(2)-Y_i(0)$ denoting the potentially heterogeneous effects of a class size reduction and introduction of a teaching aide, respectively. As before, we analyze this regression by assuming ${X}_i$ is conditionally independent of the potential achievement outcomes $Y_i(d)$ given the school indicator $W_i$: $\left(Y_i(0), Y_i(1), Y_i(2)\right)\perp {X}_i \mid W_i$.

To analyze the coefficient on $X_{i1}$, we again use the \ac{FWL} theorem to write

where $\accentset{\approx}{X}_{i1}$ again denotes a population residual, but now from regressing $X_{i1}$ on $W_i$, a constant, and $X_{i2}$. Unlike before, this residual is uncorrelated with but not mean-independent of the remaining regressors $(W_i, X_{i2})$ because the dependence between $X_{i1}$ and $X_{i2}$ is non-linear. When $X_{i2}=1$, $X_{i1}$ must be zero regardless of the value of $W_i$ (because they are mutually exclusive) while if $X_{i2}=0$ the mean of $X_{i1}$ does depend on $W_i$ unless the treatment assignment is completely random. Thus, in general, $\accentset{\approx}{X}_{i1}\neq X_{i1}-E[X_{i1}\mid W_i, X_{i2}]$.

Because $\accentset{\approx}{X}_{i1}$ does not coincide with a conditionally de-meaned $X_{i1}$, we can not generally reduce (ref) to an expression involving only the effects of the first treatment arm, $\tau_{i1}$. It turns out that we nevertheless still have $E[\accentset{\approx}{X}_{i1}Y_{i}(0)]=0$, as in (ref), since the auxilliary regression residuals are still uncorrelated with any individual characteristic like $Y_i(0)$.\footnote{To see this, note that in the auxiliary regression $X_{i1} = \mu_{0} + \mu_{1}X_{i2} + \mu_{2}W_{i} + \accentset{\approx}{X}_{i1}$ we can partial out $W_i$ and the constant from both sides to write $\tilde{X}_{i1} = \mu_{1}\tilde{X}_{i2} + \accentset{\approx}{X}_{i1}$. Thus, $\accentset{\approx}{X}_{i1} = \tilde{X}_{i1} - \mu_{1}\tilde{X}_{i2}$ is a linear combination of residuals which, per (ref), are both uncorrelated with $Y_i(0)$. It follows that $E[\accentset{\approx}{X}_{i1}Y_i(0)]=0$.} The regression thus does not suffer from \ac{OVB}. However, we do not generally have $E[\accentset{\approx}{X}_{i1}X_{i2}\tau_{i2}]=0$. Instead, simplifying (ref) by the same steps as before leads to the expression

as a generalization of (ref). Here $\lambda_{11}(W_i)=E[\accentset{\approx}{X}_{i1}X_{i1}\mid W_i]/E[\accentset{\approx}{X}_{i1}^2]$ can be shown to be non-negative and to average to one, similar to the $\phi$ weight in (ref). Thus, if not for the second term in (ref), $\beta_1$ would similarly identify a convex average of the conditional \acp{ATE} $\tau_{1}(W_i)=E[Y_i(1)-Y_i(0)\mid W_i]$. But precisely because $\accentset{\approx}{X}_{i1}\neq X_{i1}-E[X_{i1}\mid W_i, X_{i2}]$, this second term is generally present: $\lambda_{12}(W_i)=E[\accentset{\approx}{X}_{i1}X_{i2}\mid W_i]/E[\accentset{\approx}{X}_{i1}^2]$ is generally non-zero, complicating the interpretation of $\beta_1$ by including the conditional effects of the other treatment $\tau_{2}(W_i)=E[Y_i(2)-Y_i(0)\mid W_i]$.

The second contamination bias term in (ref) arises because the residualized small class treatment $\accentset{\approx}{X}_{i1}$ is not conditionally independent of the second full-time aide treatment $X_{i2}$ within schools, despite being uncorrelated with $X_{i2}$ by construction. This can be seen by viewing $\accentset{\approx}{X}_{i1}$ as the result of an equivalent two-step residualization. First, both $X_{i1}$ and $X_{i2}$ are de-meaned within schools: $\tilde{X}_{i1} = X_{i1} - E[X_{i1}\mid W_{i}] = X_{i1} - p_{1}(W_{i})$ and $\tilde{X}_{i2} = X_{i2} - E[X_{i2}\mid W_{i}] = X_{i2} - p_{2}(W_{i})$ where $p_j(W_i)=E[X_{ij}\mid W_i]$ gives the propensity score for treatment $j$. Second, a bivariate regression of $\tilde{X}_{i1}$ on $\tilde{X}_{i2}$ is used to generate the residuals $\accentset{\approx}{X}_{i1}$. When the propensity scores vary across the schools (i.e.\ $p_j(0)\neq p_j(1)$), the relationship between these residuals varies by school, and the line of best fit between $\tilde{X}_{i1}$ and $\tilde{X}_{i2}$ averages across this relationship. As a result, the line of best fit does not isolate the conditional (i.e.\ within-school) variation in $X_{i1}$: the remaining variation in $\accentset{\approx}{X}_{i1}$ will tend to predict $X_{i2}$ within schools, making the contamination weight $\lambda_{12}(W_i)$ non-zero.

Illustration and Intuition

A simple numerical example helps make the contamination bias problem concrete. Suppose in the previous setting that school $0$ (indicated by $W_i=0$) assigned only 5 percent of the students to the small classroom treatment, with 45 percent of the students assigned to the full-time aide treatment and the rest assigned to the control group. In school $1$ (indicated by $W_i=1$), there was a substantially larger push for students to be placed into treatment groups with 45 percent of students assigned to a small classroom, 45 percent assigned to a classroom with a full-time aide, and only 10 percent assigned to the control group. Therefore, $p_1(0) = 0.05$ and $p_2(0) = 0.45$ while $p_1(1) = p_2(1) = 0.45$. Suppose that the schools have the same number of students, so that $\Pr(W_{i} = 1) = 0.5$. It then follows from the above formulas that $\lambda_{12}(0) = 99/106$ and $\lambda_{12}(1) =-99/106$.

As reasoned above, the contamination weights are non-zero here because the within-school correlation between the residualized treatments, $\tilde{X}_{i1}$ and $\tilde{X}_{i12}$, is heterogeneous: in school $0$ it is about $-0.2$, so that the value of the demeaned class aide treatment is only weakly predictive of the small classroom treatment, while in school $1$ it is highly predictive with correlation $-0.8$. (ref) in (ref) illustrates this graphically, showing that because the overall regression of $\tilde{X}_{i1}$ on $\tilde{X}_{i2}$ averages over these two correlations, the regression residuals are predictive of the value of the class aide treatment.

To illustrate the potential magnitude of bias in this example, suppose that classroom reductions have no effect on student achievement (so $\tau_1(0)=\tau_1(1)=0$), but that the effect of a teaching aide varies across schools. In school $1$ the aide is highly effective, $\tau_2(1)=1$, (which may be the reason for the higher push in this school to place students into treatment groups) but in school $0$, the aide has no effect, $\tau_2(0)=0$. By (ref), the regression coefficient on the first treatment identifies

Thus, in this example, a researcher would conclude that small classrooms have a sizable negative effect on student achievement---equal in magnitude to around half of the true teaching aide effect in school $1$---despite the true small-classroom effect being zero for all students. This treatment effect coefficient can be engineered to match an arbitrary magnitude and sign by varying the heterogeneity of the teaching aide effects across schools.

To build further intuition for (ref), it is useful to consider two cases where the contamination bias term is zero. First, note that since regression residuals are by construction uncorrelated with the included regressors, $E[\lambda_{12}(W_i)]=E[\accentset{\approx}{X}_{i1}X_{i2}]/E[\accentset{\approx}{X}_{i1}^2]=0$. Therefore, $E[\lambda_{12}(W_i)\tau_2(W_i)]=E[\lambda_{12}(W_i)\tau_2(W_i)]-E[\lambda_{12}(W_i)]E[\tau_2(W_i)]=\operatorname{cov}(\lambda_{12}(W_i), \tau_2(W_i))$. If the average effects of the teaching aide treatment are constant across the two schools, $\tau_2(1)=\tau_2(0)$, then $\tau_2(W_i)$ is constant, and this covariance is zero such that contamination bias disappears. More generally, when the average teaching aide treatment effects across schools $\tau_2(W_i)$ exhibit idiosyncratic variation, in the sense that they have a weak covariance with the contamination weights across schools, the contamination bias term will be small.

Second, consider the case where $X_{i1}$ and $X_{i2}$ are independent conditional on $W_i$---such as when the small classroom and teacher aid interventions are independently assigned within schools, in contrast to the previously assumed mutual exclusivity of these treatments. In this case the conditional expectation $E[X_{i1}\mid W_i, X_{i2}]=E[X_{i1}\mid W_i]$ will be linear, since $X_{i1}$ and $X_{i2}$ are unrelated given $W_i$, and will thus be identified by the auxiliary regression of $X_{i1}$ on $W_i$, $X_{i2}$, and a constant. Consequently, the $\accentset{\approx}{X}_{i1}$ residuals will coincide with $X_{i1}-E[X_{i1}\mid W_i]$. The coefficient on $X_{i1}$ in (ref) can therefore be shown to be equivalent to the previous (ref), identifying the same convex average of $\tau_1(w)$. This case highlights that dependence across treatments is necessary for the contamination bias to arise.

General Problem

We now derive a general characterization of the contamination bias problem, in regressions of an outcome $Y_i$ on a $K$-dimensional treatment vector $X_i$ and flexible transformations of a control vector $W_i$. We focus on the case of mutually exclusive indicators $X_{ik}=\1{D_i=k}$ for values of an underlying treatment $D_i\in\{0,\dotsc, K\}$ (with the $\1{D_i=0}$ indicator omitted). We extend the characterization to a general (i.e.\ potentially non-binary) $X_i$ in (ref).

We suppose the effects of $X_i$ on $Y_i$ are estimated by a partially linear model:

where $\beta$ and $g$ are defined as the minimizers of expected squared residuals $E[U_i^2]$:

for some linear space of functions $\mathcal{G}$. This setup nests linear covariate adjustment by setting $\mathcal{G}=\{\alpha+w^\prime\gamma\colon[\alpha, \gamma^\prime]^\prime\in\mathbb{R}^{1+\dim(W_i)}\}$, in which case (ref) gives a linear regression of $Y_i$ on ${X}_i$, $W_i$, and a constant. The setup also allows for more flexible covariate adjustments---such as by specifying $\mathcal{G}$ to be a large class of “nonparametric” functions robinson88.

Two examples highlight the generality of this setup:

(ref) nests the motivating \ac{RCT} example in (ref), allowing for an arbitrary number of experimental strata in $W_i$ and multiple treatment arms in $X_i$. (ref) shows that our setup can also nest the kind of regressions considered in a recent literature on \ac{DiD} and related regression specifications goodman2021difference,hull2018movers, sun2021estimating, de2020aer,de2020two,callaway2021difference, borusyak2021revisiting,wooldridge2021mundlak. We elaborate on the connections to this literature in (ref) by considering general \ac{TWFE} specifications with non-random treatments. These include specifications with multiple static treatment indicators, as in “mover regressions” that leverage over-time transitions, as well as dynamic event study specifications.\footnote{Some papers in this \ac{DiD} literature study issues we do not consider, such as when researchers fail to include indicators for all relevant treatment states, which will generally add bias terms to our decomposition of $\beta$, below. Similarly, we do not consider multicollinearity issues like in borusyak2021revisiting by assuming a unique solution to (ref). For event studies this means we assume some units are never treated, with $A(j)=\infty$.}

As a first step towards characterizing the treatment coefficient vector $\beta$, we solve the minimization problem in (ref). Let $\tilde{X}_i$ denote the residuals from projecting $X_i$ onto the control specification, with elements $\tilde{X}_{ik}=X_{ik}-\operatorname*{argmin}_{\tilde{g}\in\mathcal{G}} E[(X_{ik}-\tilde{g}(W_i))^{2}]$. It follows from the projection theorem vdv98 that

Applying the \ac{FWL} theorem, each treatment coefficient can be written $\beta_k=E[\accentset{\approx}{X}_{ik}Y_i]/E[\accentset{\approx}{X}_{ik}^2]$ where $\accentset{\approx}{X}_{ik}$ is the residual from regressing $\tilde{X}_{ik}$ on $\tilde{X}_{i,-k}=(\tilde{X}_{i1}, \dotsc, \tilde{X}_{i,k-1}, \tilde{X}_{i, k+1}, \dotsc, \tilde{X}_{iK})'$. Letting $E^{*}[X_{ik}\mid X_{i,-k}, W_i]$ denote the projection of $X_{ik}$ onto the space $\{X_{i,-k}'\tilde{\delta}+\tilde{g}(W_{i}) \colon \tilde{\delta}\in\mathbb{R}^{K-1}, \tilde{g}\in\mathcal{G}\}$, we may write these residuals as $\accentset{\approx}{X}_{ik}=X_{ik}-E^{*}[X_{ik}\mid X_{i,-k}, W_{i}]$.

Causal Interpretation

We now consider the interpretation of each treatment coefficient $\beta_k$ in terms of causal effects. Let $Y_i(k)$ denote the potential outcome of unit $i$ when $D_i=k$. Observed outcomes are given by $Y_i=Y_i(D_i)=Y_i(0)+X_i^\prime\tau_i$ where $\tau_i$ is a vector of treatment effects with elements $\tau_{ik}=Y_i(k)-Y_i(0)$. We denote the conditional expectation of the vector of treatment effects given the controls by $\tau(W_i)=E[\tau_i\mid W_i]$, so that $\tau_k(W_i)$ is the conditional \ac{ATE} for the $k$th treatment. We let $p(W_i)=E[X_i\mid W_i]$ denote the vector of propensity scores, so that $p_{k}(W_i) = \Pr(D_i=k\mid W_i)$. Our characterization of contamination bias doesn't require the propensity scores to be bounded away from $0$ and $1$ and in fact allows them to be degenerate, i.e.\ $p_k(w)\in\{0,1\}$ for all $w$. This is the case in (ref), since $X_i$ is a non-random function of $W_i$. We return to practical questions of propensity score support in (ref).

We make two assumptions to interpret $\beta_k$ in terms of the effects $\tau_i$. First, we assume mean-independence of the potential outcomes and treatment, conditional on the controls:

A sufficient condition for this assumption is that the treatment is randomly assigned conditional on the controls, making it conditionally independent of the potential outcomes:

Such conditional random assignment appears in (ref). In (ref), where treatment is a non-random function of the unit and time indices in $W_i$, (ref) holds trivially.

Second, we assume $\mathcal{G}$ is specified such that that one of two conditions holds:

The first condition requires the covariate adjustment to be flexible enough to capture each treatment's propensity score. For example, with a linear specification for $g$, (ref) requires the propensity scores to be linear in $W_i$ AnKr99handbook. This condition holds trivially in (ref), since $W_i$ is a vector of indicators for groups within which $X_i$ is randomly assigned. When this condition holds, the projection of the treatment onto the covariates coincides with the vector of propensity scores, and the projection residuals coincide with the conditionally demeaned treatment vector $\tilde{X}_i=X_i-p(W_i)$.

In (ref), with $X_i$ being a deterministic function of unit and time indices and $g(W_i)$ including unit and time fixed effects, (ref) fails because the propensity scores are binary---they cannot be captured by a linear combination of the \acp{TWFE}. However, (ref) is satisfied by a parallel trends assumption: that the average untreated potential outcomes $Y_i(0)$ are linear in the unit and time effects. We elaborate on this setup in (ref).\footnote{Identification based on (ref) can be seen as “design-based” in that it only restricts the treatment assignment process. Identification based on (ref) can be seen as “model-based” in that it makes no assumptions on the treatment assignment process but specifies a model for the unobserved untreated potential outcomes.}

Under either condition in (ref), the specification of controls is flexible enough to avoid \ac{OVB}. To see this formally, suppose all treatment effects are constant: $\tau_{ik}=\tau_k$ for all $k$. This restriction lets us write $Y_i=Y_i(0)+X_i^\prime\tau$, where $\tau$ is a vector collecting the constant effects. The only source of bias when regressing $Y_i$ on $X_i$ and controls is then the unobserved variation in the untreated potential outcomes $Y_i(0)$. But it follows from the expression for $\beta$ in (ref) that there is no such \ac{OVB} when (ref) holds:

Here the first equality uses the fact that $E[\tilde{{X}}_{i} X_i^\prime]=E[\tilde{X}_i\tilde{X}_i^\prime]$ because $\tilde{X}_i$ is a vector of projection residuals, and the second equality uses the law of iterated expectations and (ref). Under (ref), $E[\tilde{X}_{i} \mid W_i]=0$, so that the term in braces is zero by another application of the law of iterated expectations: $E[\tilde{X}_{i} E[Y_i(0)\mid W_i]]=E[E[\tilde{X}_{i} \mid W_i] E[Y_i(0)\mid W_i]]=0$. It is likewise zero under (ref) since $\tilde{X}_{i}$ is by definition of projection orthogonal to any function in $\mathcal{G}$ such that $E[\tilde{X}_{i} E[Y_i(0)\mid W_i]]=E[\tilde{X}_{i}\mu_{0}(W_{i})]=0$. Hence, \ac{OVB} is avoided in the constant-effects case so long as either the propensity scores or the untreated potential outcomes are spanned by the control specification. Versions of this double robustness property have been previously observed in, for instance, RoMaNe92.

When treatment effects are heterogeneous but ${X}_i$ contains a single treatment indicator, $\beta$ identifies a weighted average of the conditional effects $\tau(W_i)$. Specifically, since by the previous argument we still have $E[\tilde{X}_{i}Y_i(0)]=0$, it follows from (ref) that

where the second equality uses iterated expectations and the identity $E[\tilde{{X}}_{i}^2]=E[\tilde{X}_i{X}_i]$. Under (ref), $E[\tilde{X}_i X_i\mid W_i]=E[\tilde{X}_i^2\mid W_i]=\operatorname{var}(X_i\mid W_i)$, so the weights further simplify to $\lambda_{11}(W_i)=\frac{\operatorname{var}(X_i\mid W_i)}{E[\operatorname{var}(X_i\mid W_i)]}\ge 0$. This extends the angrist98 result to a general control specification; versions of this extension appear in, for instance, AnKr99handbook, AnPi09, and ArSa16.

This result provides a robustness rationale for estimating the effect of a single as-good-as-randomly assigned treatment with a partially linear model (ref): so long as the specification of $\mathcal{G}$ is rich enough to make (ref) hold, $\beta$ will identify a convex average of heterogeneous treatment effects. In (ref) we will derive another rationale for targeting $\beta$ in this model, showing that the weights $\lambda_{11}(W_i)$ minimize the semiparametric efficiency bound (conditional on the controls) for estimating some weighted-average treatment effect.

Our first proposition shows that with multiple treatments, the interpretation of $\beta$ becomes more complicated because of contamination bias:

(ref) shows that the coefficient on $X_{ik}$ in (ref) is a sum of two terms. The first term is a weighted average of conditional \acp{ATE} $\tau_{k}(W_i)$, with own treatment weights $\lambda_{kk}(W_i)$ that average to one---generalizing the characterization of the single-treatment case, (ref). The expression for $\lambda_{kk}$ implies that these weights are convex if the implicit linear probability model used to compute $\accentset{\approx}{X}_{ik}$ fits probabilities that lie below one, $E^{*}[X_{ik}\mid X_{i,-k}=0,W_{i}]\leq 1$. The second term is a weighted average of treatment effects for other treatments $\tau_{\ell}(W_i)$, with contamination weights $\lambda_{k\ell}(W_i)$ that average to zero. Because the contamination weights are zero on average, they must be negative for some values of the controls unless they are all identically zero.\footnote{(ref) complements an algebraic result in chatto_zubiz_2021, which shows that the regression estimator of $\beta_k$ can be written in terms of weighted sample averages of outcomes among units in different treatment arms (regardless of whether (ref) hold). In contrast, our analysis interprets regression estimands in terms of weighted averages of conditional \acp{ATE} under a broad class of identifying assumptions. In a finite-population setting, aaiw20 show that $\beta$ identifies matrix-weighted averages of individual treatment effect vectors $\tau_i$; however, they do not discuss the interpretation of the estimand.} This is the case when the implicit linear probability model correctly predicts that $X_{ik}=0$ if $X_{i\ell}=1$.

Hence, if the linear probability model is correctly specified, i.e. $E[X_{ik}\mid {X}_{i,-k}, W_i]=X_{i, -k}'\alpha+g_{k}(W_i)$ for some vector $\alpha$ and $g_{k}\in\mathcal{G}$, the contamination weights $\lambda_{k\ell}(W_{i})$ are zero and the own treatment weights $\lambda_{kk}(W_{i})$ are positive. This is the analog of condition (ref) if we interpret $X_{ik}$ as a binary treatment of interest and $X_{i, -k}'\alpha+g_{k}(W_i)$ as a specification for the controls. In other words, the assignment of treatment $k$ must be additively separable between $X_{i,-k}$ and $W_{i}$. However, with mutually exclusive treatments, this won't be the case unless treatment assignment is unconditionally random. In particular, since $X_{ik}$ must equal zero if the unit is assigned to one of the other treatments regardless of the value of $W_i$, under correct specification it must be the case that $\alpha_\ell=-g_k(W_i)$ for all elements $\alpha_{\ell}$ of $\alpha$. This in turn implies that the assignment of treatment $k$ doesn't depend on $W_i$, which is impossible unless the propensity score $p_k(W_i)$ is constant.

Thus, misspecification in the linear probability model will generally yield nonsensical fitted probabilities $E^{*}[X_{ik}\mid X_{i\ell}=1,W_{i}]\neq 0$ that generate non-zero contamination weights $\lambda_{k\ell}(W_{i})$. Furthermore, if the misspecification also yields fitted probabilities $E^{*}[X_{ik}\mid X_{i,-k}=0,W_{i}]> 1$, we will have negative own treatment weights. The last part of (ref) shows that such nonsensible predictions are ruled out if (ref) holds.

We make four further remarks on our general characterization of contamination bias:

Implications

(ref) shows that treatment effect heterogeneity can induce two conceptually distinct issues in flexible regression estimates of treatment effects. First, with either single or multiple treatments, there is a negative weighting of a treatment's own effects when projecting the treatment indicator onto other treatment indicators and covariates yields fitted values exceeding one, i.e.\ when $E^{*}[X_{ik}\mid X_{i,-k}=0,W_{i}]>1$. This issue is relevant in various \ac{DiD} regressions and related approaches which rely on a model of untreated potential outcomes that ensures (ref) holds (e.g.\ parallel trends assumptions) but which potentially misspecify the assignment model in (ref). Although the recent \ac{DiD} literature focuses on \ac{TWFE} regressions, (ref) shows such negative weighing can arise more generally---such as when researchers allow for linear trends, interacted fixed effects, or other extensions of the basic parallel trends model. None of these alternative specifications for $g$ are in general flexible enough to capture the degenerate propensity scores and hence ensure that $E^{*}[X_{ik}\mid X_{i,-k}=0,W_{i}]\leq 1$.

Second, in the multiple treatment case, there is a potential for contamination bias from other treatment effects---regardless of which condition in (ref) holds. This form of bias is relevant whenever one uses an additive covariate adjustment, no matter how flexibly the covariates are specified. Versions of this problem have been noted in, for example, the sun2021estimating analysis of \ac{DiD} regressions with treatment leads and lags or the hull2018movers analysis of mover regressions (see (ref)).\footnote{The negative weights issue raised in de2020aer (when $K=1$), and the related issue that own-treatment weights may be negative in sun2021estimating and de2020two (when $K>1$), arise because the treatment probability is not linear in the unit and time effects. If (ref) holds with $K=1$, (ref) shows $\beta$ estimates a convex combination of treatment effects. This covers the setting considered in Theorem 1(iv) in athey2022design. In their Comment 2, athey2022design say that “the sum of the weights [used in Theorem 1(iv)] is one, although some of the weights may be negative”. (ref) shows these weights are, in fact, non-negative.} (ref) shows such contamination bias arises much more broadly, however.

The characterization in (ref) also relates to concerns in interpreting multiple-treatment \ac{IV} estimates with heterogeneous effects behaghel2013robustness,kirkeboen2016field,kline2016headstart,hull2018isolateing,leeMultivalued,bhuller_sigstad. This connection comes from viewing (ref) as the second stage of an \ac{IV} model estimated by a control function approach; in the linear \ac{IV} case, for example, $g(W_i)$ can be interpreted as giving the residuals from a first-stage regression of $X_i$ on a vector of valid instruments $Z_i$. In the single-treatment case, the resulting $\beta$ coefficient has an interpretation of a weighted average of conditional local average treatment effects under the appropriate first-stage monotonicity condition ia94. But as in (ref) this interpretation fails to generalize when $X_i$ includes multiple mutually-exclusive treatment indicators: each $\beta_k$ combines the local effects of treatment $k$ with a non-convex average of the effects of other treatments.

Finally, (ref) has implications for single-treatment \ac{IV} estimation with multiple instruments and flexible controls if the first stage has the form of (ref), where now $Y_i$ is interpreted as the treatment and $X_i$ gives the vector of instruments. (ref) shows that the first-stage coefficients on the instruments $\beta_k$ will not generally be convex weighted average of the true first-stage effects $\tau_{ik}$. Because of this non-convexity, the regression specification may fail to satisfy the effective monotonicity condition even when $\tau_{ik}$ is always positive: the cross-instrument contamination of causal effects may cause monotonicity violations, even when specifications with individual instruments do not. This issue is distinct from previous concerns over monotonicity failures in multiple-instrument designs mueller2015criminal, frandsen2019judging, norris2019examiner, mogstad2021identification, which are generally also present in such just-identified specifications. It is also distinct from concerns about insufficient flexibility in the control specification when monotonicity holds unconditionally blandhol2022tsls.

This new monotonicity concern may be especially important in “examiner” \ac{IV} designs, which exploit the conditional random assignment to multiple decision-makers. Many studies leverage such variation by computing average examiner decision rates, often with a leave-one-out correction, and use this “leniency” measure as a single instrument with linear controls. These \ac{IV} estimators can be thought of as implementing versions of a jackknife \ac{IV} estimator aik99, based on a first stage that uses examiner indicators as instruments, similar to (ref). (ref) thus raises a new concern with these \ac{IV} analyses when controls (such as time fixed effects) are needed to ensure ignorable treatment assignment.

Solutions

We now discuss three solutions to the contamination bias problem raised by (ref), each targeting a distinct causal parameter. First, in (ref), we discuss estimation of unweighted \acp{ATE}. The other two solutions target weighted averages of individual treatment effects using an \acf{EW} scheme in that the weights minimize the semiparametric efficiency bound for estimating weighted \acp{ATE} under homoskedasticity. In the second solution, the weights are allowed to vary across treatments, while in the third, they are constrained to be common across treatments. In (ref) we characterize these estimation targets, while in (ref) we discuss how to estimate them; we also outline our proposed guidance to researchers in measuring contamination bias.

Implementing the first solution requires strong overlap (i.e.\ that treatment propensity scores are bounded away from zero and one) while the other two solutions require nonempty overlap, ruling out fully degenerate propensity scores. Solutions allowing for degenerate propensity scores require either targeting subpopulations of the treated or adding substantive restrictions on conditional means of treated potential outcomes (beyond (ref), which only restricts untreated potential outcomes). We refer readers to de2020two, sun2021estimating, callaway2021difference,borusyak2021revisiting,wooldridge2021mundlak for such solutions in the context of \ac{DiD} regressions.

Estimating Average Treatment Effects

Many estimators exist for the \ac{ATE} of binary treatments---see imbens2009recent and abadie2018econometric for reviews. Several of these approaches extend naturally to multiple treatments: including matching on covariates or the propensity score, inverse propensity score weighting, balancing weights, interacted regression, or doubly-robust methods (see, among others, cattaneo10, ReZu20, ChNeSi21, and GrPi22). Here we summarize the last two approaches.

For the interacted regression solution, we adapt the implementation for the binary treatment case discussed in imbens2009recent to multiple treatments. Specifically, consider the specification:

where $q_{k}\in\mathcal{G}$, $k=0,\dotsc, K$ and we continue to define $\beta$ and the functions $q_k$ as minimizers of $E[\dot{U}_i^2]$. When $\mathcal{G}$ consists of linear functions, (ref) specifies a linear regression of $Y_i$ on $X_i$, $W_i$, a constant, and the interactions between each treatment indicator $X_{ik}$ and the demeaned control vector $W_i-E[W_i]$. Define $\mu_{k}(w)=E[Y_{i}(k)\mid W_{i}=w]$ for $k=0,\dotsc, K$, so that $\tau_{k}(w)=\mu_{k}(w)-\mu_{0}(w)$. If (ref) holds and $\mathcal{G}$ is furthermore rich enough to ensure $\mu_{k}\in\mathcal{G}$ for $k=0,\dotsc, K$ then $\beta=\tau$. Moreover, $q_{k}(w)=\tau_{k}(w)$ for $k=1,\dotsc, K$, such that the regression identifies both the unconditional and conditional \acp{ATE}.

The added interactions in (ref) ensure that each treatment coefficient $\beta_{k}$ is determined only by the outcomes in treatment arms with $D_{i}=0$ and $D_{i}=k$, avoiding the contamination bias in (ref). Demeaning the $q_k(W_i)$ in the interactions ensures they are appropriately centered to interpret the coefficients on the uninteracted $X_{ik}$ as \acp{ATE}.

Estimation of (ref) is conceptually straightforward for parametric $q_{k}$. In particular, if $\mathcal{G}$ consists of linear functions, one simply estimates

by \ac{OLS}, where $\mkern 1.5mu\overline{\mkern-1.5muW\mkern-1.5mu}\mkern 1.5mu=\frac{1}{N}\sum_i W_i$ is the sample average of the covariate vector. More generally, to increase the plausibility of the key assumption that $\mu_{k}\in\mathcal{G}$, one may constrain $\mathcal{G}$ only by nonparametric smoothness assumptions. Given a sequence of basis functions $\{b_{j}(W_{i})\}_{j=1}^{\infty}$, such as polynomials or splines, one then approximates $q_{k}$ with a linear combination of the first $J$ terms, with $J$ increasing with the sample size, thus tailoring the model complexity to data availability. Given a choice of $J$, estimation and inference can proceed as in the parametric case; the only difference is that the baseline covariates $W_{i}$ in (ref) are replaced by the basis vector $(b_{1}(W_{i}), \dotsc, b_{J}(W_{i}))'$ and $\mkern 1.5mu\overline{\mkern-1.5muW\mkern-1.5mu}\mkern 1.5mu$ is replaced by the sample average of this expansion. This estimator has been studied in the binary treatment case by ChHoTa08 and ImNeRi07, with the latter providing a detailed analysis of how to choose $J$ and the former showing that this sieve estimator achieves the semiparametric efficiency bound under strong overlap: it is impossible to construct another regular estimator of the \ac{ATE} with smaller asymptotic variance.

An attractive alternative approach combines the interacted regression with inverse propensity score weighting. Instead of using \ac{OLS} to estimate (ref) one uses weighted least squares, weighting observations by the inverse of some estimate $\hat{p}_{D_{i}}(W_{i})$ of the propensity score (see, e.g., robins1994estimation, wooldridge2007inverse, sloczynski2018general). An advantage of this approach is that it is doubly-robust: the estimator is consistent so long as either the propensity score estimator is consistent or the outcome model is correct (i.e. $\mu_{k}\in\mathcal{G}$). A recent literature shows how the double robustness property, when combined with cross-fitting, reduces the sensitivity of the ATE estimate to overfitting or regularization bias in estimating the nuisance functions $p_{k}$ and $\mu_{k}$. Cross-fitting also allows for using more flexible methods to approximate $p_{k}$ and $\mu_{k}$, including modern machine learning methods chernozhukov2018double,ceinr22,ChNeSi21.

Either approach should work reliably in stratified \acp{RCT} and other settings with strong overlap. But under weak overlap, when propensity scores are not bounded away from zero and one, all of these \ac{ATE} estimators may be imprecise and have poor finite-sample behavior. This is not a shortcoming of the specific estimator; indeed, KhTa10 show that under weak overlap, $\sqrt{N}$-estimation of the \ac{ATE} is not possible. Furthermore, if some propensity scores attain values of zero or one, the \ac{ATE} is not even point-identified. These results formalize the intuition that it is difficult or impossible to estimate the counterfactual outcomes for units with extreme propensity scores.\footnote{One approach to limited overlap is trimming: i.e., dropping observations with extreme propensity scores crump2006dealing,crump2009dealing, yang2016propensity. As with the estimators we derive next, trimming estimators shift the estimand from \ac{ATE} to easier-to-estimate weighted averages of conditional \acp{ATE}.} Such extreme propensity scores are common in observational settings. The solutions we discuss next downweight these difficult-to-estimate counterfactuals to address this practical challenge.

Easiest-to-Estimate Averages of Treatment Effects

Suppose in a sample of observations $i=1, \dots, N$ we wish to estimate a weighted average of conditional potential outcome contrasts $\sum_{i=1}^{N}\lambda(W_{i}) \sum_{k=0}^{K}c_{k}\mu_{k}(W_{i})/\sum_{i=1}^{N}\lambda(W_{i})$, where $\mu_{k}(W_{i})=E[Y_{i}(k)\mid W_{i}]$, $c$ is a $(K+1)$-dimensional contrast vector with elements $c_{k}$, and $\lambda(W_i)$ is some weighting scheme.\footnote{In a slight abuse of notation relative to (ref), the weights $\lambda$ here are not required to average to one. Instead, we scale the estimand by the sum of the weights, $\sum_{i=1}^{N}\lambda(W_{i})$.} We focus on two specifications for the contrast vector, leading to two alternatives the \ac{ATE} target. First, for separately estimating the effect of each treatment $k$, we set $c_{k}=1$, $c_{0}=-1$ and set the remaining entries of $c$ to $0$. The contrast of interest then becomes $\sum_{i=1}^{N}\lambda(W_{i})\tau_{k}(W_{i})/\sum_{i=1}^{N}\lambda(W_{i})$, the weighted \ac{ATE} of treatment $k$. Second, we specify $c$ so as to allow us to simultaneously contrast the effects of all $K$ treatments---we discuss this further below. For each contract vector $c$, we characterize in this section the \acf{EW} scheme $\lambda(W_i)$ that leads to the smallest possible standard errors under homoskedasticity. We discuss estimation of the corresponding estimands in (ref).

This optimization problem has four motivations. First, there is a robustness motivation: a researcher would like to estimate a given contrast as precisely as possible, at least under the benchmark of constant treatment effects, while being robust to the possibility that the effects are heterogeneous. While the optimization problem does not impose convexity, it turns out that the \ac{EW} scheme is convex. Hence, the resulting estimand identifies a convex average of conditional contrasts under heterogeneous treatment effects, and avoids any contamination bias. Such a robustness property presumably underlies the popularity of regression as a tool for estimating the effect of a binary treatment: the regression estimator is efficient under homoskedasticity and constant treatment effects while, by the angrist98 result, retaining a causal interpretation under heterogeneous effects.\footnote{There are several motivations for the interest in convex weights. First, $\lambda(W_i)\ge 0$ ensures the estimand captures average effects for some well-defined (and characterizable) subpopulation. Second, it prevents what strlb17 call a sign-reversal: if $\tau_k(w)$ has the same sign for all $w$ ($+,0$ or $-$), then the estimand will also have this sign. blandhol2022tsls call such estimands “weakly causal.” Finally, the estimand satisfies a population version of what rslr07 call boundedness: the estimand lies in the support of $\tau_k(w)$.}

Second, the \ac{EW} scheme gives a bound on the information available in the data: if the scheme yields overly large standard errors, inference on other treatment effects (such as the unweighted \ac{ATE}) must be at least as uninformative. Computing the \ac{EW} standard errors thus reveals whether informative conclusions for any treatment effect estimand are only possible under additional assumptions or with the aid of additional data. In fact, we show below that in the binary treatment case the \ac{EW} scheme is exactly the same as that used by regression. Recall that in the binary treatment case, the regression treatment weights are proportional to the conditional variance of treatment, $\operatorname{var}(D_i\mid W_i)=p_{1}(W_i)(1-p_{1}(W_i))$. Because these weights tend to zero as $p_{1}(W_{i})$ tends to zero or one, regression downweights observations with extreme propensity scores where the estimation of counterfactual outcomes is difficult, avoiding the poor finite-sample behavior of \ac{ATE} estimators under weak overlap and allowing for informative inference even when one cannot precisely estimate the unweighted \ac{ATE}.

Third, the \ac{EW} scheme can be viewed as offering an intermediate point along a particular robustness-precision “possibility frontier.” The \ac{ATE} estimator based on the interacted specification in (ref) lies on one end of this frontier, being the most robust to treatment effect heterogeneity (i.e.\ retaining a clear interpretation regardless of the form of $\tau(w)$ or how it relates to the propensity scores). But this robustness comes at the cost of imprecision and non-standard inference under weak overlap. The regression estimator based on (ref) lies on the other end of the frontier: it is likely to be precise even when overlap is weak (and is efficient under homoskedasticity if the partly linear model in (ref) is correct, such that treatment effects are constant). But this precision comes at the cost of contamination bias under heterogeneous treatment effects. The \ac{EW} scheme lies in between these extremes, purging contamination bias and retaining good performance under weak overlap by giving up explicit control over the treatment effect weighting, letting it be data-determined.\footnote{There are other approaches to resolving the robustness-precision tradeoff, such as seeking precise estimates subject to the weights $\lambda$ remaining “close” to one, or placing some restrictions on the form of effect heterogeneity, in contrast to leaving it completely unrestricted as we do here (see mst18 for an example of this approach in an \ac{IV} setting). We leave these alternatives to future research.}

Finally, while the derivation of the \ac{EW} scheme is motivated by statistical precision concerns, the resulting estimand can be seen as identifying the impact of a policy that manipulates the treatment via a particular incremental propensity score intervention. We discuss this interpretation in (ref) below.

We derive the \ac{EW} scheme in two steps. First, we establish a precision benchmark---a semiparametric efficiency bound---for estimation of a given weighted average of treatment effects under the idealized scenario that the propensity score is known. Second, we determine which weights $\lambda$ minimize the bound.

The following \namecref{theorem:seb} establishes the first step of our derivation:

As formalized in the (ref) proof, $\mathcal{V}_{\lambda, c}$ establishes the lower bound on the asymptotic variance of any regular estimator of $\theta_{\lambda, c}$ under the idealized case of known propensity scores.\footnote{The efficiency bound for the population analog $\theta^*_{\lambda, c}=E[\lambda(W_i)\sum_{k=0}^{K}c_{k}\mu_{k}(W_{i})]/E[\lambda(W_i)]$ has an additional term, $E[\lambda(W_{i})^{2}(\sum_{k=0}^{K}c_{k}\mu_{k}(W_{i})-\theta_{\lambda, c}^{*})^{2}]/E[\lambda(W_{i})]^{2}$, reflecting the variability of the conditional average contrast. The variance-minimizing weights for $\theta^*_{\lambda, c}$ thus depend on the nature of treatment effect heterogeneity. By focusing on $\theta_{\lambda, c}$, we avoid this term, which allows us give the characterization in (ref) without any assumptions about heterogeneity in treatment effects.}

To establish the second step, we minimize (ref) over $\lambda$. Simple algebra shows that the \ac{EW} scheme is (up to an arbitrary constant) given by

Observe that this scheme delivers convex weights, $\lambda^*_{c}\geq 0$, even though convexity was not imposed in the optimization. Hence, there is no cost in precision if we restrict attention to convex weighted averages of conditional \acp{ATE}.

When the contrast vector is selected to estimate the weighted average effect of a particular treatment $k$, a corollary to (ref) is that regression weights are the easiest-to-estimate:

Per (ref), the weighting $\lambda^k$ coincides with the weighting of conditional \acp{ATE} from the partially linear model (ref) when it is fit only on observations with $D_{i}\in\{0,k\}$, provided $p_{k}/(p_{k}+p_{0})\in\mathcal{G}$.\footnote{This follows since the propensity score in the subsample is given by $\Pr(D_{i}=k\mid W_{i}, D_{i}\in\{0,k\})=\frac{p_{k}(W_{i})}{p_{0}(W_{i})+p_{k}(W_{i})}$, so that $\lambda^k(W_i)$ in (ref) equals the conditional variance of the treatment indicator times the probability of being in the subsample.} (ref) thus gives a precision justification for estimating the effect of any given treatment $k$ by a partially linear regression in the subsample with $D_{i}\in\{0,k\}$ under a homoskedasticity benchmark, complementing the robustness motivation discussed earlier.\footnote{As usual, homoskedasticity is a tractable baseline: the arguments in favor of \ac{OLS} following (ref) can be extended to favor a (feasible) weighted least squares regression when $\sigma^2(W_i)$ is consistently estimable.} To estimate the effects of all treatments one can run $K$ such one-treatment-at-a-time regressions, one for each treatment arm. Plugging (ref) into (ref) reveals that the asymptotic variance is bounded so long as the overlap between the covariate distribution in each treatment arm is nonempty, i.e. $P(p_{k}(W_{i})>\varepsilon\cap p_{0}(W_{i})>\varepsilon)>\varepsilon$ for some $\varepsilon$.

For binary treatments, crump2006dealing and liBalancing show that the weighting $p_1(W_i)(1-p_{1}(W_i))$ minimizes the asymptotic variance of a particular class of inverse propensity score weighted estimators. Our (ref) extends the property to all regular estimators, and to multiple treatments.

A shortcoming of the \ac{EW} scheme in (ref) is that it is treatment-specific, precluding comparisons of the weighted-average effects across treatments.\footnote{Formally, for treatments $1$ and $2$, we estimate the weighted averages $\sum_{i}\lambda^{1}(W_i)\tau_{1}(W_i)/\sum_{i}\lambda^{1}(W_i)$ and $\sum_{i}\lambda^{2}(W_i)\tau_{2}(W_i)/\sum_{i}\lambda^{2}(W_i)$. Because the weights $\lambda^1$ and $\lambda^2$ differ, the difference between these estimands cannot generally be written as a convex combination of conditional treatment effects $\tau_{1}(W_i)-\tau_{2}(W_i)$. This critique also applies to the own-treatment weights in (ref). Thus even without contamination bias one may find the implicit multiple-treatment regression weighting deficient.} This issue is especially salient when the control group is arbitrarily chosen, such as in teacher \ac{VAM} regressions which omit an arbitrary teacher from estimation and seek causal comparisons across all teachers.

We thus turn to the question of how (ref) can be used to select a weighting scheme which allows for simultaneous comparisons across all treatment arms. Suppose that the contrast of interest is drawn at random from a given marginal treatment distribution $\Pr(D_i=k)=\pi_k$, so that $c_{j}=1$ with probability $\pi_{j}(1-\pi_{j})/(1-\sum_{k=0}^{K}\pi_{k}^{2})$ and $c_{j}=-1$ with the same probability.\footnote{Formally, we draw two treatments at random from the given marginal distribution, discarding the draw if the two treatments are equal.} Let $F_\pi$ denote this distribution over the (now random) contrasts. If the researcher wishes to report an accurate contrast estimate but needs to commit to a weighting scheme before knowing the contrast of interest, it is optimal to minimize the expected variance

Minimizing this expression over $\lambda$ is equivalent to minimizing (ref) with $c_{k}^{2}=2\pi_k(1-\pi_k)$, which yields (ref) with this contrast specification as the optimal weighting. Thus, the optimal weights are proportional to $\left(\sum_{k=0}^{K}\frac{\pi_{k}(1-\pi_{k})\sigma^{2}_{k}(W_{i})}{p_{k}(W_{i})}\right)^{-1}$. Specializing to the homoskedastic case leads to the following result:

The \ac{CW} scheme $\lambda^{\textnormal{CW}}$ generalizes the intuition behind the single binary treatment ((ref)), placing lower weight on strata with extreme propensity scores. When the treatment is binary, $K=1$, the $\pi_{k}$'s do not matter and the \ac{CW} scheme reduces to that in (ref): $\lambda^{\textnormal{CW}}(W_{i})=\lambda^{1}(W_{i})=\lambda^{0}(W_{i})=p_{1}(W_{i})p_{0}(W_{i})$. With multiple treatments, however, the weights $\lambda^{\textnormal{CW}}$ remain the same for every treatment---allowing for simultaneous comparisons across all treatment pairs $(k, \ell)$.

There are two natural choices for the marginal treatment probabilities $\pi$. First, when equally interested in all contrasts, one can set $\pi_k=1/(K+1)$. This weighting scheme was previously proposed by lili19; our characterization of it in terms of optimizing a semiparametric efficiency bound is, to our knowledge, novel. Second, if more common treatments are of greater interest, we may set $\pi_k$ to the empirical treatment probabilities $N^{-1}\sum_i X_{ik}$. This weighting targets precise estimation of contrasts involving more common treatments at the expense of contrasts involving less common treatments. We use this choice in our empirical applications in (ref). For either choice of weights, the resulting asymptotic variance in (ref) remains bounded so long as the overlap between covariate distributions in each treatment arm is not empty: $P(\cap_{k=0}^{K}p_{k}(W_{i})>\varepsilon)>\varepsilon$ for some $\varepsilon$. Non-empty overlap is a substantially weaker assumption than strong overlap, needed for $\sqrt{N}$-estimation of the unweighted \ac{ATE}, which requires this probability to equal one. For instance, in the nine empirical applications below, non-empty overlap always holds, but strong overlap fails in six.

Practical Guidance in Measuring and Avoiding Contamination Bias

A researcher interested in estimating the effects of multiple mutually exclusive treatments with regression can use (ref) to measure the extent of contamination bias in their estimates. When the propensity score is not fully degenerate, they can further estimate one of the alternative estimation targets discussed in the previous subsections. Here we provide practical guidance on both procedures, which we illustrate empirically in the next section.

For simplicity, we focus on the case where $g$ is linear and (ref) is estimated by \ac{OLS}. We suppose (ref) and both conditions in (ref) hold, such that all propensity scores $p_k$ and potential outcome conditional expectation functions $\mu_k$ are linearly spanned by the controls $W_i$. These conditions hold, for example, when $W_i$ contains a set of mutually exclusive group indicators. When $\mathcal{G}$ is unrestricted, the recommendations in this section would require non-parametric approximations for $g$ analogous to those discussed in (ref).

Under this setup, we can decompose the \ac{OLS} estimator $\hat{\beta}$ from the uninteracted regression

to obtain a sample analog of the decomposition in (ref). To this end, note that the own-treatment and contamination bias weights in (ref) are identified by the linear regression of $X_i$ on the residuals $\tilde{X}_i$. Specifically, $\lambda_{k\ell}(W_{i})$ is given by the $(k, \ell)$th element of the $K\times K$ matrix $\Lambda(W_{i})= E[\tilde{X}_{i}\tilde{X}_{i}']^{-1}E[\tilde{X}_{i}X_{i}'\mid W_{i}]$, which can be estimated by its sample analog $\hat{\Lambda}_{i}=(\dot{X}'\dot{X})^{-1}\dot{X}_{i}X_{i}',$ where $\dot{X}_{i}$ is the sample residual from an \ac{OLS} regression of $X_{i}$ on $W_{i}$ and a constant and $\dot{X}$ is a matrix collecting these sample residuals. The $(k, \ell)$th element of $\hat\Lambda_i$ estimates the weight that observation $i$ puts on the $\ell$th treatment effect in the $k$th treatment coefficient. For $k=\ell$ this is an estimate of the own-treatment weight in (ref); for $k\neq\ell$ this is an estimate of a contamination weight.

Under linearity, the $k$th conditional \ac{ATE} may be written as $\tau_k(W_{i})=\gamma_{0,k}+W_{i}^\prime\gamma_{W, k}$, where $\gamma_{0,k}$ and $\gamma_{W, k}$ are coefficients in the interacted regression specification

Estimating (ref) by \ac{OLS} yields estimates $\hat{\tau}_k(W_i)=\hat\gamma_{0,k}+W_{i}'\hat\gamma_{W, k}$. For each observation $i$, we stack the set of conditional \ac{ATE} estimates in a $K\times 1$ vector $\hat\tau (W_i)$.

Using the \ac{OLS} normal equations, we then obtain a sample analog of the population decomposition in (ref):

The first term estimates the own-treatment effect components, $E[\lambda_{kk}(W_{i})\tau_{k}(W_{i})]$, while the second term estimates the contamination bias components, $\sum_{\ell\neq k}E[\lambda_{k\ell}(W_{i})\tau_{\ell}(W_{i})]$. If the contamination bias term is large for some $\hat\beta_k$, it suggests the estimate of the $k$th treatment effect is substantially impacted by the effects of other treatments. Researchers can also compare the first term of (ref) to other weighted averages of own-treatment effects, including the ones discussed next, to gauge the impact of the regression weighting $\operatorname{diag}(\hat{\Lambda}_{i})$.\footnote{When the covariates are not saturated, it is possible that the estimated weighting function $\hat{\Lambda}(w)=\frac{1}{N}\sum_{i=1}^{N}\1{W_{i}=w}\hat{\Lambda}_{i}$ is not positive-definite for some or all $w$. In particular, the diagonal elements of $\hat{\Lambda}(w)$ need not all be positive. However, it is guaranteed that the diagonal of $\hat{\Lambda}(w)$ sums to one and the non-diagonal weights sum to zero, since $\sum_{i=1}^{N}\hat{\Lambda}_{i}=I_{k}$.}

Further analysis of the estimated weights $\hat\lambda_{k\ell}(w)=\frac{\sum_{i=1}^{N}\1{W_i=w}\hat{\Lambda}_{i, k\ell}}{\sum_{i=1}^{N}\1{W_i=w}}$ can shed more light on the regression estimates in $\hat\beta$. For example, the contamination weights for $\ell\neq k$ can be plotted against the treatment effect estimates $\hat\tau_\ell(W_i)$ to visually assess the sources of contamination bias. Low bias may arise from limited treatment effect heterogeneity, small contamination weights, or a low correlation between the two.

Estimation of the unweighted \ac{ATE} and the \ac{EW} and \ac{CW} schemes is also straightforward under the linearity assumptions. First, estimating (ref) by \ac{OLS} yields estimates of the unweighted \acp{ATE} $\tau_k=E[\tau_k(W_i)]$. The estimates are numerically equivalent to $\hat{\tau}_k=\hat\gamma_{0,k}+\mkern 1.5mu\overline{\mkern-1.5muW\mkern-1.5mu}\mkern 1.5mu^\prime\hat\gamma_{W, k}$, where $\hat\gamma_{0,k}$ and $\hat\gamma_{W, k}$ are \ac{OLS} estimates of (ref).

Second, the \ac{EW} scheme from (ref) can be estimated using the uninteracted one-treatment-at-a-time regression

where we only use observations assigned either to treatment $k$ or the control group.

The third solution is to estimate the \ac{CW} scheme $\lambda^{\textnormal{CW}}$ from (ref). We use inverse propensity score weighting in our applications below: we regress $Y_{i}$ onto $X_{i}$ and a constant, weighting each observation by $\hat{\lambda}^{\textnormal{CW}}(W_i)/\hat{p}_{D_i}(W_i)$ where $\hat{p}_{k}(W_{i})$ denotes estimated propensity scores from a multinomial logit model and

is an estimate of $\lambda^{\textnormal{CW}}$. When the weights $\pi$ are uniform, this estimator reduces to the estimator studied in lili19. The resulting estimator can be written as

When the treatment is binary and $\hat{p}_{k}$ is obtained via a linear regression, this weighted regression estimator coincides with the usual (unweighted) regression estimator that regresses $Y_{i}$ onto $D_{i}$ and $W_{i}$.\footnote{To see this, note that in this case $\hat{\lambda}(W_{i})=\hat{p}_{1}(W_{i})\hat{p}_{0}(W_{i})$, so that $\hat{\beta}_{\hat{\lambda}^{\textnormal{CW}}, 1}= \frac{\sum_{i=1}^{N}(1-\hat{p}_{1}(W_{i}))D_{i}Y_{i}}{ \sum_{i=1}^{N}(1-\hat{p}_{1}(W_{i}))D_{i}} - \frac{\sum_{i=1}^{N}\hat{p}_{1}(W_{i})(1-D_{i})Y_{i}}{\sum_{i=1}^{N}\hat{p}_{1}(W_{i}) (1-D_{i})} = \frac{\sum_{i=1}^{N}(D_{i}-\hat{p}_{1}(W_{i}))Y_{i}}{ \sum_{i=1}^{N}(D_{{i}}-\hat{p}_{1}(W_{i}))^{2}} $, where the second equality uses the least-squares normal equations $\sum_{i=1}^{N}X_{i1}=\sum_{i=1}^{N}\hat{p}_1(W_i)$ and $\sum_{i}X_{i1}\hat{p}_1(W_i)=\sum_{i=1}^{N}\hat{p}_1(W_i)^2$.} (ref) in (ref) shows that the estimator $\hat{\beta}_{\hat{\lambda}^{\textnormal{CW}}}$ is efficient in the sense that it achieves the semiparametric efficiency bound for estimating $\beta_{\lambda^{\textnormal{CW}}}=\sum_i \lambda^{\textnormal{CW}}(W_i)\tau(W_i)/\sum_i \lambda^{\textnormal{CW}}(W_i)$.

Applications

Project STAR Application

We first illustrate our framework for analyzing and addressing contamination bias with data from Project STAR, as studied in krueger1999experimental. The Project STAR \ac{RCT} randomized 11,600 students in 79 public Tennessee elementary schools to one of three types of classes: regular-sized (20--25 students), small (target size 13--17 students), or regular-sized with a teaching aide. The proportion of students randomized to the small class size and teaching aide treatment varied over schools, due to school size and other constraints on classroom organization. Students entering kindergarten in the 1985--1986 school year participated in the experiment through the third grade. Other students entering a participating school in grades 1--3 during these years were similarly randomized between the three class types. We focus on kindergarten effects, where differential attrition and other complications with the experimental analysis are minimal.\footnote{Students in regular-sized classes were randomly reassigned between classrooms with and without a teaching aide after kindergarten, complicating the interpretation of the aide effect in later grades. The randomization of students entering the sample after kindergarten was also complicated by the uneven availability of slots in small and regular-sized classes krueger1999experimental.}

Column 1 of Panel A in (ref) reports estimates of kindergarten treatment effects in a sample of 5,868 students initially randomized to the small class size and teaching aide treatments. Specifically, we estimate the partially linear regression ((ref)) where $Y_i$ is student $i$'s test score achievement at the end of kindergarten, $X_{i}=(X_{i1}, X_{i2})$ are indicators for the initial experimental assignment to a small kindergarten class and a regular-sized class with a teaching aide, respectively, and $W_i$ is a vector of school fixed effects. We follow krueger1999experimental in computing $Y_i$ as the average percentile of student $i$'s math, reading, and word recognition score on the Stanford Achievement Test in the experimental sample. As in the original analysis krueger1999experimental, we obtain a small class size effect of 5.36 with a heteroskedasticity-robust standard error of 0.78 and a teaching aide effect of 0.18 (standard error: 0.72).\footnote{Our sample and estimates are very similar to---but not exactly the same as---those in krueger1999experimental. We use heteroskedasticity-robust (non-clustered) standard errors throughout this analysis, since the randomization of students to classrooms is at the individual level.}

As discussed in (ref), treatment assignment probabilities vary across the schools indicated by the fixed effects in $W_i$. If treatment effects also vary across schools in a way that covaries with the contamination weights $\lambda_{k\ell}(W_i)$, we expect the estimated effect of small class sizes to be partly contaminated by the effect of teaching aides (and vice versa). Panel B reports the contamination bias part of the decomposition in (ref), which appears minimal for both treatment arms.

It is useful to decompose the contamination bias further into the standard deviation of the school-specific treatment effect $\tau_{\ell}(W_{i})$, standard deviation of the contamination weights, and their correlation, as discussed in (ref). (ref) in (ref) does this graphically, plotting estimates of the school-specific treatment effects $\tau_\ell(W_i)$ against the contamination weights $\lambda_{k\ell}(W_i)$ for $\ell\neq k$. As can be seen from (ref), the variability of school-specific treatment effects is substantial: Adjusting for estimation error, we estimate the standard deviation of $\tau_k(W_i)$ to be 11.0 for the small class treatment and of 9.1 for the aide treatment.\footnote{We adjust for estimation error by subtracting the average squared standard error from the empirical variance of the treatment effect estimates and taking the square root.} Both standard deviations are an order of magnitude larger than the standard errors in (ref). On the other hand, the standard deviations for the contamination weights for the small class and aide treatment are only moderate: $0.14$ and $0.11$, respectively. Moreover, the correlation between the conditional treatment effects and the contamination weights is weak: $0.10$ for the small class effect estimate and $-0.13$ for the aide effect estimate. The moderate variation in the contamination weights coupled with weak correlation between the weights and the treatment effects explains why the contamination bias is small, even though the treatment effects vary substantially across schools.

Had the experimental design been such that the contamination weights strongly correlate with the treatment effects, sizable contamination bias could have resulted. To illustrate this, we compute worst-case (positive and negative) weighted averages of the estimated $\tau_\ell(W_i)$ by re-ordering them across the computed cross-treatment weights $\lambda_{k\ell}(W_i)$. This exercise highlights potential scenarios in which the randomization strata happened to have been highly correlated with the effect heterogeneity. Columns 2 and 3 in panel B of (ref) show that both bounds on possible contamination bias are an order of magnitude larger than the actual contamination bias: $[-1.65,1.67]$ for the small class size treatment and $[-1.53,1.53]$ for the teaching aide treatment.\footnote{The point estimates and standard errors in columns 4 and 5 in (ref) do not account for the fact that the re-ordering is based on estimates of $\tau_{k}(W_i)$ rather than the true treatment effects. This biases the reported estimates away from zero, so that they give an upper bound for the worst-case contamination bias.} Overall, for both treatments, the underlying heterogeneity in this setting makes substantial contamination bias possible even though actual contamination bias turns out to be relatively small.

Columns 2--5 of panel A report four treatment effect estimates that are free of contamination bias. Column 2 gives the own-treatment effect component of the decomposition in (ref), netting out the contamination bias estimate from column 1. This doubles the teaching aide effect estimate, from 0.18 to 0.36, but the estimate remains statistically insignificant with standard errors of around 0.71; the small classroom estimate moves very little. The remaining columns report the three solutions to contamination bias discussed in (ref). Column 3 estimates the unweighted \acp{ATE} of the small class size and teaching aide treatment, by estimating the interacted regression specification in (ref). Column 4 estimates the one-treatment-at-a-time regressions in (ref) for $k=1,2$. Finally, column 5 runs a weighted regression of $Y_i$ onto $X_i$ using the \ac{CW} scheme in (ref).

There turns out to be little difference between these alternative estimates. The small class size effect varies between 5.2 and 5.6, which is close to the original estimate. The teaching aide effect varies between 0.01 and 0.26. To understand this lack of variation, recall from (ref) that the difference between the unweighted \ac{ATE} and an estimand that uses weights $\lambda(W_{i})$ is given by the covariance between $\lambda(W_{i})$ and the conditional \acp{ATE} $\tau_k(W_{i})$. Given the sizable variability in the treatment effect estimates, the covariance will be small only if the correlation between the weights and the treatment effects is small and if the weights display limited variability. This turns out to be the case here, as depicted graphically in (ref) in (ref). The figure shows that the correlations fall below $0.25$ in absolute value for all weighting schemes, and that the weights only vary between 0.7 and 1.2.

As a consequence of strong overlap, the standard errors are similar across the columns. Indeed, the efficiency gain of the \ac{EW} scheme relative to the \ac{ATE} based on an efficiency bound comparison using (ref) with $\lambda=\lambda^k$ vs $\lambda=1$ is less than $1.6\%$ for both treatments under homoskedasticity; the gain is even smaller under the \ac{CW} scheme. The reported standard errors, which allow for heteroskedasticity and don't assume known propensity scores, align with this prediction.\footnote{The standard errors reported in parentheses in Panel B are valid for the population analogs $\beta_k$ and $\beta_{\lambda^{\textnormal{CW}}}$, i.e. $E[\lambda^k(W_i)\tau_k(W_i)]/E[\lambda^k(W_i)]$ and $E[\lambda^{\textnormal{CW}}(W_{i})\tau_{k}(W_{i})]/E[\lambda^{\textnormal{CW}}(W_{i})]$. Since these standard errors are potentially conservative when viewed as standard errors for $\beta_k$ and $\beta_{\lambda^{\textnormal{CW}}}$, the standard error comparison gives an upper bound on the cost to estimating the weights.} As discussed in (ref) in (ref), these standard errors are affected by the assumption of known propensity scores, used to derive the weighting schemes underlying the estimates in columns 2 and 3. To gauge the impact of this assumption, we also report a version of the standard errors computed under the assumption that the sample treatment probabilities in each school match the true propensity scores. This changes the standard errors little, showing that there is minimal cost to estimating the weights.

Further Applications

We next study the broader relevance of contamination bias using data from eight additional studies with multiple-treatment regressions. These studies were identified by a systematic search of papers in the AEA Data and Code Repository from 2013--2022 (see (ref) for details). Five studies are experiments like Project STAR\@; the remaining three use observational regressions to estimate racial disparities across multiple race groups (which we interpret as descriptive, following (ref)).\footnote{We focus on observational studies of racial disparities as they often include regressions on multiple minority race “treatments,” use publicly available data, and are easily identifiable by a keyword search.} We replicate a single representative specification for each paper, corresponding to the first relevant regression discussed in the paper's introduction.\footnote{“Relevant” here means a multiple-treatment regression specification with controls, where at least one treatment coefficient was statistically significant. The introduction in cole2013barriers did not discuss any relevant specifications; we instead pick the first specification with variation in treatment probabilities across strata where our results would be most relevant.} (ref) lists the papers and specifications.

We conduct two preliminary analyses of each study before assessing contamination bias and comparing alternative estimators. First, we ensure that the estimation sample satisfies overlap, since otherwise the decomposition in (ref) is typically not identified. If strong overlap fails, we identify a large subset of each analysis sample where it is satisfied. Columns 4 and 5 of (ref) list the number of observations in the full and overlap samples (the sample sizes are equal if the original estimation sample satisfies overlap). Second, we check for propensity score variation in each of the studies. In principle, protocol descriptions can reveal whether some regression controls are necessary (and hence generate propensity score variation) or whether the controls are just added to improve precision. In practice, however, this is not always clear from paper descriptions.\footnote{Moreover, some regression specifications are run on a non-random subsample of the full experimental population (due to, e.g., attrition, or in a susample analysis). This could generate propensity score variation even in simple experimental protocols.} Column 6 of (ref) gives a quantitative sense of the variability in the propensity scores by reporting the standard deviation of the estimated propensity score, showing that its variability in the observational studies is substantially higher; the dagger symbol indicates that a hypothesis test for non-zero variation in the population propensity scores was statistically significant. (ref) details the overlap sample construction and these tests. We replicate the analyses from (ref) for each of the eight papers in (ref); we summarize the takeaways here.

(ref) summarizes the statistical and practical significance of contamination bias in the estimated effect of each treatment for each specification (as estimated in the overlap sample). Column A shows the absolute value of the contamination bias $t$-statistics for each regression coefficient, obtained from the decomposition in (ref). In both columns, we sort treatments within papers by this absolute $t$-statistic and sort papers by the maximum absolute $t$-statistic across treatments. Column B shows a normalized version of the decomposition that divides each term by the standard error of the regression coefficient. The darker bar shows the own-treatment effect component of the decomposition, while the lighter bar denotes the contamination bias component (which can be of the same or opposite sign).

The figure shows economically and statistically meaningful contamination bias in two of the three observational studies while showing no evidence for bias in any of the experimental studies. This aligns with the intuition that the large propensity score variability in observational studies generates much larger variability in the contamination weights. Specifications from both the demel2013 and drexler2014 experiments have some of the smallest contamination bias and also smallest propensity score variation, consistent with the theoretical results that contamination bias requires variation in the contamination weights which in turn requires variation in the propensity scores. On the other hand, the two studies with statistically significant contamination bias (fryerlevitt2013 and weisburst2019) also display the greatest variation in propensity scores. Broadly, these results highlight the importance of testing for contamination bias---especially in observational settings where the included covariates are likely to drive sizable variation in propensity scores and hence contamination weights.

(ref) plots estimates of the treatment effects for each estimator from (ref), again normalizing by the standard error of the regression coefficient. We include a line between the estimates from OLS regression and from the common-weights (CW) estimator we propose. Among observational studies, we see substantial variation across the different estimates and a much larger difference between the OLS estimator and the CW estimator. In the experimental papers, the difference is much smaller.\footnote{The same pattern arises when comparing the estimates in the full sample; see (ref).} This is consistent with the larger propensity score variability in observational studies magnifying the impact of the choice of weighting scheme.

Conclusion

Regressions with multiple treatments and flexible controls are common across a wide range of empirical settings in economics. We show that such regressions generally fail to estimate a convex weighted average of treatment effects: coefficients on each treatment are generally contaminated by the effects of other treatments. We provide intuition for why the influential result of angrist98 fails to generalize to multiple treatments, and show how the contamination bias problem connects to a recent literature studying \ac{DiD} regressions. We then discuss three alternative estimators that are free of this bias.

Our analysis of nine empirical applications finds economically and statistically meaningful contamination bias in observational studies. Contamination bias in experimental studies is more limited, even in papers that display statistically significant variation in the propensity scores. We also find that the choice among alternative estimators that are free of contamination bias matters more in the observational studies. Overall, our analysis highlights the importance of testing the empirical relevance of theoretical concerns with how regression combines heterogeneous effects---particularly in observational studies.

\setstretch{1.1} \printbibliography

\setstretch{1.25}