CI

    confidence interval
OLS

    ordinary least squares
CLT

    central limit theorem
IV

    instrumental variables
ATE

    average treatment effect
RCT

    randomized control trial
VAM

    value-added model
LAN

    locally asymptotically normal
DiD

    difference-in-differences
OVB

    omitted variables bias
FWL

    Frisch-Waugh-Lovell

# Contamination Bias in Linear Regressions††thanks: Contact: paul.goldsmith-
<EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>We thank
Jason Abaluck, Isaiah Andrews, Josh Angrist, Tim Armstrong, Kirill Borusyak,
Kyle Butts, Clément de Chaisemartin, Peng Ding, Jin Hahn, Xavier
D’Haultfœuille, Simon Lee, Bernard Salanié, Tymon Słoczyński, Jonathan Roth,
Jacob Wallace, and numerous seminar participants for helpful comments. Hull
acknowledges support from National Science Foundation Grant SES-2049250.
Kolesár acknowledges support by the Sloan Research Fellowship and by the
National Science Foundation Grant SES-22049356. Mauricio Cáceres Bravo, Jerray
Chang, and Dwaipayan Saha provided expert research assistance. An earlier
draft of this paper circulated under the title “On Estimating Multiple
Treatment Effects with Regression.”

Paul Goldsmith-Pinkham
Yale University  Peter Hull
Brown University  Michal Kolesár
Princeton University

0.1cm0.1cm

We study the interpretation of regressions with multiple treatments and
flexible controls. Such regressions are often used to analyze stratified
randomized control trials with multiple intervention arms, to estimate value-
added (for, e.g., teachers) with observational data, and to leverage the
quasi-random assignment of decision-makers (e.g. bail judges). We show that
these regressions generally fail to estimate convex averages of heterogeneous
treatment effects, even when the treatments are conditionally randomly
assigned and the controls are sufficiently flexible to avoid omitted variables
bias. Instead, estimates of each treatment’s effects are generally
contaminated by a non-convex average of the effects of other treatments. Thus,
recent concerns about heterogeneity-induced bias in regressions leveraging
potential outcome restrictions (e.g. parallel trends assumptions) also arise
with “design-based” identification strategies. We discuss solutions to the
contamination bias and propose a new class of efficient estimators of weighted
average effects that avoid bias. In a re-analysis of the Project STAR trial,
we find minimal bias because treatment effect heterogeneity is largely
idiosyncratic. But sizeable contamination bias arises when effect
heterogeneity becomes correlated with treatment propensity scores.

##  1 Introduction

Consider a linear regression of an outcome $Y_{i}$ on a vector of mutually
exclusive treatment indicators $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 randomized
control trial (RCT), with the randomization protocol varying across some
experimental strata indicators in $W_{i}$. Or, in an education value-added
model (VAM), $X_{i}$ might indicate the matching of students $i$ to different
teachers or schools with $W_{i}$ including measures of student demographics,
lagged achievement, or other controls which give a credible selection-on-
observables assumption. The regression might also be the first stage of an
instrumental variables (IV) regression, perhaps leveraging the as-good-as-
random assignment of multiple decision-makers (e.g. bail judges or benefit
administrators) indicated in $X_{i}$, conditional on some controls $W_{i}$.
These sorts of regressions are widely used across many fields in
economics.111Prominent RCT examples include Project STAR [45] and the RAND
Health Insurance Experiment [50]. Prominent VAM examples include studies of
teachers [39, 19], schools [8, 3, 53], and healthcare institutions [34, 2,
29]. Prominent “judge IV” examples include [43, 49, 24].

This paper shows that such regressions generally fail to identify convex
weighted averages of heterogeneous treatment effects, and discusses solutions
to this problem. The problem may be surprising given an influential result in
[4], showing that regressions on a single binary treatment $D_{i}$ and
flexible controls $W_{i}$ estimate a convex weighted 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. Despite a
set of treatments being completely randomly assigned within groups, as in a
stratified multi-armed RCT, a regression on treatment and strata indicators
generally fails to yield causally interpretable regression coefficients.
Instead, regression estimates of each treatment’s effect are generally
contaminated by a non-convex average of the effects of other treatments: the
regression coefficient for a given RCT treatment arm generally incorporates
the effects of _all_ arms.

We first derive a general characterization of this “contamination bias” in
multiple-treatment regressions. To separate the problem from the well-known
challenge of omitted variables bias (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. We
show that the regression coefficient on each treatment identifies a convex
weighted average of its causal effects, plus a contamination bias term that is
generally non-zero. The bias term is given by a linear combination of the
causal effects of other treatments, with weights that sum to zero. As a
result, each treatment effect estimate will generally incorporate the effects
of other treatments, unless the effects are uncorrelated with the
contamination weights. Since the contamination 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 OVB. To see this intuition clearly, consider the most
flexible specification of controls as a set of strata indicators. OVB is
avoided when the treatments are as good as randomly assigned within strata.
But because the treatments enter the regression linearly, the [4] 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.

We show that contamination bias also arises under an alternative “model-based”
identifying assumption that, instead of making assumptions on the treatment’s
“design” (i.e. propensity scores), assumes the regression parametrization of
covariates is flexible enough to include the conditional mean of the potential
outcome under no treatment, $Y_{i}(0)$. In a linear model with two-way unit
and time fixed effects, this reduces to the parallel trends restriction used
in difference-in-differences (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.222Alternatively $X_{i}$ may indicate multiple contemporaneous
treatments, as in certain “mover” regressions. We show that replacing the
restriction on propensity scores with an assumption on $Y_{i}(0)$ generates an
additional source of bias: the own-treatment effect weights are no longer
guaranteed to be positive. This result shows that the negative weighting and
contamination bias issues documented previously in the context of two-way
fixed effects regressions [[, e.g.,]]goodman2021difference, sun2021estimating,
de2020aer, de2020two, callaway2021difference,
borusyak2021revisiting,wooldridge2021mundlak,hull2018movers are more
general—and conceptually distinct—problems.333Our analysis also relates to
issues with interpreting multiple-treatment IV estimates [11, 41, 42, 36, 47,
13]. Negative weighting arises because regressions leveraging model-based
restrictions on $Y_{i}(0)$ are generally not robust to treatment effect
heterogeneity. Contamination bias arises because linear regression fails to
account for the non-linear dependence across multiple dependent treatments and
controls.

Even more broadly, contamination bias can arise in descriptive regressions
which seek to estimate averages of certain conditional group contrasts without
assuming a causal framework—as in studies of treatment or outcome disparities
across multiple racial or ethnic groups, studies of regional variation in
healthcare utilization or outcomes, or studies of industry wage
gaps.444Prominent examples of such analyses respectively include [28], [60],
and [46]. Our analysis shows that in such regressions the coefficient on a
given group or region averages the conditional contrasts across all other
groups or regions, with non-convex weights.

Our bias characterization also has implications for IV regressions leveraging
multiple correlated instruments, such as indicators for as good as randomly
assigned judges. Contamination bias in the first-stage regression of treatment
on multiple instruments and flexible controls (e.g. courtroom fixed effects)
can generate violations of the effective first-stage monotonicity restriction,
even when conventional first-stage monotonicity is satisfied unconditionally.
We show how this problem is distinct from previous concerns over the
monotonicity assumption in judge IV designs [54, 25, 56, 52] and over
insufficient flexibility in the control parameterization [14].555The
contamination bias issue is also distinct from the [26, 27] critique of
regression to analyze randomized trials, which concerns estimation not
identification.

We then discuss three solutions to the contamination bias problem, and their
trade-offs, in the baseline case of conditionally ignorable treatments. One
conceptually principled solution is to adapt approaches to estimating the
average treatment effect (ATE) of a conditionally ignorable binary treatment
[[, see, e.g.]for a review]imbens2009recent to the multiple treatment case [[,
e.g.]]cattaneo10, ChNeSi21, GrPi22. For example, one could run an expanded
regression that includes interactions between the treatments and (demeaned)
controls.666In the judge IV case, the analogous solution interacts judge
indicators with courtroom fixed effects. Such ATE estimators achieve the
semiparametric efficiency bound under an assumption of strong overlap of the
covariate distribution for units in each treatment arm. But this approach may
be infeasible or yield imprecise estimates under limited overlap—a common
scenario in practice [21].

This practical consideration motivates an alternative approach: estimating a
convex weigh-ed average of treatment effects, as regression does in the binary
treatment case, while avoiding the contamination bias of multiple-treatment
regressions. We derive the weights that are “easiest” to estimate, in that
they minimize a semiparametric efficiency bound under homoskedasticity. These
optimal weights coincide with the implicit linear regression weights when the
treatment is binary (i.e. the [4] case), formalizing a virtue of regression
adjustment. In the multiple treatment case, the optimal weights for a given
treatment-control contrast are similarly given by a linear regression which
restricts estimation to the individuals who are either in the control group or
the treatment group of interest. For estimating effects that are directly
comparable across multiple treatments, our optimal weight characterization
leads to a new estimator. We give guidance for how applied researchers can
gauge the extent of contamination bias in practice and apply the different
solutions with help from a new Stata package, multe.777The Stata package is
available at https://github.com/gphk-metrics/stata-multe.

We illustrate the contamination bias problem and solutions in an application
to the Project STAR trial, which randomized students within schools to either
a small classroom treatment, a teaching aide treatment, or control conditions.
We find the potential for sizeable bias in estimates of both treatment
effects, from regressions with school fixed effects, due to significant
treatment effect heterogeneity. Nevertheless, we show that the actual
contamination is likely to be minimal, as the effect heterogeneity turns out
to be largely uncorrelated with the contamination weights. The application
thus highlights the importance of testing the empirical relevance of
theoretical concerns with how regression combines heterogeneous effects.

We structure the rest of the paper as follows. Section 2 illustrates
contamination bias in a simple example with two mutually exclusive treatment
indicators and one binary control. Section 3 characterizes the general problem
in regressions with multiple treatments and flexible controls, and discusses
connections to previous analyses. Section 4 discusses the robustness and
efficiency properties of three solutions. Section 5 gives guidance for
measuring and avoiding contamination bias in practice, and illustrates these
tools in the Project STAR experiment. Section 6 concludes. All proofs and
extensions are given in Appendix A.

##  2 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 in our application in Section 5.2.

###  2.1 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 a
constant:

$Y_{i}=\alpha+\beta D_{i}+\gamma W_{i}+U_{i}.$ (1)

By definition, $U_{i}$ is a mean-zero regression residual that is uncorrelated
with $D_{i}$ and $W_{i}$. [45], for example, primarily studied the effect of
small class size $D_{i}$ on the test scores $Y_{i}$ of middle school students
indexed by $i$. Project STAR randomized students to classes within schools
with at least three classes per grade. The number of students assigned to each
intervention thus varied both by the number of students in a school and the
relative classroom size. To account for this non-random treatment variation,
[45] followed earlier analyses of Project STAR in estimating regressions with
school (and sometimes school-by-period) fixed effects as controls. Such
specifications are often found in stratified RCTs 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 eq. 1 corresponds to a
stylized two-school version of a Project STAR regression.

We wish to interpret the regression 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_{i}=Y_{i}(1)-Y_{i}(0)$, and we can write realized achievement as
$Y_{i}=Y_{i}(0)+\tau_{i}D_{i}$. To formalize the random assignment of
treatment within schools, we assume that $D_{i}$ is conditionally independent
of potential outcomes given the control $W_{i}$:

$\left(Y_{i}(0),Y_{i}(1)\right)\perp D_{i}\mid W_{i}.$ (2)

[4] showed that regression coefficients like $\beta$ identify a weighted
average of within-strata ATEs, with convex weights.888See [61] for an
alternative representation of this estimand, in terms of conditional average
effects on the treated and untreated, under slightly different assumptions. In
our stylized Project STAR regression, this result shows:

$\beta=\phi\tau(0)+(1-\phi)\tau(1),\quad\text{where}\quad\phi=\frac{\operatorname{var}(D_{i}\mid
W_{i}=0)\Pr(W_{i}=0)}{\sum_{w=0}^{1}\operatorname{var}(D_{i}\mid
W_{i}=w)\Pr(W_{i}=w)}\in[0,1]$ (3)

gives a convex weighting scheme and $\tau(w)=E[Y_{i}(1)-Y_{i}(0)\mid W_{i}=w]$
is the 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(w)$ across the two schools.

Equation 3 can be derived by applying the Frisch-Waugh-Lovell (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:

$\beta=\frac{E[\tilde{D}_{i}Y_{i}]}{E[\tilde{D}_{i}^{2}]}=\frac{E[\tilde{D}_{i}Y_{i}(0)]}{E[\tilde{D}_{i}^{2}]}+\frac{E[\tilde{D}_{i}D_{i}\tau_{i}]}{E[\tilde{D}_{i}^{2}]},$
(4)

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,

$E[\tilde{D}_{i}Y_{i}(0)]=E[E[\tilde{D}_{i}Y_{i}(0)\mid
W_{i}]]=E[E[\tilde{D}_{i}\mid W_{i}]E[Y_{i}(0)\mid W_{i}]]=0.$ (5)

The first equality in eq. 5 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 eq. 4 is zero. Analogous arguments show that

$E[\tilde{D}_{i}D_{i}\tau_{i}]=E[E[\tilde{D}_{i}D_{i}\tau_{i}\mid
W_{i}]]=E[E[\tilde{D}_{i}D_{i}\mid W_{i}]E[\tau_{i}\mid
W_{i}]]=E[\operatorname{var}(D_{i}\mid W_{i})\tau(W_{i})],$

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 eq. 4 as

$\beta=\frac{E[\operatorname{var}(D_{i}\mid
W_{i})\tau(W_{i})]}{E[\operatorname{var}(D_{i}\mid
W_{i})]}=\phi\tau(0)+(1-\phi)\tau(1),$

proving the representation of $\beta$ in eq. 3.

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 FWL theorem, treatment
coefficients like $\beta$ can always be represented as in eq. 4 even without
this property. We next show, however, that the remaining steps in the
derivation of eq. 3 fail when an additional treatment arm is included. This
failure can be attributed to the fact that the auxiliary FWL regression
delivers a treatment residual that is uncorrelated, but not mean-independent
of the other regressors, leading to an additional bias term in the expression
for the regression coefficient.

###  2.2 Contamination Bias with Two Randomized Treatments

In reality, as noted above, Project STAR randomized two mutually exclusive
interventions within schools: a reduction in class size ($D_{i}=1$) and the
introduction of 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})^{\prime}$, where the first element
$X_{i1}=\operatorname{\mathbbm{1}}\\{D_{i}=1\\}$ indicates assignment to a
small class and the second element
$X_{i2}=\operatorname{\mathbbm{1}}\\{D_{i}=2\\}$ indicates assignment to a
class with a full-time aide. We continue to include a constant and the school
indicator $W_{i}$ as controls, yielding the regression

$Y_{i}=\alpha+\beta_{1}X_{i1}+\beta_{2}X_{i2}+\gamma W_{i}+U_{i}.$ (6)

To account for the second treatment, 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 FWL theorem to write

$\beta_{1}=\frac{E[\accentset{\approx}{X}_{i1}Y_{i}]}{E[\accentset{\approx}{X}_{i1}^{2}]}=\frac{E[\accentset{\approx}{X}_{i1}Y_{i}(0)]}{E[\accentset{\approx}{X}_{i1}^{2}]}+\frac{E[\accentset{\approx}{X}_{i1}X_{i1}\tau_{i1}]}{E[\accentset{\approx}{X}_{i1}^{2}]}+\frac{E[\accentset{\approx}{X}_{i1}X_{i2}\tau_{i2}]}{E[\accentset{\approx}{X}_{i1}^{2}]},$
(7)

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 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 eq. 7 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 eq. 5, since the auxilliary
regression residuals are still uncorrelated with any individual characteristic
like $Y_{i}(0)$.999To 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 eq. 5, are both uncorrelated with
$Y_{i}(0)$. It follows that $E[\accentset{\approx}{X}_{i1}Y_{i}(0)]=0$. In
this sense, the regression does not suffer from OVB. However, we do not
generally have $E[\accentset{\approx}{X}_{i1}X_{i2}\tau_{i2}]=0$. Instead,
simplifying eq. 7 by the same steps as before leads to the expression

$\beta_{1}=E[\lambda_{11}(W_{i})\tau_{1}(W_{i})]+E[\lambda_{12}(W_{i})\tau_{2}(W_{i})]$
(8)

as a generalization of eq. 3. 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 eq. 3. Thus, if not for the
second term in eq. 8, $\beta_{1}$ would similarly identify a convex average of
the conditional ATEs $\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 _contamination bias_ term in eq. 8 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.

###  2.3 Illustration and Intuition

A simple numerical example helps make the contamination bias problem concrete.
Suppose, in the previous setting, 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 a classroom with a full-time aide 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, such that 45 percent of students were assigned to a small classroom,
45 percent were assigned to a classroom with a full-time aide, and only 10
percent were assigned to the control group. Therefore, $p_{1}(0)=0.05$,
$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 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$, while
in school $1$ it is $-0.8$.101010Here the conditional correlation is
$\operatorname{corr}(\tilde{X}_{i1},\tilde{X}_{i12}\mid
W_{i})=-\sqrt{p_{1}(W_{i})/(1-p_{1}(W_{i}))}\sqrt{p_{2}(W_{i})/(1-p_{2}(W_{i}))}$.
The overall regression of $\tilde{X}_{i1}$ on $\tilde{X}_{i2}$ will average
over these two correlations, leading to a misspecified residual
$\accentset{\approx}{X}_{i1}$ that is correlated with $X_{i2}$ within each
school. Figure 1 illustrates this averaging by plotting the different
potential pairs of the two demeaned treatments ($\tilde{X}_{i1}$,
$\tilde{X}_{i2}$), with the two school strata in different colors and shapes.
The figure shows how within the first school, the value of the demeaned class
aide treatment is only weakly predictive of the small classroom treatment,
while it is highly predictive in the second school. The overall regression
line in black averages over these different relationships, leading to
residuals that are predictive of the value of the class aide treatment.

Figure 1: Regression of Small Classroom Treatment on Class Aide Treatment

_Notes:_ This figure plots values of the demeaned class aide treatment
($\tilde{X}_{2i}$, the x-axis) against values of the demeaned small classroom
treatment ($\tilde{X}_{1i}$, the y-axis) in our numerical example. The size of
the points corresponds to the density of observations. The solid red and blue
lines mark the within-school regression of the two residualized treatments,
while the dashed black line is the overall regression line. The residuals from
this line give $\accentset{\approx}{X}_{i1}$.

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 the 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 the school $0$, the aide has no
effect, $\tau_{2}(0)=0$. Equation 8 then shows that the regression coefficient
on the first treatment identifies

$\beta_{1}=E[\lambda_{11}(W_{i})\cdot
0]+E[\lambda_{12}(W_{i})\tau_{2}(W_{i})]=0+(-99/106\times 1+99/106\times
0)/2\approx-0.47.$

Thus, in this example, a researcher would conclude that small classrooms have
a sizeable 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 made arbitrarily large or small (and positive or negative),
depending on the heterogeneity of the teaching aide effects across schools.

To build further intuition for eq. 8, it is useful to consider two cases where
the contamination bias term is zero. First, suppose the average effects of the
teaching aide treatment are constant across the two schools:
$\tau_{2}(0)=\tau_{2}(1)\equiv\tau_{2}$. 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$.
Thus, the contamination weights on the second treatment effects average to
zero, and the contamination bias disappears:
$E[\lambda_{12}(W_{i})\tau_{2}(W_{i})]=E[\lambda_{12}(W_{i})]\tau_{2}=0$. More
generally, the contamination bias will be small when the variation in average
teacher’s aide treatment effects across schools $\tau_{2}(W_{i})$ is small, or
when this treatment effect heterogeneity is only weakly correlated with the
contamination weights across schools.

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 again coincide with $X_{i1}-E[X_{i1}\mid W_{i}]$. The
coefficient on $X_{i1}$ in eq. 6 can therefore be shown to be equivalent to
the previous eq. 3, identifying the same convex average of $\tau_{1}(w)$. This
case highlights that dependence across treatments is necessary for the
contamination bias to arise.

Before proceeding to a general characterization of contamination bias, we note
that the above intuition about the non-linear conditional expectation
$E[X_{i1}\mid W_{i},X_{i2}]$ also suggests a simple solution to the problem.
By including interactions of $W_{i}$ and $X_{i2}$ in eq. 6, the auxiliary
regression of $X_{i1}$ on the other regressors will be saturated and thus will
capture the inherent nonlinearity in $E[X_{i1}\mid W_{i},X_{i2}]$. We show
below how such interacted regressions can obviate contamination bias. In
particular, we show how a particular interacted regression specification gives
an efficient estimator of (unweighted) ATEs that is immune to the bias of the
simpler specification. We then propose a new class of estimators which—as with
the [4] result for binary treatments—identify a convex average of conditional
ATEs. These estimators may yield smaller standard errors, while still being
free from bias.

##  3 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}=\operatorname{\mathbbm{1}}\\{D_{i}=k\\}$ for values of an underlying
treatment $D_{i}\in\\{0,\dotsc,K\\}$ (with the
$\operatorname{\mathbbm{1}}\\{D_{i}=0\\}$ indicator omitted). We extend the
characterization to a general treatment vector in Section A.1.

We suppose the effects of $X_{i}$ on $Y_{i}$ are estimated by a partially
linear model:

$Y_{i}=X_{i}^{\prime}\beta+g(W_{i})+U_{i},$ (9)

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

$(\beta,g)=\operatorname*{argmin}_{\tilde{\beta}\in\mathbb{R}^{K},\tilde{g}\in\mathcal{G}}E[(Y_{i}-{X}_{i}^{\prime}\tilde{\beta}-\tilde{g}(W_{i}))^{2}]$
(10)

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 eq. 9 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 [[, e.g.]]robinson88.

Two examples highlight the generality of this setup and are useful for
developing our characterization of contamination bias below:

######  Example 1 (_Multi-Armed RCT_).

$W_{i}$ is a vector of mutually-exclusive indicators for experimental strata,
within which $X_{i}$ is randomly assigned to individuals $i$. $g$ is linear.

######  Example 2 (_Two-Way Fixed Effects_).

$i=(j,t)$ indexes panel data, with units $j$ observed over a set of periods
$t$, and $W_{i}$ contains unit and time indicators. $X_{i}$ contains
indicators for leads and lags relative to a unit’s deterministic adoption date
$A(j)$, $\operatorname{\mathbbm{1}}\\{A(j)=t-r\\}$, with the indicator
$\operatorname{\mathbbm{1}}\\{A(j)=-1\\}$ excluded to normalize the treatment
effects. $g$ is linear. The expectation in eq. 10 is over the set of units and
time periods in the population from which $i=(j,t)$ is drawn.

Example 1 nests the motivating RCT example in Section 2, allowing for an
arbitrary number of experimental strata in $W_{i}$ and random treatment arms
in $X_{i}$. Example 2 shows that our setup also nests the kind of regressions
considered in an existing literature on various two-way fixed effect and
difference-in-differences regressions [[, e.g.]]goodman2021difference,
sun2021estimating, de2020aer,de2020two,callaway2021difference,
borusyak2021revisiting,wooldridge2021mundlak.111111Some papers in this
literature study issues we do not consider, such as when researchers fail to
include indicators for all relevant treatment states. This specification of
$X_{i}$ will generally add bias terms to our decomposition of $\beta$, below.
Similarly, we do not consider multicollinearity issues like in [15], by
implicitly assuming a unique solution to eq. 10. In Example 2 this means we
assume some units are “never treated,” with $A(j)=\infty$. See [59] for a
recent review of these and other issues. Here $X_{i}$ is a non-random function
of the unit and time indicators in $W_{i}$. Instead of (or in addition to)
indicators for the leads and lags of each unit’s adoption date $A(j,t)$ we
might imagine $X_{i}$ including multiple static treatment indicators, as in
“mover regressions” that leverage over-time transitions [35].

As a first step towards characterizing the $\beta$ treatment coefficient
vector, we solve the minimization problem in eq. 10. 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 [[, e.g.]Theorem 11.1]vdv98 that

$\beta=E[\tilde{{X}}_{i}\tilde{{X}}_{i}^{\prime}]^{-1}E[\tilde{X}_{i}Y_{i}].$
(11)

A further application of the FWL theorem allows us to write each treatment
coefficient as

$\beta_{k}=\frac{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})^{\prime}$.

###  3.1 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 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 Example 2, since $X_{i}$ is a non-random function of
$W_{i}$. We return to practical questions of propensity score support in
Section 4.

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:

######  Assumption 1.

$E[Y_{i}(k)\mid D_{i},W_{i}]=E[Y_{i}(k)\mid W_{i}]$ for all $k$.

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:

$\left(Y_{i}(0),\dotsc,Y_{i}(K)\right)\perp{D}_{i}\mid W_{i}.$ (12)

Such conditional random assignment appears in Example 1. In Example 2, where
treatment is a non-random function of the two-way fixed effects in $W_{i}$,
Assumption 1 holds trivially.

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

######  Assumption 2.

Let $\mu_{0}(w)=E[Y_{i}(0)\mid W_{i}=w]$ and recall $p_{k}(w)=E[X_{ik}\mid
W_{i}=w]$. Either

$p_{k}\in\mathcal{G}$ (13)

for all $k$, or

$\mu_{0}\in\mathcal{G}.$ (14)

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$, eq. 13 requires the propensity scores to be linear in
$W_{i}$ [[, cf. eq. (30) in]]AnKr99handbook. This condition holds trivially in
Example 1, 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 are given by $\tilde{X}_{i}=X_{i}-p(W_{i})$.

In Example 2, with $X_{i}$ being a deterministic function of unit and time
indices and a linear specification for $g$, eq. 13 fails because the
propensity scores are binary—they cannot be captured by a linear function of
the unit and time indicators in $W_{i}$. However, eq. 14 can still be
satisfied by a parallel trends assumption: that the average untreated
potential outcomes $Y_{i}(0)$ are linear in the unit and time
effects.121212Identification based on eq. 13 can be seen as “design-based” in
that it leverages only the conditional random assignment of $D_{i}$ and
specifies the treatment assignment process. Identification based on eq. 14 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 Assumption 2, the specification of controls is
flexible enough to avoid 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 definition of
$\beta$ in eq. 11 that there is no such OVB when Assumption 2 holds; the
coefficient vector identifies the constant effects:

$\beta=E[\tilde{{X}}_{i}\tilde{{X}}_{i}^{\prime}]^{-1}E[\tilde{{X}}_{i}Y_{i}]=E[\tilde{{X}}_{i}\tilde{{X}}_{i}^{\prime}]^{-1}({E[\tilde{{X}}_{i}Y_{i}(0)]}+E[\tilde{{X}}_{i}\tilde{X}_{i}^{\prime}]\tau)\\\
=E[\tilde{{X}}_{i}\tilde{{X}}_{i}^{\prime}]^{-1}\underbrace{E[\tilde{X}_{i}E[Y_{i}(0)\mid
W_{i}]]}_{=0}+\tau=\tau.$

Here the first line 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
line uses the law of iterated expectations and Assumption 1. Under eq. 13,
$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 eq. 14 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, 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
robustness property have been previously observed in, for instance, [57].

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 we still have
$E[\tilde{X}_{i}Y_{i}(0)]=0$ under Assumptions 1 and 2, it follows from eq. 11
that

$\beta=\frac{E[\tilde{X}_{i}X_{i}\tau_{i}]}{E[\tilde{X}_{i}^{2}]}=E[\lambda_{11}(W_{i})\tau(W_{i})],\quad\text{with}\quad\lambda_{11}(W_{i})=\frac{E[\tilde{X}_{i}X_{i}\mid
W_{i}]}{E[\tilde{X}_{i}X_{i}]},$ (15)

where the second equality uses iterated expectations and the fact that
$E[\tilde{{X}}_{i}^{2}]=E[\tilde{X}_{i}{X}_{i}]$. Under eq. 13,
$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})]}\geq 0$. This extends the [4]
result to a general control specification; versions of this extension appear
in, for instance, [6], [7, Chapter 3.3], and [9]. The result provides a
rationale for estimating the effect of a scalar as good as randomly assigned
treatment using a partially linear model: so long as the specification of
$\mathcal{G}$ is rich enough so that eq. 13 holds, this model will identify a
convex average of heterogeneous treatment effects. Moreover, as we will show
in Section 4, the weights $\lambda_{11}(W_{i})$ are efficient in that they
minimize the semiparametric efficiency bound (conditional on the controls) for
estimating a weighted-average treatment effect. This result makes the
partially linear specification (9) especially appealing with a single binary
treatment. On the other hand, when eq. 14 holds but eq. 13 does not, the
weights $\lambda_{11}(W_{i})$ need not be positive. We return to this point in
Section 3.2.

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

######  Proposition 1.

Under Assumptions 1 and 2, the treatment coefficients in the partially linear
model (9) identify

$\beta_{k}=E[\lambda_{kk}(W_{i})\tau_{k}(W_{i})]+\sum_{\ell\neq
k}E[\lambda_{k\ell}(W_{i})\tau_{\ell}(W_{i})],$ (16)

where

$\lambda_{kk}(W_{i})=\frac{E[\accentset{\approx}{X}_{ik}X_{ik}\mid
W_{i}]}{E[\accentset{\approx}{X}_{ik}^{2}]}\qquad\text{and}\qquad\lambda_{k\ell}(W_{i})=\frac{E[\accentset{\approx}{X}_{ik}X_{i\ell}\mid
W_{i}]}{E[\accentset{\approx}{X}_{ik}^{2}]}.$

These weights satisfy $E[\lambda_{kk}(W_{i})]=1$ and
$E[\lambda_{k\ell}(W_{i})]=0$. Furthermore, if eq. 13 holds,
$\lambda_{kk}(W_{i})\geq 0$ for each $k$.

Proposition 1 shows that the coefficient on $X_{ik}$ in eq. 9 is a sum of two
terms. The first term is a weighted average of conditional ATEs
$\tau_{k}(W_{i})$, with weights $\lambda_{kk}(W_{i})$ that average to one and
are guaranteed to be convex when eq. 13 holds. This term generalizes the
characterization of the single-treatment case, eq. 15. The second term is a
weighted average of treatment effects for _other_ treatments
$\tau_{\ell}(W_{i})$, with weights $\lambda_{k\ell}(W_{i})$ that average to
zero. Because these contamination weights are zero on average, they must be
negative for some values of the controls unless they are all identically zero.

Each treatment coefficient ${\beta}_{k}$ thus generally suffers from
contamination bias. Two exceptions are when $\lambda_{k\ell}(W_{i})=0$ almost-
surely for all $\ell\neq k$, and when the conditional effects of these other
treatments are homogeneous such that $\tau_{\ell}(W_{i})=\tau_{\ell}$. In the
second case
$E[\lambda_{k\ell}(W_{i})\tau_{\ell}(W_{i})]=\tau_{\ell}E[\lambda_{k\ell}(W_{i})]=0$,
so there is no contamination bias term. By the law of iterated expectations
the first case holds if $E[\accentset{\approx}{X}_{ik}\mid X_{i,-k},W_{i}]=0$,
or, equivalently, if the conditional expectation of $X_{ik}$ given $X_{i,-k}$
and $W_{i}$ is partially linear (i.e.
$E[X_{ik}\mid{X}_{i,-k},W_{i}]=X_{i,-k}^{\prime}\alpha+g_{k}(W_{i})$ for some
vector $\alpha$ and $g_{k}\in\mathcal{G}$). In other words, it holds when the
assignment of treatment $k$ depends linearly on the other treatment indicators
and a flexible function of the controls. This condition is the analog of
condition (13) if we interpret $X_{ik}$ as a binary treatment of interest, and
$X_{i,-k}^{\prime}\alpha+g_{k}(W_{i})$ as a specification for the controls.
However, with mutually exclusive treatments, it cannot hold 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}$, we have $\alpha_{\ell}=-g_{k}(W_{i})$ for all elements
$\alpha_{\ell}$ of $\alpha$. This in turn implies the assignment of treatment
$k$ doesn’t depend on $W_{i}$, which can’t be the case unless the propensity
score $p_{k}(W_{i})$ is constant.

A third weaker case of no contamination bias is when the weights
$\lambda_{k\ell}(W_{i})$ and conditional ATEs $\tau_{\ell}(W_{i})$ vary, but
are uncorrelated with each other. More generally, contamination bias may be
minimal when the contamination weights $\lambda_{k\ell}(W_{i})$ and the
conditional ATEs $\tau_{\ell}(W_{i})$ are weakly correlated: that is, when the
factors influencing treatment effect heterogeneity are largely unrelated to
the factors influencing the treatment assignment process. We return to this
possibility in our empirical application (Section 5.2).

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

######  Remark 1.

Since the weights in eq. 16 are functions of the variances
$E[\accentset{\approx}{X}_{ik}^{2}]$ and covariances
$E[\accentset{\approx}{X}_{ik}X_{i\ell}]$ and
$E[\accentset{\approx}{X}_{ik}X_{ik}]$, they are identified and can be used to
further characterize each $\beta_{k}$ coefficient. For example, the
contamination bias term can be bounded by the identified contamination weights
$\lambda_{k\ell}(W_{i})$ and bounds on the heterogeneity in conditional ATEs
$\tau_{\ell}(W_{i})$. We illustrate such an approach in our empirical
application.

######  Remark 2.

The results in Proposition 1 are stated for the case when $X_{i}$ are mutually
exclusive treatment indicators. In Section A.1 we relax this assumption to
allow for combinations of non-mutually exclusive treatments (either discrete
or continuous). In this case, the own-treatment weights $\lambda_{kk}(W_{i})$
may be negative even if eq. 13 holds.

######  Remark 3.

While we derived Proposition 1 in the context of a causal model, analogous
results follow for descriptive regressions that do not assume the existence of
potential outcomes or impose Assumption 1. Consider, specifically, the goal of
estimating a convex average of conditional group contrasts $E[Y_{i}\mid
D_{i}=k,W_{i}=w]-E[Y_{i}\mid D_{i}=0,W_{i}=w]$ with a partially linear model
eq. 9 and replace condition (14) with an assumption that $E[Y_{i}\mid
D_{i}=0,W_{i}=w]\in\mathcal{G}$. The same steps that lead to Proposition 1
then show that such regressions also generally suffer from contamination bias:
the coefficient on a given group indicator averages the conditional contrasts
across all other groups, with non-convex weights. Furthermore, the weights on
own-group conditional contrasts are not necessarily positive. These sorts of
conditional contrast comparisons are therefore not generally robust to
misspecification of the conditional mean, $E[Y_{i}\mid D_{i},W_{i}]$.

###  3.2 Implications

Proposition 1 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
potential for negative weighting of a treatment’s _own_ effects when condition
(14) holds but condition (13) fails. This negative weighting issue is relevant
in various DiD regressions and related estimators which rely on such models
for untreated potential outcomes (via parallel trends assumptions) while
conditioning on treatment assignment. Although the recent DiD literature
focuses on two-way fixed effect regressions, Proposition 1 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 eq. 13 holds.131313More broadly, negative weighting issues arise whenever
the covariate specification is not flexible enough for eq. 13 to hold. For
example, consider a scalar covariate $W_{i}$ with a uniform distribution on
$[0,1]$, and a binary treatment with non-linear propensity score
$p(W_{i})=\min(2W_{i},0.9)$. Suppose that $Y_{i}(0)=0$, so that (14) holds.
Then $\lambda_{11}(W_{i})$ is negative for $W_{i}\geq 2911/3402\approx 0.855$.
If, say $Y_{i}(1)=W_{i}$, so that $\tau(W_{i})=W_{i}\geq 0$, the regression
coefficient on the treatment is negative despite uniformly non-negative
treatment effects.

Second, in the multiple treatment case, there is a potential for contamination
bias from _other_ treatment effects regardless of which condition in
Assumption 2 holds. This form of bias is thus relevant whenever one uses an
additive covariate adjustment, regardless of how flexibly the covariates are
specified. Versions of this problem have been noted in, for example, the [62]
analysis of DiD regressions with treatment leads and lags or the [35] analysis
of mover regressions.141414The negative weights issue raised in [22] (when
$K=1$), and the related issue that own-treatment weights may be negative in
[62] and [23] (when $K>1$), arise because the treatment probability is not
linear in the unit and time effects. If eq. 13 holds with $K=1$, Proposition 1
shows $\beta$ estimates a convex combination of treatment effects. This covers
the setting considered in Theorem 1(iv) in [10]. In their Comment 2, [10] say
that “the sum of the weights [used in Theorem 1(iv)] is one, although some of
the weights may be negative.” Proposition 1 shows these weights are, in fact,
non-negative.

The characterization in Proposition 1 also relates to concerns in interpreting
multiple-treatment IV estimates with heterogeneous treatment effects [11, 41,
42, 36, 47, 13]. This connection comes from viewing eq. 9 as the second stage
of an IV model estimated by a control function approach; in the linear 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 [37]. But as
in Proposition 1 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, Proposition 1 has implications for single-treatment IV estimation
with multiple instruments and flexible controls. The first stage of such IV
regressions will tend to have the form of eq. 9, where now $Y_{i}$ is
interpreted as the treatment and $X_{i}$ gives the vector of instruments.
Proposition 1 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 the true effects are always positive. In other words, the cross-
instrument contamination of causal effects may cause monotonicity violations,
even when specifications with individual instruments would be appropriate.
This issue is distinct from previous concerns over monotonicity failures in
multiple-instrument designs [54, 25, 56, 52], which are generally also present
in such just-identified specifications. It is also distinct from some concerns
about insufficient flexibility in the control specification when monotonicity
holds unconditionally [14].

This new monotonicity concern may be especially important in “judge” 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 IV estimators can
be thought of as implementing versions of a jackknife IV estimator [5], based
on a first stage that uses examiner indicators as instruments, similar to eq.
9. Proposition 1 thus raises a new concern with these IV analyses when the
controls (such as courtroom or time fixed effects) are needed to ensure
ignorable treatment assignment.151515As we discuss in Section 4, one solution
to this problem is to interact the examiner instruments with the controls,
which would amount to computing “leniency” separately within location and time
cells. This may greatly increase the effective number of instruments,
heightening concerns of many-instrument bias in finite samples as well as the
importance of appropriate leave-one-out corrections [[,
e.g.,]]kolesar2013estimation.

##  4 Solutions

We now discuss solutions to the contamination bias problem raised by
Proposition 1. We focus on the case of conditionally ignorable treatment
assignment (in the sense that eq. 12 holds, and the propensity scores are not
degenerate) since solutions that allow for degenerate propensity scores are
generally different and have been previously explored in the literature in the
context of DiD regressions. We refer readers to [23, 62, 16, 15, 67] for such
solutions.

We propose three solutions, each identifying a distinct causal parameter.
First, in Section 4.1, we discuss estimation of ATEs. The other two solutions,
discussed in Sections 4.2 and 4.3, estimate weighted averages of individual
treatment effects using weights that are easiest to estimate in that they
minimize the semiparametric efficiency bound for estimating weighted ATEs
under homoskedasticity. If the weights are allowed to vary across treatments,
it is optimal to estimate the effect of each $k$ using the partially linear
model in eq. 9, but in a sample restricted to individuals in the control group
and to those receiving treatment $k$. If the weights are constrained to be
common across treatments, this leads to a new weighted regression estimator.

###  4.1 Estimating Average Treatment Effects

Many estimators exist for the ATE of binary treatments (see [38] for a
review). A number of these approaches extend naturally to multiple treatments,
including matching, inverse propensity score weighting, regressions with
interactions, or doubly-robust combinations of these methods (see, among
others, [17], [18], and [31]).

Rather than reviewing all of these approaches, we briefly outline a simple
implementation of one method which follows the intuition given at the end of
Section 2. Namely, one may estimate the ATE vector $\tau$ by expanding the
partially linear model in eq. 9 to include treatment interaction terms. This
generalizes the implementation in the binary treatment case discussed in [38,
Section 5.3]. Consider the expanded model

$Y_{i}=X_{i}^{\prime}\beta+q_{0}(W_{i})+\sum_{k=1}^{K}X_{ik}\left(q_{k}(W_{i})-E[q_{k}(W_{i})]\right)+\dot{U}_{i},$
(17)

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, eq. 17 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)$. When Assumption 1
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 ATEs.

Following the intuition at the end of Section 2, the added interactions in eq.
17 ensure that each treatment coefficients $\beta_{k}$ is determined only by
the outcomes in treatment arms with $D_{i}=0$ and $D_{i}=k$, avoiding the
other-treatment contamination bias in Proposition 1. Demeaning the
$q_{k}(W_{i})$ in the interactions ensures they are appropriately centered to
interpret the coefficients on the uninteracted $X_{ik}$ as ATEs.

Estimation of eq. 17 by least squares is conceptually straightforward, with
sample averages replacing expectations. Furthermore, it can be shown that the
resulting estimator achieves the semiparametric efficiency bound under strong
overlap (i.e. when the propensity score is bounded away from zero and one)
when implemented as a series estimator: it is impossible to construct another
regular estimator of the ATE with smaller asymptotic variance.

Nonetheless, under weak overlap, the estimator may be imprecise, with poor
finite-sample behavior. This is not a shortcoming of the specific estimator:
[40] shows that identification of the ATE is irregular under weak overlap, and
it is not possible to estimate it at a $\sqrt{N}$-rate. These results
formalize the intuition that it is difficult to reliably estimate the
counterfactual outcomes for observations with extreme propensity scores.
Overlap concerns tend to be more severe with multiple treatments, because some
propensity scores necessarily become closer to zero or one as more treatment
arms are added. We next turn to the problem of estimating weighted averages of
conditional ATEs that downweight these difficult-to-estimate counterfactuals.

###  4.2 Efficient Weighted 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.161616In a slight abuse of notation relative to Section 3, 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 to
estimating the ATE based on eq. 17. 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})\sum_{k=0}^{K}\tau_{k}(W_{i})/\sum_{i=1}^{N}\lambda(W_{i})$,
the weighted ATE of treatment $k$ across different strata. Second, we specify
$c$ so as to allow us to simultaneously contrast the effects of all $K$
treatments—we discuss this further below.

Given the contrast vector $c$, we consider the problem of finding the
weighting scheme $\lambda(W_{i})$ that is the “easiest” to estimate in that it
leads to the smallest possible standard errors. This objective has two
motivations. First is a robustness concern: a researcher would like to
estimate a given contrast as efficiently as possible, at least under the
benchmark of constant treatment effects, while being robust to the possibility
that the effects are heterogeneous. Under constant effects, the weighting
$\lambda(W_{i})$ is immaterial. But the robustness property ensures that the
estimand retains a causal interpretation as a convex average of conditional
contrasts under weak conditions, avoiding the contamination bias displayed by
the regression estimator per Proposition 1. Such a motivation presumably
underlies the popularity of regression for estimating the effect of a binary
treatment: the regression estimator is efficient under homoskedasticity and
constant treatment effects, while, by the [4] result, retaining a causal
interpretation under heterogeneous effects.

The second motivation is that the easiest-to-estimate weighting scheme gives a
bound on the information available in the data: if these weights nonetheless
yield overly large standard errors, inference on other treatment effects (such
as the unweighted ATE) will be at least as uninformative. Computing standard
errors for this efficient weighted average of treatment effects can be useful
as it reveals whether informative conclusions—regardless of how one specifies
the treatment effect of interest—are only possible under additional
assumptions or with the aid of additional data. If the easiest-to-estimate
weighting scheme yields small standard errors even though the standard errors
for the unweighted ATE are large, it can be concluded that the data is
informative about _some_ treatment effects—even if it is not informative about
the unweighted average.

In fact, our solution below shows that in the binary treatment case the
easiest-to-estimate weighting scheme is exactly the same as the weights used
by regression. This special case illustrates the second motivation: the
optimal binary 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}))$, which tend to zero as $p_{1}(W_{i})$
tends to zero or one. Regression thus downweights observations with extreme
propensity scores where the estimation of counterfactual outcomes is
difficult, avoiding the poor finite-sample behavior of ATE estimators under
weak overlap and allowing regression to be informative even in cases when it
is not possible to precisely estimate the unweighted ATE. More generally,
since regression solves the efficient binary treatment weighting scheme,
regression estimates establish the extent to which internally valid and
informative inference for _any_ causal effect are possible with the data at
hand.

We derive the easiest-to-estimate weighting scheme for multiple treatments in
two steps. First, we establish an efficiency 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 weighted average is the “easiest” to estimate, in
that it minimizes the semiparametric efficiency bound over the choice of
$\lambda(W_{i})$. When the contrast vector is specified to allow simultaneous
comparison of all treatments, estimation of this efficient weighted average
leads to a new estimator; we discuss its implementation when the propensity
score is not known in Section 4.3.

The following proposition establishes the first step of our derivation:

######  Proposition 2.

Suppose eq. 12 holds in an i.i.d. sample of size $N$, with known non-
degenerate propensity scores $p_{k}(W_{i})$. Let
$\sigma^{2}_{k}(W_{i})=\operatorname{var}(Y_{i}(k)\mid W_{i})$. Consider the
problem of estimating the weighted average of contrasts

$\theta_{\lambda,c}=\frac{1}{\sum_{i=1}^{N}\lambda(W_{i})}\sum_{i=1}^{N}\lambda(W_{i})\sum_{k=0}^{K}c_{k}\mu_{k}(W_{i}),$

where the weighting function $\lambda$ and contrast vector $c$ are both known.
Suppose the weighting function satisfies $E[\lambda(W_{i})]\neq 0$, and that
the second moments of $\lambda(W_{i})$ and $\mu(W_{i})$ are bounded. Then,
conditional on the controls $W_{1},\dots,W_{N}$, the semiparametric efficiency
bound is almost surely given by

$\mathcal{V}_{\lambda,c}=\frac{1}{E[\lambda(W_{i})]^{2}}E\left[\sum_{k=0}^{K}\frac{\lambda(W_{i})^{2}c_{k}^{2}\sigma_{k}^{2}(W_{i})}{p_{k}(W_{i})}\right].$
(18)

As formalized in the proof (see Section A.2), the efficiency bound
$\mathcal{V}_{\lambda,c}$ establishes a lower bound on the asymptotic variance
of any regular estimator of $\theta_{\lambda,c}$ under the idealized situation
of known propensity scores.171717The 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 optimal
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 to characterize the optimal weights in eq. 19 while remaining
completely agnostic about heterogeneity in treatment effects. See [20] for
additional discussion in the context of a binary treatment.

To establish the second step, we choose $\lambda$ to minimize eq. 18 subject
to the constraint that these weights are non-negative and average to one:
$E[\lambda(W_{i})]=1$. Simple algebra shows that this variance-minimizing
weighting scheme is given by

$\lambda^{*}_{c}(W_{i})=\left(\sum_{k=0}^{K}\frac{c_{k}^{2}\sigma^{2}_{k}(W_{i})}{p_{k}(W_{i})}\right)^{-1}.$
(19)

The asymptotic variance of this easiest-to-estimate weighting,

$\mathcal{V}_{\lambda^{*}_{c},c}=E\left[\left(\sum_{k=0}^{K}\frac{c_{k}^{2}\sigma^{2}_{k}(W_{i})}{p_{k}(W_{i})}\right)^{-1}\right]^{-1},$

is the harmonic mean of
$\sum_{k=0}^{K}\frac{c_{k}^{2}\sigma^{2}_{k}(W_{i})}{p_{k}(W_{i})}$. In
contrast, the efficiency bound for the unweighted contrast is given by the
arithmetic mean
$E\left[\left(\sum_{k=0}^{K}\frac{c_{k}^{2}\sigma^{2}_{k}(W_{i})}{p_{k}(W_{i})}\right)\right]$,
which can be considerably bigger when the propensity scores are not bounded
away from zero or one.

When the contrast vector $c$ is selected to estimate the weighted average
effect of a particular treatment $k$, Proposition 2 implies that the
regression weights are efficient:

######  Corollary 1.

For some $k\geq 1$, let $c^{k}$ be a vector with elements $c^{k}_{j}=1$ if
$j=k$, $c^{k}_{j}=-1$ if $j=0$, and $c^{k}_{j}=0$ otherwise. Suppose that the
conditional variance of relevant potential outcomes is homoskedastic:
$\sigma^{2}_{k}(W_{i})=\sigma^{2}_{0}(W_{i})=\sigma^{2}$. Then the variance-
minimizing weighting scheme is given by $\lambda^{*}_{c^{k}}=\lambda^{k}$,
where

$\lambda^{k}(W_{i})=\frac{p_{0}(W_{i})p_{k}(W_{i})}{p_{0}(W_{i})+p_{k}(W_{i})},$
(20)

with the semiparametric efficiency bound given by

$\mathcal{V}_{\lambda^{k},c^{k}}=\sigma^{2}E\left[\frac{p_{0}(W_{i})p_{k}(W_{i})}{p_{0}(W_{i})+p_{k}(W_{i})}\right]^{-1},$
(21)

where $p_{0}(W_{i})=\Pr(D_{i}=0\mid W_{i})=1-\sum_{k=1}^{K}p_{k}(W_{i})$.

Per eq. 15, the optimal weighting $\lambda^{k}$ coincides with the implicit
weighting of conditional ATEs from the partially linear model (9) when it is
fit only on observations with $D_{i}\in\\{0,k\\}$. 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})}$,
making $\lambda^{k}(W_{i})$ in eq. 20 the conditional variance of the
treatment indicator. Moreover, it follows by standard arguments that
regressing $Y_{i}$ onto $X_{ik}$ and $g(W_{i})$ in the subsample with
$D_{i}\in\\{0,k\\}$ efficiently estimates this weighted average effect
provided $g$ is sufficiently flexible (such as when $g$ is linear and $W_{i}$
consists of group indicators).181818As we discuss in the next subsection, when
the propensity score is unknown, the semiparametric efficiency bound for
estimating $\theta_{\lambda^{k},c^{k}}$ has an additional term relative to eq.
21 arising from the estimation of the optimal weights eq. 20. Thus, while the
regression estimator is semiparametrically efficient, it does not generally
attain the efficiency bound derived in Proposition 2 that assumes a known
propensity score. When the treatment $D_{i}$ is binary, this simply amounts to
running a regression on the binary treatment indicator, with an additive
covariate adjustment.

Corollary 1 thus gives justification for estimating the effect of any given
treatment $k$ by a partially linear regression with an additive covariate
adjustment in the subsample with $D_{i}\in\\{0,k\\}$. To estimate the effects
of all treatments, one runs $K$ such regressions, restricting the sample to
one treatment arm and the control group. Such one-treatment-at-a-time
regressions are simple to implement and do not require explicitly estimating
the propensity score. The regression coefficients are causally interpretable
as weighted averages of conditional treatment effects $\tau_{k}(W_{i})$, so
long as $p_{k}/(p_{0}+p_{k})\in\mathcal{G}$. Moreover, the weighted averages
are locally efficient in the sense of Corollary 1.191919As usual,
homoskedasticity is a tractable baseline: the arguments in favor of ordinary
least squares regression following Corollary 1 can be extended to favor a
(feasible) generalized least squares regression when $\sigma^{2}(W_{i})$ is
known or consistently estimable.

While the robustness property of the one-treatment-at-a-time regression is
well-established, by [4] and subsequent extensions, our efficiency
characterization appears novel. It builds on earlier results in [20, Corollary
5.2] (a working paper version of [21]), and [48, Corollary 1], who show that
the weighting $p_{1}(W_{i})(1-p_{1}(W_{i}))$ is optimal for estimating the
effect of a binary treatment in that the weighting minimizes the asymptotic
variance of a particular class of inverse propensity score weighted
estimators. Corollary 1 above extends the optimality of this weighting to all
regular estimators, and to multiple treatments. Importantly, this result
formalizes a common motivation for using regression to estimate the effects of
a single treatment instead of more involved unconditional ATE estimators: when
treatment effect heterogeneity is minimal or only weakly correlated with the
$\lambda^{*}_{k}(W_{i})$ weights, the regression’s weighted-average effect
will be close to the ATE while being more precisely estimated.

A shortcoming of the optimal weighting scheme in Corollary 1 is that it is
treatment-specific, so comparisons of the “one-at-a-time” weighted-average
effects across treatments are generally not causally
interpretable.202020Formally, 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 issue is
especially salient when the control group is arbitrarily chosen, such as in
teacher VAM regressions which omit an arbitrary teacher from estimation and
seek to make causal comparisons across all teachers.212121Note that this
critique also applies to the own-treatment weights in Proposition 1. Thus even
without contamination bias one may find the implicit multiple-treatment
regression weighting unsatisfying.

We thus turn to the question of how Proposition 2 can be used to select an
efficient weighting scheme that allows for simultaneous comparisons across all
treatment arms. We are interested in reporting estimates of a vector
$\beta_{\lambda^{\text{C}}}$ of $K$ coefficients with elements
$\beta_{\lambda^{\text{C}},k}=\sum_{i=1}^{N}\lambda^{\text{C}}(W_{i})\tau_{k}(W_{i})/\sum_{i=1}^{N}\lambda^{\text{C}}(W_{i})$,
where the weights $\lambda^{\text{C}}$ are common across treatment arms. If we
are equally interested in all $K(K+1)$ contrasts (that is, weighted averages
$\mu_{j}(W_{i})-\mu_{k}(W_{i})$, for all $j\neq k$, $j,k=0,\dotsc,K$), a
natural approach is to choose the weighting scheme $\lambda^{C}$ that
minimizes the average variance across all contrasts:

$\int\mathcal{V}_{\lambda,c}dF(c)=\frac{1}{E[\lambda(W_{i})]^{2}}\sum_{k=0}^{K}\frac{2}{K+1}E\left[\frac{\lambda(W_{i})^{2}\sigma_{k}^{2}(W_{i})}{p_{k}(W_{i})}\right],$

where $F$ gives the uniform distribution over the possible (now random)
contrasts $c$, so that $c_{j}=1$ with probability $1/(K+1)$ and $-1$ with
probability $1/(K+1)$. Minimizing this expression over $\lambda$ is equivalent
to minimizing eq. 18 with $c_{k}^{2}=2/(K+1)$, which leads to the following
result:

######  Corollary 2.

Let $F$ denote the uniform distribution over the possible contrast vectors.
Suppose that $\sigma^{2}_{k}(W_{i})=\sigma^{2}$ for all $k$. Then the
weighting scheme minimizing the average variance bound
$\int\mathcal{V}_{\lambda,c}dF(c)$ is given by

$\lambda^{\text{C}}(W_{i})=1\Big{/}\sum_{k=0}^{K}p_{k}(W_{i})^{-1}.$ (22)

The weights $\lambda^{\text{C}}$ generalize the intuition behind the single
binary treatment (Corollary 1), placing higher weight on covariate strata
where the treatments are evenly distributed, and putting less weight on strata
with limited overlap. When the treatment is binary, $K=1$, the weights reduce
to the one-at-a-time weights in Corollary 1,
$\lambda^{C}(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^{C}$ remain the same
for every treatment, allowing for simultaneous comparisons across all
treatment pairs $(k,\ell)$. We next consider estimating $\beta_{\lambda^{C}}$
using a weighted regression approach.

###  4.3 Estimating Efficiently Weighted Average Effects

If the propensity scores $p(W_{i})$ were known, one could estimate
$\beta_{\lambda^{\text{C}}}$ by a weighted regression of $Y_{i}$ onto $X_{i}$
and a constant, with each observation weighted by
${\lambda}^{\text{C}}(W_{i})/p_{D_{i}}(W_{i})$. When the treatment is binary,
this estimator reduces to the estimator studied in [48]. Since the propensity
score is unknown, we replace the infeasible weights with feasible weights
$\hat{\lambda}^{\text{C}}(W_{i})/\hat{p}_{D_{i}}(W_{i})$, where
$\hat{p}_{k}(W_{i})$ is a feasible estimate of the propensity score and
$\hat{\lambda}^{\text{C}}(W_{i})=1/\sum_{k=0}^{K}1/\hat{p}_{k}(W_{i})$. When
$\mathcal{G}$ is finite-dimensional, we may use the regression estimator that
projects $X_{ik}$ onto $g(W_{i})$:

$\hat{p}_{k}(W_{i})=\arg\min_{\tilde{p}\in\mathcal{G}}\sum_{i=1}^{N}\left(X_{ik}-\tilde{p}(W_{i})\right)^{2}.$

With linear $\mathcal{G}$, for example, $\hat{p}_{k}(W_{i})$ is simply the
fitted value from a linear regression of $X_{ik}$ on $W_{i}$ and a constant.
The resulting estimator can be written as

$\hat{\beta}_{\hat{\lambda}^{\text{C}},k}=\frac{1}{\sum_{i=1}^{N}\frac{\hat{\lambda}^{\text{C}}(W_{i})}{\hat{p}_{k}(W_{i})}X_{ik}}\sum_{i=1}^{N}\frac{\hat{\lambda}^{\text{C}}(W_{i})}{\hat{p}_{k}(W_{i})}X_{ik}Y_{i}-\frac{1}{\sum_{i=1}^{N}\frac{\hat{\lambda}^{\text{C}}(W_{i})}{\hat{p}_{0}(W_{i})}X_{i0}}\sum_{i=1}^{N}\frac{\hat{\lambda}^{\text{C}}(W_{i})}{\hat{p}_{0}(W_{i})}X_{i0}Y_{i}.$
(23)

When the treatment is binary and $\mathcal{G}$ is linear, this weighted
regression estimator coincides with the usual (unweighted) regression
estimator that regresses $Y_{i}$ onto $D_{i}$ and $W_{i}$.222222To 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}^{\text{C}},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}$.

We now show that the estimator $\hat{\beta}_{\hat{\lambda}^{\text{C}}}$ is
efficient in the sense that is achieves the semiparametric efficiency bound
for estimating $\beta_{\lambda^{\text{C}}}$:

######  Proposition 3.

Suppose eq. 12 holds in an i.i.d. sample of size $N$, with known non-
degenerate propensity scores $p_{k}(W_{i})$. Let
$\beta^{*}_{\lambda^{C},k}=E[\lambda^{\text{C}}(W_{i})\tau_{k}(W_{i})]/E[\lambda^{\text{C}}(W_{i})]$,
and
$\alpha^{*}_{k}=\beta^{*}_{\lambda^{C},k}+E[\lambda^{\text{C}}(W_{i})\mu_{0}(W_{i})]/E[\lambda^{\text{C}}(W_{i})]$.
Suppose that the fourth moments of $\lambda^{C}(W_{i})$ and $\mu(W_{i})$ are
bounded, and that $p_{k}\in\mathcal{G}$,
$(\mu_{k}(W_{i})-\alpha^{*}_{k})\frac{\lambda^{\text{C}}(W_{i})^{2}}{p_{k^{\prime}}(W_{i})^{2}}\in\mathcal{G}$,
and
$(\mu_{k}(W_{i})-\alpha^{*}_{k})\frac{\lambda^{\text{C}}(W_{i})}{p_{k}(W_{i})}\in\mathcal{G}$
for all $k,k^{\prime}$. Then, provided it is asymptotically linear and
regular, $\hat{\beta}_{\hat{\lambda}^{\text{C}}}$ achieves the semiparametric
efficiency bound for estimating ${\beta}_{\lambda^{\text{C}}}$, with diagonal
elements of its asymptotic variance of:

$\frac{1}{E[\lambda^{\text{C}}(W_{i})]^{2}}E\left[\frac{\lambda^{\text{C}}(W_{i})^{2}\sigma_{0}^{2}(W_{i})}{p_{0}(W_{i})}+\frac{\lambda^{\text{C}}(W_{i})^{2}\sigma_{k}^{2}(W_{i})}{p_{k}(W_{i})}\right.\\\
\left.+{\lambda^{\text{C}}}(W_{i})^{2}(\tau_{k}(W_{i})-\beta^{*}_{\lambda^{\text{C}},k})^{2}\left(\sum_{k^{\prime}=0}^{K}\frac{{\lambda^{\text{C}}}(W_{i})^{2}}{p_{k}(W_{i})^{3}}-1\right)\right].$

This efficiency result doesn’t rely on homoskedasticity: under
heteroskedasticity, the estimator $\hat{\beta}_{\hat{\lambda}^{\text{C}}}$ is
still efficient for ${\beta}_{\lambda^{\text{C}}}$ (although the weighting
$\lambda^{\text{C}}(W_{i})$ is itself not optimal under heteroskedasticity).
Proposition 3 is stated under the high-level condition that
$\hat{\beta}_{\hat{\lambda}^{\text{C}}}$ is regular; the proof uses the
calculations in [55] to verify that the estimator achieves the efficiency
bound. Primitive regularity conditions will depend on the form of
$\mathcal{G}$ and are omitted for brevity.

######  Remark 4.

The asymptotic variance of the estimator $\hat{\beta}_{\lambda^{\text{C}}}$ is
larger than the asymptotic variance of the infeasible estimator that replaces
the estimated weights $\hat{\lambda}^{\text{C}}(W_{i})/\hat{p}_{D_{i}}(W_{i})$
in eq. 23 with the infeasible weights
${\lambda}^{\text{C}}(W_{i})/p_{D_{i}}(W_{i})$. The latter achieves the
asymptotic variance implied by Corollary 2,

$\frac{1}{E[\lambda^{\text{C}}(W_{i})]^{2}}E\left[\frac{\lambda^{\text{C}}(W_{i})^{2}\sigma_{0}^{2}(W_{i})}{p_{0}(W_{i})}+\frac{\lambda^{\text{C}}(W_{i})^{2}\sigma_{k}^{2}(W_{i})}{p_{k}(W_{i})}\right].$
(24)

The extra term of the asymptotic variance in Proposition 3 relative to eq. 24
reflects the cost of having to estimate the weights.232323The extra term shows
this cost is zero if either there is no treatment effect heterogeneity, so
that $\tau_{k}(W_{i})=\beta^{*}_{\lambda^{C},k}$, or if the treatment
assignment is completely randomized so that $p_{k}(W_{i})=1/(K+1)$. In the
latter case $\lambda^{*}(W_{i})=1/(K+1)^{2}$ so
$\sum_{k=0}^{K}\lambda^{\text{C}}(W_{i})^{2}/p(W_{i})^{3}=1$. The extra term
can be avoided altogether if we interpret
$\hat{\beta}_{\hat{\lambda}^{\text{C}}}$ as an estimator of
${\beta}_{\hat{\lambda}^{\text{C}}}$. This follows from arguments in [20,
Lemma B.6].

##  5 Practical Guidance and Application

###  5.1 Measuring and Avoiding Contamination Bias

A researcher interested in estimating the effects of multiple dependent
treatments with regression can use Proposition 1 to measure the extent of
contamination bias in her estimates. When the treatment assignment is
conditionally ignorable, she can further compute one of the three alternative
estimators discussed in Section 4.242424We again refer readers to [23, 62, 16,
15, 67] for solutions under a parallel trends assumption. Here we provide
practical guidance on both procedures, which we illustrate in an application
in the next subsection.

For simplicity, we focus on the case where $g$ is linear and eq. 9 is
estimated by ordinary least squares (OLS). We assume Assumption 1 and both
conditions in Assumption 2 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.

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

$Y_{i}=\alpha+\sum_{k=1}^{K}X_{ik}\beta_{k}+W_{i}^{\prime}\gamma+U_{i},$ (25)

and obtain a sample analog of the decomposition in Proposition 1. To this end,
note that the own-treatment and contamination bias weights in Proposition 1
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}^{\prime}]^{-1}E[\tilde{X}_{i}X_{i}^{\prime}\mid
W_{i}].$

An estimate of this weight matrix is given by the sample analog:

$\hat{\Lambda}_{i}=(\dot{X}^{\prime}\dot{X})^{-1}\dot{X}_{i}X_{i}^{\prime},$

where $\dot{X}_{i}$ is the sample residual from an 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 Proposition 1; for $k\neq\ell$ this is an estimate of a
contamination weight.

Under linearity, the $k$th conditional ATEs 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

$Y_{i}=\alpha_{0}+\sum_{k=1}^{K}X_{ik}\gamma_{0,k}+W_{i}^{\prime}\alpha_{W,0}+\sum_{k=1}^{K}X_{ik}W_{i}^{\prime}\gamma_{W,k}+\dot{U}_{i}.$
(26)

Estimating eq. 26 by OLS yields estimates
$\hat{\tau}_{k}(W_{i})=\hat{\gamma}_{0,k}+W_{i}^{\prime}\hat{\gamma}_{W,k}$.
For each observation $i$, we stack the set of conditional ATE estimates in a
$K\times 1$ vector $\hat{\tau}(W_{i})$.

Using OLS normal equations, we then obtain the exact decomposition

$\hat{\beta}=\sum_{i=1}^{N}\operatorname{diag}(\hat{\Lambda}_{i})\hat{\tau}(W_{i})+\sum_{i=1}^{N}[\hat{\Lambda}_{i}-\operatorname{diag}(\hat{\Lambda}_{i})]\hat{\tau}(W_{i}),$
(27)

which is the sample analog of the population decomposition in Proposition 1.
The first term in this decomposition 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 eq. 27 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})$.252525When the covariates are not
saturated, it is possible that the estimated weighting function
$\hat{\Lambda}(w)=\frac{1}{N}\sum_{i=1}^{N}\operatorname{\mathbbm{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}\operatorname{\mathbbm{1}}\\{W_{i}=w\\}\hat{\Lambda}_{i,k\ell}}{\sum_{i=1}^{N}\operatorname{\mathbbm{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 or a low correlation between such heterogeneity and the
contamination weights.

Implementing the alternative estimators from Section 4 is also straightforward
under the linearity assumptions. For the first solution, estimating the
interacted regression

$Y_{i}=\alpha_{0}+\sum_{k=1}^{K}X_{ik}\tau_{k}+W_{i}^{\prime}\alpha_{W,0}+\sum_{k=1}^{K}X_{ik}(W_{i}-\mkern
1.5mu\overline{\mkern-1.5muW\mkern-1.5mu}\mkern
1.5mu)^{\prime}\gamma_{W,k}+\dot{U}_{i}.$ (28)

by OLS yields estimates of the unweighted ATEs $\tau_{k}=E[\tau_{k}(W_{i})]$.
Here $\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.
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 OLS estimates of eq. 26.

The second solution is to estimate the uninteracted regression,

$Y_{i}=\ddot{\alpha}_{k}+X_{ik}\ddot{\beta}_{k}+W_{i}^{\prime}\ddot{\gamma}_{k}+\ddot{U}_{ik}$
(29)

among observations assigned either to treatment $k$ or the control group,
$D_{i}\in\\{0,k\\}$, for each of the treatments $k=1,\dotsc,K$. These one-
treatment-at-a-time regressions estimate convex weighted averages of treatment
effects, with weights that are efficient under homoskedasticity (in the sense
of corollary 1) but which will generally vary across the different treatments.
This can make comparisons across treatment arms difficult.

The third solution is to estimate an efficiently weighted average of the
conditional ATEs, with weights that are constrained to be common across
treatments. Under linearity, we can estimate the common weights
$\lambda^{\text{C}}$ as

$\hat{\lambda}^{\text{C}}(W_{i})=\left(\sum_{k=0}^{K}\hat{p}_{k}(W_{i})^{-1}\right)^{-1},$
(30)

where $\hat{p}_{k}(W_{i})=X_{ik}-\dot{X}_{ik}$ denote estimated propensity
scores. We then regress $Y_{i}$ on $X_{i}$, weighting each observation by
$\hat{\lambda}^{C}(W_{i})/\hat{p}_{D_{i}}(W_{i})$.

While the second and third solutions may yield more precise estimates than the
equal-weighted ATE estimates, the gains in precision are achieved by changing
the estimand to a different convex average of conditional treatment effects.
In particular, covariate values $w$ where the propensity score $p_{k}(w)$ is
close to zero for some $k$ will be effectively discarded.

On the other hand, if the conditional treatment effects $\tau(W_{i})$ are
approximately independent of the propensity scores $p(W_{i})$, the weighting
scheme may have little effect on the estimands, even if the treatment effect
heterogeneity is substantial. In such cases, we also expect the contamination
bias to be small, since the contamination weights are a function of the
propensity scores. We next investigate this possibility in our application.

###  5.2 Application

We illustrate the potential for contamination bias using data from Project
STAR, as analyzed in [45]. The Project STAR 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.262626Students 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 [45]. All analyses in this section are conducted with our Stata
package, multe, which researchers can use to gauge the extent of contamination
bias in similar applications.

Table 1: Project STAR Contamination Bias and Treatment Effect Estimates

| A. Contamination Bias Estimates
---|---
| Regression | Own | Bias | Worst-Case Bias
| Coefficient | Effect | Negative |  Positive
| (1) | (2) | (3) | (4) |  (5)
Small Class Size | 5.357 | 5.202 | 0.155 | -1.654 |  1.670
| (0.778) | (0.778) | (0.160) | (0.185) |  (0.187)
Teaching Aide | 0.177 | 0.360 | -0.183 | -1.529 |  1.530
| (0.720) | (0.714) | (0.149) | (0.176) |  (0.177)
| B. Treatment Effect Estimates
|  | Unweighted | Efficiently-Weighted |
|  | (ATE) | One-at-a-time | Common |
|  | (1) | (2) | (3) |
Small Class Size |  | 5.561 | 5.295 | 5.563 |
|  | (0.763) | (0.775) | (0.764) |
|  | [0.744] | [0.743] | [0.742] |
Teaching Aide |  | 0.070 | 0.263 | -0.003 |
|  | (0.708) | (0.715) | (0.712) |
|  | [0.694] | [0.691] | [0.695] |

_Notes:_ Panel A estimates the contamination bias and range of potential
contamination bias in regression estimates of small class and teaching aide
treatment effects for the Project STAR kindergarten analysis. The analysis
sample includes 5,868 students. Column 1 reports estimates from a partially
linear model in eq. 25. Columns 2 and 3 estimate the own- and cross-treatment
decomposition of this estimate in eq. 27. Columns 4 and 5 reports the smallest
(largest) possible contamination bias from reordering the conditional ATEs to
be as negatively (positively) correlated with the cross-treatment weights as
possible. Panel B summarizes estimates of small class and teaching aide
treatment effects from different specifications in the kindergarten sample of
Project STAR. Column 1 reports estimates of small class and teaching aide
treatment effects from the interacted model in eq. 28. Column 2 reports
estimates from the treatment-specific regressions in eq. 29. Column 3 reports
estimates from the efficiently weighted specification, using the estimated
weights in eq. 30. Robust standard errors are reported in parentheses.
Standard errors that assume the propensity scores are known are reported in
square brackets.

Column 1 of Panel A in Table 1 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 uninteracted
regression in eq. 25, 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 [45] 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 [45, column 6 of Table V, panel A], 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).272727Our sample and
estimates are very similar to—but not exactly the same as—those in [45]. We
use robust (non-clustered) standard errors throughout this analysis, since the
randomization of students to classrooms is at the individual level [1].
Results are similar when we cluster by classroom.

As discussed in Section 2, treatment assignment probabilities vary across the
schools, indicated by the fixed effects in $W_{i}$. If treatment effects also
vary across schools, and if this variation is correlated 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). Net of any contamination bias, we expect each $\beta_{k}$ to
identify a weighted average of own treatment effects $\tau_{k}(W_{i})$, with
convex weights given by $\lambda_{kk}(W_{i})$.

Columns 2 and 3 of Table 1 apply the decomposition in eq. 27 to the regression
coefficients in column 1. The contamination bias appears to be minimal. The
small class size regression estimate of 5.36 is composed of a weighted average
of small class size treatment effects equalling 5.20 and a weighted average of
teaching aide treatment effects equalling 0.16. Similarly, the teaching aide
regression coefficient of 0.18 decomposes into a weighted average of teaching
aide treatment effects equalling 0.36 and a weighted average of small class
size effects equalling $-0.18$. Netting out the contamination bias estimate
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 lack of meaningful bias in the regression estimates of Project STAR
effects is due to a weak correlation between the conditional treatment effects
$\tau(W_{i})$ and the contamination weights. These correlations are shown in
Figure 2, which plots estimates of the school-specific treatment effects
$\tau_{k}(W_{i})$ against the own-treatment and contamination weights
$\lambda_{kk}(W_{i})$ and $\lambda_{k\ell}(W_{i})$ for $\ell\neq k$. Panels A
and C show that the own-treatment weight correlation is negative for the small
class size treatment ($-0.19$) and positive for the aide treatment ($0.25$).
The partially linear regression model’s estimate of own-treatment effects
(column 2 of Table 1) thus understates the average small class size effect and
overstates the average aide effect, relative to the ATE. Panels B and D
further show that correlation between estimated cross-treatment effects and
weights is positive for the small class effect estimate (0.10) and negative
for the aide effect estimate ($-0.13$). There is thus positive contamination
bias in the partially linear regression model’s estimate of small class size
effects and negative contamination bias in the regression’s aide effect
estimate, as shown in column 2 of Table 1. But neither set of correlations is
strong enough to meaningfully bias the estimates.

Panel A: Small Class
Own-Treatment Weight

Panel B: Aide
Cross-Treatment Weight

Panel C: Aide
Own-Treatment Weight

Panel D: Small Class
Cross-Treatment Weight

Figure 2: Project STAR Treatment Effects and Regression Weights

_Note:_ This figure shows correlations between estimated school-specific
treatment effects and the implicit school-specific regression weights in
column 1 of Table 1. Panels A and B show correlations for the decomposition of
the small class treatment effect estimate in columns 2 and 3 of Table 1. Panel
A plots the estimated small class treatment effects by school against the
estimated own-treatment weights, while Panel B plots the estimated teaching
aide treatment effects by school against the estimated cross-treatment
weights. Panels C and D show analogous correlations for the decomposition of
the teaching aide treatment effect estimate in columns 2 and 3 of Table 1.
Panel C plots the estimated teaching aide treatment effects by school against
the estimated own-treatment weights, while Panel D plots the estimated small
class treatment effects by school against the estimated cross-treatment
weights. Correlations and lines of best fit are reported on each panel. The
size of the points is proportional to the number of students enrolled in each
school.

Importantly, Figure 2 shows that the lack of contamination bias is not due to
a lack of treatment effect heterogeneity across schools. There is considerable
variation along the $y$-axis of each plot. Adjusting for estimation error, we
find a standard deviation of $\tau_{k}(W_{i})$ across the schools indexed by
$W_{i}$ of 12.7 for the small class treatment and of 10.9 for the aide
treatment.282828We 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 Table 1. Thus, had the
experimental design been such that the contamination weights strongly
correlate with this variation, sizeable contamination bias could have
resulted. In practice, the variation in school-specific propensity scores
$p_{k}(W_{i})$ appears to have been largely unrelated to school-specific
treatment effects.292929The own-treatment weights in Figure 2 are highly
correlated with the respective treatment propensity score. For the small class
size (teaching aide) treatment this correlation is 0.92 (0.73).

To illustrate the potential for contamination bias in this setting, we compute
worst-case (positive and negative) weighted averages of the estimated
$\tau_{k}(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 heterogeneity in treatment effects. Columns 4 and 5 in Table 1 shows that
both bounds on possible contamination bias are large relative to the effect
estimates: $[-1.654,1.670]$ for the small class size treatment and
$[-1.529,1.530]$ for the teaching aide treatment. Thus, while actual
contamination bias for both treatments is an order of magnitude smaller, the
underlying heterogeneity in this setting makes severe contamination bias
possible.

Panel B of Table 1 illustrates the three solutions to the contamination bias
problem discussed in Section 4. Column 1 estimates the unweighted ATEs of the
small class size and teaching aide treatment, by estimating the interacted
regression specification in eq. 28. Column 2 estimates the one-treatment-at-a-
time regressions in eq. 29 for $k=1,2$. Finally, column 3 estimates the
efficiently-weighted ATEs of each treatment, by running a weighted regression
of $Y_{i}$ onto $X_{i}$, using the common weight estimates in eq. 30.

As discussed in Remark 4, the optimal weighting schemes underlying the
estimates in columns 2 and 3 of Panel B are derived under the assumption that
the propensity scores are known. To gauge the relative importance of this
assumption, Panel B also reports a version of the standard errors computed
under the assumption that the sample treatment probabilities in each school
match the true propensity scores.303030This is the case under stratified block
randomization, where a fixed proportion of students is assigned to the two
treatments. In contrast, sample treatment proportions need not match the true
propensity scores under a Bernoulli trial where each student is assigned to
treatments according to a coin toss. This changes the standard errors little,
showing that there is minimal cost to estimating the optimal weights.313131The
standard errors reported in parentheses in Panel B are valid for the
population analogs $\beta_{k}$ and $\beta_{\lambda^{\text{C}}}$, i.e.
$E[\lambda^{k}(W_{i})\tau_{k}(W_{i})]/E[\lambda^{k}(W_{i})]$ and
$E[\lambda^{\text{C}}(W_{i})\tau_{k}(W_{i})]/E[\lambda^{\text{C}}(W_{i})]$.
Since these standard errors are potentially conservative when viewed as
standard errors for $\beta_{k}$ and $\beta_{\lambda^{\text{C}}}$, the standard
error comparison gives an upper bound on the cost to estimating the optimal
weights.

There turns out to be little difference between the partially linear model
estimates of Project STAR treatment effects and these alternative estimates.
In columns 1 and 2 of Panel B we estimate a small class size effect of 5.56
and 5.30, which are close to the 5.36 estimate in column 1 of Panel A.
Teaching aide effect estimates are also similar: 0.07 and 0.26 in columns 1
and 2 of Panel B, compared to 0.18 in column 1 of Panel A. The efficiently
weighted estimates in column 3 of Panel B are again similar: 5.56 for the
small class size treatment and $0.00$ for the teaching aide treatment.
Interestingly, the standard errors are roughly constant across the columns,
regardless of whether the propensity score is treated as known.

##  6 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 heterogeneous effects,
with coefficients on each treatment generally contaminated by the effects of
other treatments. We provide intuition for why the influential result of [4]
fails to generalize to multiple treatments, and show how the contamination
bias problem connects to a recent literature studying DiD regressions and
related estimators. We discuss three alternative estimators that are free of
this bias, including a new estimator that efficiently weights conditional
average treatment effects. The analysis underling this estimator also
formalizes a virtue of regression adjustment in the binary treatment case: the
weighting that it implicitly uses to combine heterogeneous treatment effects
minimizes the semiparametric efficiency bound for convex weighted averages of
ATEs.

Our application to Project STAR shows that significant contamination bias
could arise in RCTs when there is significant treatment effect heterogeneity.
Whether the bias _does_ arise, however, depends on the correlation between
effect heterogeneity and the contamination weights we derive in our
theoretical analysis. Researchers can estimate this correlation, and report it
alongside the alternative estimates that are free of contamination bias. Such
investigation reveals whether, as in our application, the results based on
alternative estimators are more similar than the worst-case bounds implied by
the theory. Broadly, our analysis highlights the importance of testing the
empirical relevance of theoretical concerns with how regression combines
heterogeneous effects.

We expect the tools in this paper to be especially relevant in modern RCT
designs that generate substantial variation in treatment propensity scores to
maximize efficiency [[, e.g.]]tabord-meehan. Propensity scores are also likely
to vary dramatically in quasi-experimental analyses, such as with teacher
VAMs, where a large number of covariates are needed to make the conditionally
ignorability of treatment plausible. Contamination bias diagnostics can be a
useful tool for ensuring the reliability and robustness of regression
estimates in such settings.

## References

  * [1] Alberto Abadie et al.  “When Should You Adjust Standard Errors for Clustering?”, 2017  DOI: 10.3386/w24003
  * [2] Jason Abaluck et al.  “Mortality effects and choice across private health insurance plans”  In _The Quarterly Journal of Economics_ 136.3, 2021, pp. 1557–1610  DOI: 10.1093/qje/qjab017
  * [3] Joshua Angrist et al.  “Credible School Value-Added with Undersubscribed School Lotteries”  In _The Review of Economics and Statistics_ , 2021  DOI: 10.1162/rest_a_01149
  * [4] Joshua D. Angrist  “Estimating the Labor Market Impact of Voluntary Military Service Using Social Security Data on Military Applicants”  In _Econometrica_ 66.2, 1998, pp. 249–288  DOI: 10.2307/2998558
  * [5] Joshua D. Angrist, Guido W. Imbens and Alan B. Krueger  “Jackknife Instrumental Variables Estimation”  In _Journal of Applied Econometrics_ 14.1, 1999, pp. 57–67  DOI: 10.1002/(SICI)1099-1255(199901/02)14:1<57::AID-JAE501>3.0.CO;2-G
  * [6] Joshua D. Angrist and Alan B. Krueger  “Empirical Strategies in Labor Economics”  In _Handbook of Labor Economics_ 3A Amsterdam: Elsevier, 1999, pp. 1277–1366  DOI: 10.1016/S1573-4463(99)03004-7
  * [7] Joshua D. Angrist and Jorn-Steffen Pischke  “Mostly Harmless Econometrics: An Empiricist’s Companion”  Princeton, NJ: Princeton University Press, 2009  DOI: 10.2307/j.ctvcm4j72
  * [8] Joshua D. Angrist et al.  “Leveraging lotteries for school value-added: Testing and estimation”  In _The Quarterly Journal of Economics_ 132.2, 2017, pp. 871–919  DOI: 10.1093/qje/qjx001
  * [9] Peter M. Aronow and Cyrus Samii  “Does Regression Produce Representative Estimates of Causal Effects?”  In _American Journal of Political Science_ 60.1, 2016, pp. 250–267  DOI: 10.1111/ajps.12185
  * [10] Susan Athey and Guido W. Imbens  “Design-based analysis in difference-in-differences settings with staggered adoption”  In _Journal of Econometrics_ 226.1, 2022, pp. 62–79  DOI: 10.1016/j.jeconom.2020.10.012
  * [11] Luc Behaghel, Bruno Crépon and Marc Gurgand  “Robustness of the encouragement design in a two-treatment randomized control trial”, 2013
  * [12] Abraham Berman and Robert J. Plemmons  “Nonnegative Matrices in the Mathematical Sciences”, Classics in Applied Mathematics 9  Philadelphia, PA: Society for IndustrialApplied Mathematics, 1994  DOI: 10.1137/1.9781611971262
  * [13] Manudeep Bhuller and Henrik Sigstad  “2SLS with Multiple Treatments”, 2022  arXiv:2205.07836
  * [14] Christine Blandhol et al.  “When is TSLS Actually LATE?”, 2022  DOI: 10.2139/ssrn.4014707
  * [15] Kirill Borusyak, Xavier Jaravel and Jann Spiess  “Revisiting Event Study Designs: Robust and Efficient Estimation”, 2022  arXiv:2108.12419
  * [16] Brantly Callaway and Pedro H.C. Sant’Anna  “Difference-in-differences with multiple time periods”  In _Journal of Econometrics_ 225.2, 2021, pp. 200–230  DOI: 10.1016/j.jeconom.2020.12.001
  * [17] Matias D. Cattaneo  “Efficient Semiparametric Estimation of Multi-Valued Treatment Effects under Ignorability”  In _Journal of Econometrics_ 155.2, 2010, pp. 138–154  DOI: 10.1016/j.jeconom.2009.09.023
  * [18] Victor Chernozhukov, Whitney K. Newey and Rahul Singh  “Automatic Debiased Machine Learning of Causal and Structural Effects”, 2021  arXiv:1809.05224
  * [19] Raj Chetty, John N. Friedman and Jonah E. Rockoff  “Measuring the impacts of teachers I: Evaluating bias in teacher value-added estimates”  In _American Economic Review_ 104.9, 2014, pp. 2593–2632  DOI: 10.1257/aer.104.9.2593
  * [20] Richard K. Crump et al.  “Moving the Goalposts: Addressing Limited Overlap in the Estimation of Average Treatment Effects by Changing the Estimand”, 2006  DOI: 10.3386/t0330
  * [21] Richard K. Crump et al.  “Dealing with limited overlap in estimation of average treatment effects”  In _Biometrika_ 96.1, 2009, pp. 187–199  DOI: 10.1093/biomet/asn055
  * [22] Clément de Chaisemartin and Xavier D’Haultfœuille  “Two-Way Fixed Effects Estimators with Heterogeneous Treatment Effects”  In _American Economic Review_ 110.9, 2020, pp. 2964–2996  DOI: 10.1257/aer.20181169
  * [23] Clément de Chaisemartin and Xavier D’Haultfœuille  “Two-way Fixed Effects Regressions with Several Treatments”, 2022  arXiv:2012.10077
  * [24] Will Dobbie and Jae Song  “Debt relief and debtor outcomes: Measuring the effects of consumer bankruptcy protection”  In _American Economic Review_ 105.3, 2015, pp. 1272–1311  DOI: 10.1257/aer.20130612
  * [25] Brigham Frandsen, Lars Lefgren and Emily Leslie  “Judging judge fixed effects”, 2019  DOI: 10.3386/w25528
  * [26] David A. Freedman  “On Regression Adjustments in Experiments with Several Treatments”  In _The Annals of Applied Statistics_ 2.1, 2008, pp. 176–196  DOI: 10.1214/07-AOAS143
  * [27] David A. Freedman  “On Regression Adjustments to Experimental Data”  In _Advances in Applied Mathematics_ 40.2, 2008, pp. 180–193  DOI: 10.1016/j.aam.2006.12.003
  * [28] Roland G. Fryer and Steven D. Levitt  “Testing for Racial Differences in the Mental Ability of Young Children”  In _American Economic Review_ 103.2, 2013, pp. 981–1005  DOI: 10.1257/aer.103.2.981
  * [29] Michael Geruso, Timothy Layton and Jacob Wallace  “Are all managed care plans created equal? Evidence from random plan assignment in Medicaid”, 2020  DOI: 10.3386/w27762
  * [30] Andrew Goodman-Bacon  “Difference-in-Differences with Variation in Treatment Timing”  In _Journal of Econometrics_ 225.2, 2021, pp. 254–277  DOI: 10.1016/j.jeconom.2021.03.014
  * [31] Bryan S. Graham and Cristine Campos de Xavier Pinto  “Semiparametrically Efficient Estimation of the Average Linear Regression Function”  In _Journal of Econometrics_ 226.1, 2022, pp. 115–138  DOI: 10.1016/j.jeconom.2021.07.008
  * [32] Jinyong Hahn  “On the Role of the Propensity Score in Efficient Semiparametric Estimation of Average Treatment Effects”  In _Econometrica_ 66.2, 1998, pp. 315–331  DOI: 10.2307/2998560
  * [33] Keisuke Hirano, Guido W. Imbens and Geert Ridder  “Efficient Estimation of Average Treatment Effects Using the Estimated Propensity Score”  In _Econometrica_ 71.4, 2003, pp. 1161–1189  DOI: 10.1111/1468-0262.00442
  * [34] Peter D. Hull  “Estimating hospital quality with quasi-experimental data”, 2018  DOI: 10.2139/ssrn.3118358
  * [35] Peter D. Hull  “Estimating Treatment Effects in Mover Designs”, 2018  arXiv:1804.06721
  * [36] Peter D. Hull  “IsoLATEing: Identifying Counterfactual-Specific Treatment Effects with Cross-Stratum Comparisons”, 2018  DOI: 10.2139/ssrn.2705108
  * [37] Guido W. Imbens and Joshua D. Angrist  “Identification and Estimation of Local Average Treatment Effects”  In _Econometrica_ 62.2, 1994, pp. 467–475  DOI: 10.2307/2951620
  * [38] Guido W. Imbens and Jeffrey Wooldridge  “Recent developments in the econometrics of program evaluation”  In _Journal of Economic Literature_ 47.1, 2009, pp. 5–86  DOI: 10.1257/jel.47.1.5
  * [39] Thomas J. Kane and Douglas O. Staiger  “Estimating Teacher Impacts on Student Achievement: An Experimental Evaluation”, 2008  DOI: 10.3386/w14607
  * [40] Shakeeb Khan and Elie Tamer  “Irregular Identification, Support Conditions, and Inverse Weight Estimation”  In _Econometrica_ 78.6, 2010, pp. 2021–2042  DOI: 10.3982/ECTA7372
  * [41] Lars J. Kirkeboen, Edwin Leuven and Magne Mogstad  “Field of study, earnings, and self-selection”  In _The Quarterly Journal of Economics_ 131.3, 2016, pp. 1057–1111  DOI: 10.1093/qje/qjw019
  * [42] Patrick Kline and Christopher R. Walters  “Evaluating Public Programs with Close Substitutes: The Case of Head Start”  In _Quarterly Journal of Economics_ 131.4, 2016, pp. 1795–1848  DOI: 10.1093/qje/qjw027
  * [43] Jeffrey R. Kling  “Incarceration length, employment, and earnings”  In _American Economic Review_ 96.3, 2006, pp. 863–876  DOI: 10.1257/aer.96.3.863
  * [44] Michal Kolesár  “Estimation in an instrumental variables model with treatment effect heterogeneity”  In _Unpublished Working Paper_ , 2013
  * [45] Alan B. Krueger  “Experimental estimates of education production functions”  In _The Quarterly Journal of Economics_ 114.2, 1999, pp. 497–532  DOI: 10.1162/003355399556052
  * [46] Alan B. Krueger and Lawrence H. Summers  “Efficiency Wages and the Inter-Industry Wage Structure”  In _Econometrica_ 56.2, 1988, pp. 259–293  DOI: 10.2307/1911072
  * [47] Sokbae Lee and Bernard Salanié  “Identifying Effects of Multivalued Treatments”  In _Econometrica_ 86.6, 2018, pp. 1939–1963  DOI: 10.3982/ECTA14269
  * [48] Fan Li, Kari Lock Morgan and Alan M. Zaslavsky  “Balancing covariates via propenstiy score weighting”  In _Journal of the American Statistical Association_ 113.521, 2018, pp. 390–400  DOI: 10.1080/01621459.2016.1260466
  * [49] Nicole Maestas, Kathleen J. Mullen and Alexander Strand  “Does disability insurance receipt discourage work? Using examiner assignment to estimate causal effects of SSDI receipt”  In _American Economic Review_ 103.5, 2013, pp. 1797–1829  DOI: 10.1257/aer.103.5.1797
  * [50] Willard G. Manning et al.  “Health insurance and the demand for medical care: evidence from a randomized experiment”  In _American Economic Review_ 77.3, 1987, pp. 251–277
  * [51] Brad McNeney and Jon A. Wellner  “Application of Convolution Theorems in Semiparametric Models with Non-i.i.d. Data”  In _Journal of Statistical Planning and Inference_ 91.2, 2000, pp. 441–480  DOI: 10.1016/S0378-3758(00)00193-2
  * [52] Magne Mogstad, Alexander Torgovitsky and Christopher R. Walters  “The Causal Interpretation of Two-Stage Least Squares with Multiple Instrumental Variables”  In _American Economic Review_ 111.11, 2021, pp. 3663–3698  DOI: 10.1257/aer.20190221
  * [53] Jack Mountjoy and Brent Hickman  “The Returns to College(s): Estimating Value-Added and Match Effects in Higher Education”, 2020  DOI: 10.2139/ssrn.3537773
  * [54] Michael Mueller-Smith  “The criminal and labor market impacts of incarceration” Unpubished manuscript, University of Michigan, 2015
  * [55] Whitney K. Newey  “The Asymptotic Variance of Semiparametric Estimators”  In _Econometrica_ 62.6, 1994, pp. 1349–1382  DOI: 10.2307/2951752
  * [56] Samuel Norris  “Examiner inconsistency: Evidence from refugee appeals”, 2019  DOI: 10.2139/ssrn.3267611
  * [57] James M. Robins, Steven D. Mark and Whitney K. Newey  “Estimating Exposure Effects by Modelling the Expectation of Exposure Conditional on Confounders”  In _Biometrics_ 48.2, 1992, pp. 479  DOI: 10.2307/2532304
  * [58] Peter Michael Robinson  “Root-$N$-Consistent Semiparametric Regression”  In _Econometrica_ 56.4, 1988, pp. 931  DOI: 10.2307/1912705
  * [59] Jonathan Roth et al.  “What’s Trending in Difference-in-Differencees? A Synthesis of the Recent Econometrics Literature”, 2022  arXiv:2201.01194
  * [60] Jonathan Skinner  “Chapter Two - Causes and Consequences of Regional Variations in Health Care”  In _Handbook of Health Economics_ 2, Handbook of Health Economics  Elsevier, 2011, pp. 45–93  DOI: 10.1016/B978-0-444-53592-4.00002-5
  * [61] Tymon Słoczyński  “Interpreting OLS Estimands When Treatment Effects are Heterogeneous: Smaller Groups Get Larger Weights”  In _Review of Economics and Statistics_ , 2022, pp. 1–9  DOI: 10.1162/rest_a_00953
  * [62] Liyang Sun and Sarah Abraham  “Estimating Dynamic Treatment Effects in Event Studies with Heterogeneous Treatment Effects”  In _Journal of Econometrics_ 225.2, 2021, pp. 175–199  DOI: 10.1016/j.jeconom.2020.09.006
  * [63] Max Tabord-Meehan  “Stratification Trees for Adaptive Randomization in Randomized Controlled Trials”, 2021
  * [64] Aad van der Vaart  “Asymptotic Statistics”  New York, NY: Cambridge University Press, 1998  DOI: 10.1017/CBO9780511802256
  * [65] Aad van der Vaart and Jon A. Wellner  “Prohorov and Continuous Mapping Theorems in the Hoffmann-Jørgensen Weak Convergence Theory, with Applications to Convolution and Asymptotic Minimax Theorems” Unpublished manuscript, Unversity of Seattle, 1989
  * [66] Aad van der Vaart and Jon A. Wellner  “Weak Convergence and Empirical Processes”  New York, NY: Springer, 1996  DOI: 10.1007/978-1-4757-2545-2
  * [67] Jeffrey M. Wooldridge  “Two-Way Fixed Effects, the Two-Way Mundlak Regression, and Difference-in-Differences Estimators”, 2021  DOI: 10.2139/ssrn.3906345

##  Appendix A Proofs

###  A.1 Proof of Proposition 1

We prove a generalization of the Proposition 1 which allows any vector of
treatments $X_{i}$ (which may not be binary or mutually exclusive). We
continue to consider the partially linear model in eq. 9, and maintain
Assumption 2, as well as conditional mean-independence of the potential
outcomes $E[Y_{i}(x)\mid X_{i},W_{i}]=E[Y_{i}(x)\mid W_{i}]$, which extends
Assumption 1. We also assume that the potential outcomes $Y_{i}(x)$ are linear
in $x$, conditional on $W_{i}$:

$E[Y_{i}(x)\mid W_{i}=w]=E[Y_{i}(0)\mid W_{i}=w]+x^{\prime}\tau(w),$

for some function $\tau$. This condition holds trivially in the main-text
discussion of mutually exclusive binary treatments. More generally,
$\tau_{k}(w)$ corresponds to the conditional average effect of increasing
$X_{ik}$ by one unit among observations with $W_{i}=w$. Although this
assumption is not essential, it considerably simplifies the derivations. We
continue to define $\tau=E[\tau(W_{i})]$ as the average vector of per-unit
effects.

We now prove that under these assumptions $\beta_{k}$ is given by the
expression in eq. 16. We further prove that $E[\lambda_{kk}(W_{i})]=1$ and
$E[\lambda_{k\ell}(W_{i})]=0$ for $\ell\neq k$ in general, and that
$\lambda_{kk}(W_{i})\geq 0$ in the case of mutually exclusive treatment
indicators.

First note that by iterated expectations and conditional mean-independence,
$E[\accentset{\approx}{X}_{ik}Y_{i}]=E[E[\accentset{\approx}{X}_{ik}Y_{i}\mid
X_{i},W_{i}]]=E[\accentset{\approx}{X}_{ik}E[Y_{i}(0)\mid
W_{i}]]+E[\accentset{\approx}{X}_{ik}X_{i}^{\prime}\tau(W_{i})]$. By
definition of projection, $E[\tilde{X}_{i}g(W_{i})]=0$ for all
$g\in\mathcal{G}$ [64, Theorem 11.1]; thus if eq. 14 holds
$E[\accentset{\approx}{X}_{ik}E[Y_{i}(0)\mid W_{i}]]=0$. Similarly, under eq.
13, $E[\accentset{\approx}{X}_{ik}\mid W_{i}]=0$, so by iterated expectations,
$E[\accentset{\approx}{X}_{ik}E[Y_{i}(0)\mid
W_{i}]]=E[E[\accentset{\approx}{X}_{ik}\mid W_{i}]E[Y_{i}(0)\mid W_{i}]]=0$.
Thus,

$\beta_{k}=\frac{E[\accentset{\approx}{X}_{ik}X_{i}^{\prime}\tau(W_{i})]}{E[\accentset{\approx}{X}_{ik}^{2}]}=\frac{E[\accentset{\approx}{X}_{ik}X_{ik}\tau_{k}(W_{i})]}{E[\accentset{\approx}{X}_{ik}^{2}]}+\frac{\sum_{\ell\neq
k}E[\accentset{\approx}{X}_{ik}X_{i\ell}\tau_{\ell}(W_{i})]}{E[\accentset{\approx}{X}_{ik}^{2}]}.$

This proves eq. 16.

To show that $E[\lambda_{kk}(W_{i})]=1$ and $E[\lambda_{k\ell}(W_{i})]=0$ for
$\ell\neq k$ in general, note that

$E[\lambda_{kk}(W_{i})]=\frac{E[\accentset{\approx}{X}_{ik}X_{ik}]}{E[\accentset{\approx}{X}_{ik}^{2}]}=1,$

since $\accentset{\approx}{X}_{i,k}$ is a residual from projecting $X_{ik}$
onto the space spanned by functions of the form
$\tilde{g}(W_{i})+X_{i,-k}^{\prime}\tilde{\beta}_{-k}$, so that
$E[\accentset{\approx}{X}_{ik}X_{ik}]=E[\accentset{\approx}{X}_{ik}^{2}]$.
Furthermore, $\accentset{\approx}{X}_{i,k}$ must also be orthogonal to
$X_{i,-k}$ by definition of projection, so that
$E[\lambda_{k\ell}(W_{i})]=E[\accentset{\approx}{X}_{ik}X_{i\ell}]/E[\accentset{\approx}{X}_{ik}^{2}]=0$.

Finally, we show that $\lambda_{kk}(W_{i})\geq 0$ if eq. 13 holds and $X_{i}$
consists of mutually exclusive indicators. To that end, observe that
$\lambda_{k\ell}(W_{i})$ is given by the $(k,\ell)$ element of

$\Lambda(W_{i})=E[\tilde{X}_{i}\tilde{X}_{i}^{\prime}]^{-1}E[\tilde{X}_{i}X_{i}^{\prime}\mid
W_{i}]$

If Equation 13 holds, then we can write this as
$\Lambda(W_{i})=E[v(W_{i})]^{-1}v(W_{i})$ where
$v(W_{i})=E[\tilde{X}_{i}\tilde{X}_{i}^{\prime}\mid W_{i}]$. If $X$ is a
vector of mutually exclusive indicators, then
$v(W_{i})=\operatorname{diag}(p(W_{i}))-p(W_{i})p(W_{i})^{\prime}$. Let
$v_{-k}(W_{i})$ denote the submatrix with the $k$th row and column removed,
and for any vector $a$, let $a_{-k}$ denote subvector with the $k$th column
removed. Then by the block matrix inverse formula,

$\lambda_{kk}(W_{i})=\frac{p_{k}(W_{i})(1-p_{k}(W_{i}))-E[p_{k}(W_{i})p(W_{i})_{-k}^{\prime}]E[v_{-k}(W_{i})]^{-1}p(W_{i})_{-k}p_{k}(W_{i})}{E[p_{k}(W_{i})(1-p_{k}(W_{i}))]-E[p_{k}(W_{i})p_{-k}(W_{i})^{\prime}]E[v_{-k}(W_{i})]^{-1}E[p_{k}(W_{i})p_{-k}(W_{i})]}$

Note $p_{0}(W_{i})=1-\sum_{k=1}^{K}p_{k}(W_{i})$ and
$p_{k}(W_{i})p(W_{i})_{-k}=v_{-k}(W_{i})\iota-p_{0}(W_{i})p(W_{i})_{-k}$,
where $\iota$ denotes a $(K-1)$-vector of ones. Thus the numerator can be
written as

$p_{k}(W_{i})(1-p_{k}(W_{i}))-\iota^{\prime}p(W_{i})_{-k}p_{k}(W_{i})\\\
+E[p_{0}(W_{i})p(W_{i})_{-k}^{\prime}]E[v_{-k}(W_{i})]^{-1}p(W_{i})_{-k}p_{k}(W_{i})\\\
=p_{k}(W_{i})p_{0}(W_{i})+E[p_{0}(W_{i})p(W_{i})_{-k}^{\prime}]E[v_{-k}(W_{i})]^{-1}p(W_{i})_{-k}p_{k}(W_{i}).$

The eigenvalues of $E[v_{-k}(W_{i})]$ are positive because it is a covariance
matrix. Furthermore, since the off-diagonal elements of $E[v(W_{i})]$ are
negative, the off-diagonal elements of $E[v_{-k}(W_{i})]$ are also negative.
It therefore follows that $E[v_{-k}(W_{i})]$ is an $M$-matrix [12, property
$D_{16}$, p. 135]. Hence, all elements of $E[v_{-k}(W_{i})]^{-1}$ are positive
[12, property $N_{38}$, p. 137]. Thus, both summands in the above expression
are positive, so that $\lambda_{kk}(W_{i})\geq 0$.

###  A.2 Proof of Proposition 2

The parameter of interest $\theta_{\lambda,c}$ depends on the realizations of
the controls. We therefore derive the semiparametric efficiency bound
conditional on the controls; i.e. we show that eq. 18 is almost surely the
variance bound for estimators that are regular conditional on the controls.
Relative to the earlier results in [32] and [33], we need account for the fact
that the data are no longer i.i.d. once we condition on the controls.

To that end, we use the notion of semiparametric efficiency based on the
convolution theorem of [65, Theorem 2.1] [[, see also]Chapter 3.11]vVWe96. We
first review the result for convenience. Consider a model
$\\{P_{n,\theta}\colon\theta\in\Theta\\}$ parametrized by (a possibly
infinite-dimensional) parameter $\theta$. Let $\dot{\mathcal{P}}$ denote a
tangent space, a linear subspace of some Hilbert space with an inner product
$\langle\cdot,\cdot\rangle$. Suppose that the model is locally asymptotically
normal (LAN) at $\theta$ relative to a tangent space $\dot{\mathcal{P}}$: for
each $g\in\dot{\mathcal{P}}$, there exists a sequence $\theta_{n}(g)$ such
that the likelihood ratios are asymptotically quadratic,
$dP_{n,\theta_{n}(g)}/dP_{n,\theta}=\Delta_{n,g}-\langle
g,g\rangle/2+o_{P_{n,\theta}}(1)$, where
$(\Delta_{n,g})_{g\in\dot{\mathcal{P}}}$ converges under $P_{n,\theta}$ to a
Gaussian process with covariance kernel $\langle g_{1},g_{2}\rangle$. Suppose
also that the parameter $\beta_{n}(P_{n,\theta})$ is differentiable: for each
$g$,
$\sqrt{n}(\beta_{n}(P_{n,\theta_{n}(g)})-\beta_{n}(P_{n,\theta}))\to\langle\psi,g\rangle$
for some $\psi$ that lies in the completion of $\dot{\mathcal{P}}$. Then the
semiparametric efficiency bound is given by $\langle\psi,\psi\rangle$: the
asymptotic distribution of any regular estimator of this parameter, based on a
sample $\mathbf{S}_{n}\sim P_{n,\theta}$, is given by the convolution of a
random variable $Z\sim\mathcal{N}(0,\langle\psi,\psi\rangle)$ and some other
random variable $U$ that is independent of $Z$.

To apply this result in our setting, we proceed in three steps. First, we
define the tangent space and the probability-one set over which we will prove
the efficiency bound. Next, we verify that the model is LAN. Finally, we
verify differentiability and derive the efficient influence function $\psi$.

#### Step 1

By the conditional independence assumption in eq. 12, we can write the density
of the vector $(Y_{i}(0),\dotsc,Y_{i}(K),D_{i})$ (with respect to some
$\sigma$-finite measure) conditional on $W_{i}=w$ as $f(y_{0},\dotsc,y_{K}\mid
w)\cdot\prod_{k=0}^{K}p_{k}(w)^{\operatorname{\mathbbm{1}}\\{d=k\\}}$, where
$f$ denotes the conditional density of the potential outcomes, conditional on
the controls. The density of the observed data
$\mathbf{S}_{N}=\\{(Y_{i},D_{i})\\}_{i=1}^{N}$ conditional on
$(W_{1},\dotsc,W_{N})=(w_{1},\dotsc,w_{N})$ is given by
$\prod_{i=1}^{N}\prod_{k=0}^{K}(f_{k}(y_{i}\mid
w_{i})p_{k}(w_{i}))^{\operatorname{\mathbbm{1}}\\{d_{i}=k\\}}$, where
$f_{k}(y\mid w)=\int f(y_{k},y_{-k}\mid w)dy_{-k}$.

Since the propensity scores are known, the model is parametrized by
$\theta=f$. Consider one-dimensional submodels of the form $f_{k}(y\mid
w;t)=f_{k}(y\mid w)(1+ts_{k}(y\mid w))$, where the function $s_{k}$ is bounded
and satisfies $\int s_{k}(y\mid w)f_{k}(y\mid w)dy=0$ for all
$w\in\mathcal{W}$ with $\mathcal{W}$ denoting the support of $W_{i}$. For
small enough $t$, we have $f_{k}(y\mid w;t)\geq 0$ by boundedness of $s_{k}$;
hence $f_{k}(y\mid w;t)$ is a well-defined density for $t$ small enough. The
joint log-likelihood, conditional on the controls, is given by

$\sum_{i=1}^{N}\sum_{k=0}^{K}\operatorname{\mathbbm{1}}\\{D_{i}=k\\}(\log
f_{k}(Y_{i}\mid w_{i};t)+\log p_{k}(w_{i})).$

The score at $t=0$ is $\sum_{i=1}^{N}s(Y_{i},D_{i}\mid w_{i})$, with
$s(Y_{i},D_{i}\mid
w_{i})=\sum_{k=0}^{K}\operatorname{\mathbbm{1}}\\{D_{i}=k\\}s_{k}(Y_{i}\mid
w_{i})$.

This result suggests defining the tangent space to consist of functions
$s(y,d\mid w)=\sum_{k=0}^{K}\operatorname{\mathbbm{1}}\\{d=k\\}s_{k}(y\mid
W_{i}=w)$, such that $s_{k}$ is bounded and satisfies $\int s_{k}(y\mid
w)f_{k}(y\mid w)dy=0$ for all $w\in\mathcal{W}$. Define the inner product on
this space by $\langle s_{1},s_{2}\rangle=E[s_{1}(Y_{i},D_{i}\mid
W_{i})s_{2}(Y_{i},D_{i}\mid W_{i})]$. Note this is a marginal (rather than a
conditional) expectation, over the unconditional distribution
$(Y_{i},D_{i},W_{i})$ of the observed data.

We will prove the efficiency bound on the event $\mathcal{E}$ that (i)
$\frac{1}{N}\sum_{i=1}^{N}E[s(Y_{i},D_{i}\mid W_{i})^{2}\mid W_{i}]\to
E[s(Y,D_{i}\mid W_{i})^{2}]$, (ii) $\frac{1}{N}\sum_{i=1}^{N}\lambda(W_{i})\to
E[\lambda(W_{i})]$, and (iii)
$\frac{1}{N}\sum_{i=1}^{N}\lambda(W_{i})\sum_{k=0}^{K}c_{k}\cdot
E[Y_{i}(k)s_{k}(Y_{i}(k)\mid W_{i})\mid
W_{i}]\to\sum_{k=0}^{K}c_{k}E[\lambda(W_{i})Y_{i}(k)s_{k}(Y_{i}(k)\mid
W_{i})]$. By assumptions of the proposition, these are all averages of
functions of $W_{i}$ with finite absolute moments. Hence, by the law of large
numbers, $\mathcal{E}$ is a probability one set.

#### Step 2